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 = a.
function_space(0)->dofmap();
105 std::shared_ptr<const fem::DofMap> dofmap1 = a.
function_space(1)->dofmap();
113 throw std::runtime_error(
"Unset constant in Form");
114 const Eigen::Array<T, Eigen::Dynamic, 1> constants =
pack_constants(a);
117 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> coeffs
122 if (needs_permutation_data)
123 mesh->topology_mutable().create_entity_permutations();
124 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info
125 = needs_permutation_data
126 ? mesh->topology().get_cell_permutation_info()
127 : Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>(num_cells);
129 for (
int i = 0; i < integrals.
num_integrals(IntegralType::cell); ++i)
132 const std::vector<std::int32_t>& active_cells
134 impl::assemble_cells<T>(mat_set_values, mesh->geometry(), active_cells,
135 dofs0, dofs1, bc0, bc1, fn, coeffs, constants,
139 if (integrals.
num_integrals(IntegralType::exterior_facet) > 0
140 or integrals.
num_integrals(IntegralType::interior_facet) > 0)
143 mesh->topology_mutable().create_entities(tdim - 1);
144 mesh->topology_mutable().create_connectivity(tdim - 1, tdim);
146 const int facets_per_cell
148 const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>& perms
149 = needs_permutation_data
150 ? mesh->topology().get_facet_permutations()
151 : Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>(
152 facets_per_cell, num_cells);
154 for (
int i = 0; i < integrals.
num_integrals(IntegralType::exterior_facet);
159 const std::vector<std::int32_t>& active_facets
161 impl::assemble_exterior_facets<T>(mat_set_values, *mesh, active_facets,
162 dofs0, dofs1, bc0, bc1, fn, coeffs,
163 constants, cell_info, perms);
166 const std::vector<int> c_offsets = a.
coefficients().offsets();
167 for (
int i = 0; i < integrals.
num_integrals(IntegralType::interior_facet);
172 const std::vector<std::int32_t>& active_facets
174 impl::assemble_interior_facets<T>(
175 mat_set_values, *mesh, active_facets, *dofmap0, *dofmap1, bc0, bc1,
176 fn, coeffs, c_offsets, constants, cell_info, perms);
181 template <
typename T>
183 const std::function<
int(std::int32_t,
const std::int32_t*, std::int32_t,
184 const std::int32_t*,
const T*)>& mat_set,
186 const std::vector<std::int32_t>& active_cells,
189 const std::vector<bool>& bc0,
const std::vector<bool>& bc1,
190 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
191 const std::uint8_t*,
const std::uint32_t)>& kernel,
192 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
194 const Eigen::Array<T, Eigen::Dynamic, 1>& constants,
195 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info)
197 const int gdim = geometry.
dim();
203 const int num_dofs_g = x_dofmap.
num_links(0);
204 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
208 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
209 coordinate_dofs(num_dofs_g, gdim);
210 const int num_dofs0 = dofmap0.
links(0).size();
211 const int num_dofs1 = dofmap1.
links(0).size();
212 Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> Ae(
213 num_dofs0, num_dofs1);
216 for (std::int32_t c : active_cells)
219 auto x_dofs = x_dofmap.
links(c);
220 for (
int i = 0; i < x_dofs.rows(); ++i)
221 for (
int j = 0; j < gdim; ++j)
222 coordinate_dofs(i, j) = x_g(x_dofs[i], j);
225 std::fill(Ae.data(), Ae.data() + num_dofs0 * num_dofs1, 0);
226 kernel(Ae.data(), coeffs.row(c).data(), constants.data(),
227 coordinate_dofs.data(),
nullptr,
nullptr, cell_info[c]);
230 auto dofs0 = dofmap0.
links(c);
231 auto dofs1 = dofmap1.
links(c);
234 for (Eigen::Index i = 0; i < Ae.rows(); ++i)
242 for (Eigen::Index j = 0; j < Ae.cols(); ++j)
249 mat_set(dofs0.size(), dofs0.data(), dofs1.size(), dofs1.data(), Ae.data());
253 template <
typename T>
254 void assemble_exterior_facets(
255 const std::function<
int(std::int32_t,
const std::int32_t*, std::int32_t,
256 const std::int32_t*,
const T*)>& mat_set_values,
257 const mesh::Mesh& mesh,
const std::vector<std::int32_t>& active_facets,
260 const std::vector<bool>& bc0,
const std::vector<bool>& bc1,
261 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
262 const std::uint8_t*,
const std::uint32_t)>& kernel,
263 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
265 const Eigen::Array<T, Eigen::Dynamic, 1> constants,
266 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info,
267 const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>& perms)
276 const int num_dofs_g = x_dofmap.
num_links(0);
277 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
281 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
282 coordinate_dofs(num_dofs_g, gdim);
283 const int num_dofs0 = dofmap0.
links(0).size();
284 const int num_dofs1 = dofmap1.
links(0).size();
285 Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> Ae(
286 num_dofs0, num_dofs1);
293 for (std::int32_t f : active_facets)
295 auto cells = f_to_c->links(f);
296 assert(cells.rows() == 1);
299 auto facets = c_to_f->links(cells[0]);
300 const auto* it = std::find(facets.data(), facets.data() + facets.rows(), f);
301 assert(it != (facets.data() + facets.rows()));
302 const int local_facet = std::distance(facets.data(), it);
305 auto x_dofs = x_dofmap.
links(cells[0]);
306 for (
int i = 0; i < num_dofs_g; ++i)
307 for (
int j = 0; j < gdim; ++j)
308 coordinate_dofs(i, j) = x_g(x_dofs[i], j);
311 std::fill(Ae.data(), Ae.data() + num_dofs0 * num_dofs1, 0);
312 kernel(Ae.data(), coeffs.row(cells[0]).data(), constants.data(),
313 coordinate_dofs.data(), &local_facet, &perms(local_facet, cells[0]),
314 cell_info[cells[0]]);
317 auto dmap0 = dofmap0.
links(cells[0]);
318 auto dmap1 = dofmap1.
links(cells[0]);
321 for (Eigen::Index i = 0; i < Ae.rows(); ++i)
329 for (Eigen::Index j = 0; j < Ae.cols(); ++j)
336 mat_set_values(dmap0.size(), dmap0.data(), dmap1.size(), dmap1.data(),
341 template <
typename T>
342 void assemble_interior_facets(
343 const std::function<
int(std::int32_t,
const std::int32_t*, std::int32_t,
344 const std::int32_t*,
const T*)>& mat_set_values,
345 const mesh::Mesh& mesh,
const std::vector<std::int32_t>& active_facets,
346 const DofMap& dofmap0,
const DofMap& dofmap1,
const std::vector<bool>& bc0,
347 const std::vector<bool>& bc1,
348 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
349 const std::uint8_t*,
const std::uint32_t)>& fn,
350 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
352 const std::vector<int>& offsets,
353 const Eigen::Array<T, Eigen::Dynamic, 1>& constants,
354 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info,
355 const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>& perms)
364 const int num_dofs_g = x_dofmap.
num_links(0);
366 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
370 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
371 coordinate_dofs(2 * num_dofs_g, gdim);
372 Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> Ae;
373 Eigen::Array<T, Eigen::Dynamic, 1> coeff_array(2 * offsets.back());
374 assert(offsets.back() == coeffs.cols());
377 Eigen::Array<std::int32_t, Eigen::Dynamic, 1> dmapjoint0, dmapjoint1;
384 for (std::int32_t facet_index : active_facets)
387 auto cells = c->links(facet_index);
388 assert(cells.rows() == 2);
391 auto facets0 = c_to_f->links(cells[0]);
392 const auto* it0 = std::find(facets0.data(), facets0.data() + facets0.rows(),
394 assert(it0 != (facets0.data() + facets0.rows()));
395 const int local_facet0 = std::distance(facets0.data(), it0);
396 auto facets1 = c_to_f->links(cells[1]);
397 const auto* it1 = std::find(facets1.data(), facets1.data() + facets1.rows(),
399 assert(it1 != (facets1.data() + facets1.rows()));
400 const int local_facet1 = std::distance(facets1.data(), it1);
402 const std::array local_facet{local_facet0, local_facet1};
405 auto x_dofs0 = x_dofmap.
links(cells[0]);
406 auto x_dofs1 = x_dofmap.
links(cells[1]);
407 for (
int i = 0; i < num_dofs_g; ++i)
409 for (
int j = 0; j < gdim; ++j)
411 coordinate_dofs(i, j) = x_g(x_dofs0[i], j);
412 coordinate_dofs(i + num_dofs_g, j) = x_g(x_dofs1[i], j);
417 auto dmap0_cell0 = dofmap0.
cell_dofs(cells[0]);
418 auto dmap0_cell1 = dofmap0.
cell_dofs(cells[1]);
419 dmapjoint0.resize(dmap0_cell0.size() + dmap0_cell1.size());
420 dmapjoint0.head(dmap0_cell0.size()) = dmap0_cell0;
421 dmapjoint0.tail(dmap0_cell1.size()) = dmap0_cell1;
423 auto dmap1_cell0 = dofmap1.
cell_dofs(cells[0]);
424 auto dmap1_cell1 = dofmap1.
cell_dofs(cells[1]);
425 dmapjoint1.resize(dmap1_cell0.size() + dmap1_cell1.size());
426 dmapjoint1.head(dmap1_cell0.size()) = dmap1_cell0;
427 dmapjoint1.tail(dmap1_cell1.size()) = dmap1_cell1;
431 auto coeff_cell0 = coeffs.row(cells[0]);
432 auto coeff_cell1 = coeffs.row(cells[1]);
435 for (std::size_t i = 0; i < offsets.size() - 1; ++i)
438 const int num_entries = offsets[i + 1] - offsets[i];
439 coeff_array.segment(2 * offsets[i], num_entries)
440 = coeff_cell0.segment(offsets[i], num_entries);
441 coeff_array.segment(offsets[i + 1] + offsets[i], num_entries)
442 = coeff_cell1.segment(offsets[i], num_entries);
446 Ae.setZero(dmapjoint0.size(), dmapjoint1.size());
447 const std::array perm{perms(local_facet[0], cells[0]),
448 perms(local_facet[1], cells[1])};
449 fn(Ae.data(), coeff_array.data(), constants.data(), coordinate_dofs.data(),
450 local_facet.data(), perm.data(), cell_info[cells[0]]);
455 for (Eigen::Index i = 0; i < dmapjoint0.size(); ++i)
457 if (bc0[dmapjoint0[i]])
463 for (Eigen::Index j = 0; j < dmapjoint1.size(); ++j)
465 if (bc1[dmapjoint1[j]])
470 mat_set_values(dmapjoint0.size(), dmapjoint0.data(), dmapjoint1.size(),
471 dmapjoint1.data(), Ae.data());