13 #include <dolfinx/fem/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>
23 namespace dolfinx::fem::impl
34 const std::function<
int(std::int32_t,
const std::int32_t*, std::int32_t,
35 const std::int32_t*,
const T*)>& mat_set_values,
36 const Form<T>& a,
const std::vector<bool>& bc0,
37 const std::vector<bool>& bc1);
42 const std::function<
int(std::int32_t,
const std::int32_t*, std::int32_t,
43 const std::int32_t*,
const T*)>& mat_set_values,
45 const std::vector<std::int32_t>& active_cells,
48 const std::vector<bool>& bc0,
const std::vector<bool>& bc1,
49 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
50 const std::uint8_t*,
const std::uint32_t)>& kernel,
51 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
53 const std::vector<T>& constants,
54 const std::vector<std::uint32_t>& cell_info);
58 void assemble_exterior_facets(
59 const std::function<
int(std::int32_t,
const std::int32_t*, std::int32_t,
60 const std::int32_t*,
const T*)>& mat_set_values,
61 const mesh::Mesh& mesh,
const std::vector<std::int32_t>& active_facets,
64 const std::vector<bool>& bc0,
const std::vector<bool>& bc1,
65 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
66 const std::uint8_t*,
const std::uint32_t)>& fn,
67 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
69 const std::vector<T>& constants,
70 const std::vector<std::uint32_t>& cell_info,
71 const std::vector<std::uint8_t>& perms);
75 void assemble_interior_facets(
76 const std::function<
int(std::int32_t,
const std::int32_t*, std::int32_t,
77 const std::int32_t*,
const T*)>& mat_set_values,
78 const mesh::Mesh& mesh,
const std::vector<std::int32_t>& active_facets,
79 const DofMap& dofmap0,
int bs0,
const DofMap& dofmap1,
int bs1,
80 const std::vector<bool>& bc0,
const std::vector<bool>& bc1,
81 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
82 const std::uint8_t*,
const std::uint32_t)>& kernel,
83 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
85 const std::vector<int>& offsets,
const std::vector<T>& constants,
86 const std::vector<std::uint32_t>& cell_info,
87 const std::vector<std::uint8_t>& 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
111 const int bs0 = dofmap0->bs();
113 const int bs1 = dofmap1->bs();
119 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> coeffs
123 if (needs_permutation_data)
124 mesh->topology_mutable().create_entity_permutations();
125 const std::vector<std::uint32_t>& cell_info
126 = needs_permutation_data ? mesh->topology().get_cell_permutation_info()
127 : std::vector<std::uint32_t>(num_cells);
131 const auto& fn = a.
kernel(IntegralType::cell, i);
132 const std::vector<std::int32_t>& active_cells
133 = a.
domains(IntegralType::cell, i);
134 impl::assemble_cells<T>(mat_set_values, mesh->geometry(), active_cells,
135 dofs0, bs0, dofs1, bs1, bc0, bc1, fn, coeffs,
136 constants, cell_info);
142 mesh->topology_mutable().create_connectivity(tdim - 1, tdim);
143 mesh->topology_mutable().create_entity_permutations();
145 const std::vector<std::uint8_t>& perms
146 = mesh->topology().get_facet_permutations();
148 for (
int i : a.
integral_ids(IntegralType::exterior_facet))
150 const auto& fn = a.
kernel(IntegralType::exterior_facet, i);
151 const std::vector<std::int32_t>& active_facets
152 = a.
domains(IntegralType::exterior_facet, i);
153 impl::assemble_exterior_facets<T>(mat_set_values, *mesh, active_facets,
154 dofs0, bs0, dofs1, bs1, bc0, bc1, fn,
155 coeffs, constants, cell_info, perms);
159 for (
int i : a.
integral_ids(IntegralType::interior_facet))
161 const auto& fn = a.
kernel(IntegralType::interior_facet, i);
162 const std::vector<std::int32_t>& active_facets
163 = a.
domains(IntegralType::interior_facet, i);
164 impl::assemble_interior_facets<T>(
165 mat_set_values, *mesh, active_facets, *dofmap0, bs0, *dofmap1, bs1,
166 bc0, bc1, fn, coeffs, c_offsets, constants, cell_info, perms);
171 template <
typename T>
173 const std::function<
int(std::int32_t,
const std::int32_t*, std::int32_t,
174 const std::int32_t*,
const T*)>& mat_set,
176 const std::vector<std::int32_t>& active_cells,
179 const std::vector<bool>& bc0,
const std::vector<bool>& bc1,
180 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
181 const std::uint8_t*,
const std::uint32_t)>& kernel,
182 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
184 const std::vector<T>& constants,
185 const std::vector<std::uint32_t>& cell_info)
187 const int gdim = geometry.
dim();
193 const int num_dofs_g = x_dofmap.
num_links(0);
194 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
198 const int num_dofs0 = dofmap0.
links(0).size();
199 const int num_dofs1 = dofmap1.
links(0).size();
200 const int ndim0 = bs0 * num_dofs0;
201 const int ndim1 = bs1 * num_dofs1;
202 std::vector<T> Ae(ndim0 * ndim1);
203 std::vector<double> coordinate_dofs(num_dofs_g * gdim);
204 for (std::int32_t c : active_cells)
207 auto x_dofs = x_dofmap.
links(c);
208 for (std::size_t i = 0; i < x_dofs.size(); ++i)
210 std::copy_n(x_g.row(x_dofs[i]).data(), gdim,
211 std::next(coordinate_dofs.begin(), i * gdim));
215 std::fill(Ae.begin(), Ae.end(), 0);
216 kernel(Ae.data(), coeffs.row(c).data(), constants.data(),
217 coordinate_dofs.data(),
nullptr,
nullptr, cell_info[c]);
220 auto dofs0 = dofmap0.
links(c);
221 auto dofs1 = dofmap1.
links(c);
224 for (
int i = 0; i < num_dofs0; ++i)
226 for (
int k = 0; k < bs0; ++k)
228 if (bc0[bs0 * dofs0[i] + k])
231 const int row = bs0 * i + k;
232 std::fill_n(std::next(Ae.begin(), ndim1 * row), ndim1, 0.0);
240 for (
int j = 0; j < num_dofs1; ++j)
242 for (
int k = 0; k < bs1; ++k)
244 if (bc1[bs1 * dofs1[j] + k])
247 const int col = bs1 * j + k;
248 for (
int row = 0; row < ndim0; ++row)
249 Ae[row * ndim1 + col] = 0.0;
255 mat_set(dofs0.size(), dofs0.data(), dofs1.size(), dofs1.data(), Ae.data());
259 template <
typename T>
260 void assemble_exterior_facets(
261 const std::function<
int(std::int32_t,
const std::int32_t*, std::int32_t,
262 const std::int32_t*,
const T*)>& mat_set_values,
263 const mesh::Mesh& mesh,
const std::vector<std::int32_t>& active_facets,
266 const std::vector<bool>& bc0,
const std::vector<bool>& bc1,
267 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
268 const std::uint8_t*,
const std::uint32_t)>& kernel,
269 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
271 const std::vector<T>& constants,
272 const std::vector<std::uint32_t>& cell_info,
273 const std::vector<std::uint8_t>& perms)
282 const int num_dofs_g = x_dofmap.
num_links(0);
283 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
287 std::vector<double> coordinate_dofs(num_dofs_g * gdim);
288 const int num_dofs0 = dofmap0.
links(0).size();
289 const int num_dofs1 = dofmap1.
links(0).size();
290 const int ndim0 = bs0 * num_dofs0;
291 const int ndim1 = bs1 * num_dofs1;
292 std::vector<T> Ae(ndim0 * ndim1);
299 for (std::int32_t f : active_facets)
301 auto cells = f_to_c->links(f);
302 assert(cells.size() == 1);
305 auto facets = c_to_f->links(cells[0]);
306 auto it = std::find(facets.begin(), facets.end(), f);
307 assert(it != facets.end());
308 const int local_facet = std::distance(facets.begin(), it);
311 auto x_dofs = x_dofmap.
links(cells[0]);
312 for (std::size_t i = 0; i < x_dofs.size(); ++i)
314 std::copy_n(x_g.row(x_dofs[i]).data(), gdim,
315 std::next(coordinate_dofs.begin(), i * gdim));
319 std::fill(Ae.begin(), Ae.end(), 0);
320 kernel(Ae.data(), coeffs.row(cells[0]).data(), constants.data(),
321 coordinate_dofs.data(), &local_facet,
322 &perms[cells[0] * facets.size() + local_facet], cell_info[cells[0]]);
325 auto dofs0 = dofmap0.
links(cells[0]);
326 auto dofs1 = dofmap1.
links(cells[0]);
329 for (
int i = 0; i < num_dofs0; ++i)
331 for (
int k = 0; k < bs0; ++k)
333 if (bc0[bs0 * dofs0[i] + k])
336 const int row = bs0 * i + k;
337 std::fill_n(std::next(Ae.begin(), ndim1 * row), ndim1, 0.0);
344 for (
int j = 0; j < num_dofs1; ++j)
346 for (
int k = 0; k < bs1; ++k)
348 if (bc1[bs1 * dofs1[j] + k])
351 const int col = bs1 * j + k;
352 for (
int row = 0; row < ndim0; ++row)
353 Ae[row * ndim1 + col] = 0.0;
359 mat_set_values(dofs0.size(), dofs0.data(), dofs1.size(), dofs1.data(),
364 template <
typename T>
365 void assemble_interior_facets(
366 const std::function<
int(std::int32_t,
const std::int32_t*, std::int32_t,
367 const std::int32_t*,
const T*)>& mat_set_values,
368 const mesh::Mesh& mesh,
const std::vector<std::int32_t>& active_facets,
369 const DofMap& dofmap0,
int bs0,
const DofMap& dofmap1,
int bs1,
370 const std::vector<bool>& bc0,
const std::vector<bool>& bc1,
371 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
372 const std::uint8_t*,
const std::uint32_t)>& fn,
373 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
375 const std::vector<int>& offsets,
const std::vector<T>& constants,
376 const std::vector<std::uint32_t>& cell_info,
377 const std::vector<std::uint8_t>& perms)
386 const int num_dofs_g = x_dofmap.
num_links(0);
388 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
392 std::vector<double> coordinate_dofs(2 * num_dofs_g * gdim);
393 std::vector<T> Ae, be;
394 std::vector<T> coeff_array(2 * offsets.back());
395 assert(offsets.back() == coeffs.cols());
398 std::vector<std::int32_t> dmapjoint0, dmapjoint1;
405 const int offset_g = gdim * num_dofs_g;
406 for (std::int32_t facet_index : active_facets)
409 auto cells = c->links(facet_index);
410 assert(cells.size() == 2);
413 auto facets0 = c_to_f->links(cells[0]);
414 const auto* it0 = std::find(facets0.begin(), facets0.end(), facet_index);
415 assert(it0 != facets0.end());
416 const int local_facet0 = std::distance(facets0.begin(), it0);
417 auto facets1 = c_to_f->links(cells[1]);
418 const auto* it1 = std::find(facets1.begin(), facets1.end(), facet_index);
419 assert(it1 != facets1.end());
420 const int local_facet1 = std::distance(facets1.begin(), it1);
422 const std::array local_facet{local_facet0, local_facet1};
425 auto x_dofs0 = x_dofmap.
links(cells[0]);
426 auto x_dofs1 = x_dofmap.
links(cells[1]);
427 for (
int i = 0; i < num_dofs_g; ++i)
429 for (
int j = 0; j < gdim; ++j)
431 coordinate_dofs[i * gdim + j] = x_g(x_dofs0[i], j);
432 coordinate_dofs[offset_g + i * gdim + j] = x_g(x_dofs1[i], j);
437 tcb::span<const std::int32_t> dmap0_cell0 = dofmap0.cell_dofs(cells[0]);
438 tcb::span<const std::int32_t> dmap0_cell1 = dofmap0.cell_dofs(cells[1]);
439 dmapjoint0.resize(dmap0_cell0.size() + dmap0_cell1.size());
440 std::copy(dmap0_cell0.begin(), dmap0_cell0.end(), dmapjoint0.begin());
441 std::copy(dmap0_cell1.begin(), dmap0_cell1.end(),
442 std::next(dmapjoint0.begin(), dmap0_cell0.size()));
444 tcb::span<const std::int32_t> dmap1_cell0 = dofmap1.cell_dofs(cells[0]);
445 tcb::span<const std::int32_t> dmap1_cell1 = dofmap1.cell_dofs(cells[1]);
446 dmapjoint1.resize(dmap1_cell0.size() + dmap1_cell1.size());
447 std::copy(dmap1_cell0.begin(), dmap1_cell0.end(), dmapjoint1.begin());
448 std::copy(dmap1_cell1.begin(), dmap1_cell1.end(),
449 std::next(dmapjoint1.begin(), dmap1_cell0.size()));
453 auto coeff_cell0 = coeffs.row(cells[0]);
454 auto coeff_cell1 = coeffs.row(cells[1]);
457 for (std::size_t i = 0; i < offsets.size() - 1; ++i)
460 const int num_entries = offsets[i + 1] - offsets[i];
461 std::copy_n(coeff_cell0.data() + offsets[i], num_entries,
462 std::next(coeff_array.begin(), 2 * offsets[i]));
463 std::copy_n(coeff_cell1.data() + offsets[i], num_entries,
464 std::next(coeff_array.begin(), offsets[i + 1] + offsets[i]));
467 const int num_rows = bs0 * dmapjoint0.size();
468 const int num_cols = bs1 * dmapjoint1.size();
471 Ae.resize(num_rows * num_cols);
472 std::fill(Ae.begin(), Ae.end(), 0);
474 const int facets_per_cell = facets0.size();
475 const std::array perm{perms[cells[0] * facets_per_cell + local_facet[0]],
476 perms[cells[1] * facets_per_cell + local_facet[1]]};
477 fn(Ae.data(), coeff_array.data(), constants.data(), coordinate_dofs.data(),
478 local_facet.data(), perm.data(), cell_info[cells[0]]);
483 for (std::size_t i = 0; i < dmapjoint0.size(); ++i)
485 for (
int k = 0; k < bs0; ++k)
487 if (bc0[bs0 * dmapjoint0[i] + k])
490 std::fill_n(std::next(Ae.begin(), num_cols * (bs0 * i + k)),
498 for (std::size_t j = 0; j < dmapjoint1.size(); ++j)
500 for (
int k = 0; k < bs1; ++k)
502 if (bc1[bs1 * dmapjoint1[j] + k])
505 for (
int m = 0; m < num_rows; ++m)
506 Ae[m * num_cols + bs1 * j + k] = 0.0;
512 mat_set_values(dmapjoint0.size(), dmapjoint0.data(), dmapjoint1.size(),
513 dmapjoint1.data(), Ae.data());
Degree-of-freedom map.
Definition: DofMap.h:65
This class provides a static adjacency list data structure. It is commonly used to store directed gra...
Definition: AdjacencyList.h:46
tcb::span< T > links(int node)
Get the links (edges) for given node.
Definition: AdjacencyList.h:128
int num_links(int node) const
Number of connections for given node.
Definition: AdjacencyList.h:118
Geometry stores the geometry imposed on a mesh.
Definition: Geometry.h:39
const graph::AdjacencyList< std::int32_t > & dofmap() const
DOF map.
Definition: Geometry.cpp:22
Eigen::Array< double, Eigen::Dynamic, 3, Eigen::RowMajor > & x()
Geometry degrees-of-freedom.
Definition: Geometry.cpp:32
int dim() const
Return Euclidean dimension of coordinate system.
Definition: Geometry.cpp:20
A Mesh consists of a set of connected and numbered mesh topological entities, and geometry data.
Definition: Mesh.h:57
Geometry & geometry()
Get mesh geometry.
Definition: Mesh.cpp:134
Topology & topology()
Get mesh topology.
Definition: Mesh.cpp:128
std::shared_ptr< const graph::AdjacencyList< std::int32_t > > connectivity(int d0, int d1) const
Return connectivity from entities of dimension d0 to entities of dimension d1.
Definition: Topology.cpp:261
int dim() const
Return topological dimension.
Definition: Topology.cpp:160
std::vector< typename U::scalar_type > pack_constants(const U &u)
Pack constants of u of generic type U ready for assembly.
Definition: utils.h:455
void assemble_matrix(const std::function< int(std::int32_t, const std::int32_t *, std::int32_t, const std::int32_t *, const T *)> &mat_add, const Form< T > &a, const std::vector< std::shared_ptr< const DirichletBC< T >>> &bcs)
Assemble bilinear form into a matrix.
Definition: assembler.h:86
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