12 #include <dolfinx/fem/FunctionSpace.h>
13 #include <dolfinx/graph/AdjacencyList.h>
14 #include <dolfinx/la/utils.h>
15 #include <dolfinx/mesh/Geometry.h>
16 #include <dolfinx/mesh/Mesh.h>
17 #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 array2d<T>& coeffs,
const std::vector<T>& constants,
51 const std::vector<std::uint32_t>& cell_info);
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 array2d<T>& coeffs,
const std::vector<T>& constants,
65 const std::vector<std::uint32_t>& cell_info,
66 const std::vector<std::uint8_t>& perms);
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,
int bs0,
const DofMap& dofmap1,
int bs1,
75 const std::vector<bool>& bc0,
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 array2d<T>& coeffs,
const std::vector<int>& offsets,
79 const std::vector<T>& constants,
80 const std::vector<std::uint32_t>& cell_info,
81 const std::vector<std::uint8_t>& perms);
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();
93 const int tdim = mesh->topology().dim();
94 const std::int32_t num_cells
95 = mesh->topology().connectivity(tdim, 0)->num_nodes();
98 std::shared_ptr<const fem::DofMap> dofmap0
100 std::shared_ptr<const fem::DofMap> dofmap1
105 const int bs0 = dofmap0->bs();
107 const int bs1 = dofmap1->bs();
116 if (needs_permutation_data)
117 mesh->topology_mutable().create_entity_permutations();
118 const std::vector<std::uint32_t>& cell_info
119 = needs_permutation_data ? mesh->topology().get_cell_permutation_info()
120 : std::vector<std::uint32_t>(num_cells);
124 const auto& fn = a.
kernel(IntegralType::cell, i);
125 const std::vector<std::int32_t>& active_cells
126 = a.
domains(IntegralType::cell, i);
127 impl::assemble_cells<T>(mat_set_values, mesh->geometry(), active_cells,
128 dofs0, bs0, dofs1, bs1, bc0, bc1, fn, coeffs,
129 constants, cell_info);
135 mesh->topology_mutable().create_connectivity(tdim - 1, tdim);
136 mesh->topology_mutable().create_entity_permutations();
138 const std::vector<std::uint8_t>& perms
139 = mesh->topology().get_facet_permutations();
141 for (
int i : a.
integral_ids(IntegralType::exterior_facet))
143 const auto& fn = a.
kernel(IntegralType::exterior_facet, i);
144 const std::vector<std::int32_t>& active_facets
145 = a.
domains(IntegralType::exterior_facet, i);
146 impl::assemble_exterior_facets<T>(mat_set_values, *mesh, active_facets,
147 dofs0, bs0, dofs1, bs1, bc0, bc1, fn,
148 coeffs, constants, cell_info, perms);
152 for (
int i : a.
integral_ids(IntegralType::interior_facet))
154 const auto& fn = a.
kernel(IntegralType::interior_facet, i);
155 const std::vector<std::int32_t>& active_facets
156 = a.
domains(IntegralType::interior_facet, i);
157 impl::assemble_interior_facets<T>(
158 mat_set_values, *mesh, active_facets, *dofmap0, bs0, *dofmap1, bs1,
159 bc0, bc1, fn, coeffs, c_offsets, constants, cell_info, perms);
164 template <
typename T>
166 const std::function<
int(std::int32_t,
const std::int32_t*, std::int32_t,
167 const std::int32_t*,
const T*)>& mat_set,
169 const std::vector<std::int32_t>& active_cells,
172 const std::vector<bool>& bc0,
const std::vector<bool>& bc1,
173 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
174 const std::uint8_t*,
const std::uint32_t)>& kernel,
175 const array2d<T>& coeffs,
const std::vector<T>& constants,
176 const std::vector<std::uint32_t>& cell_info)
178 const int gdim = geometry.
dim();
184 const int num_dofs_g = x_dofmap.
num_links(0);
188 const int num_dofs0 = dofmap0.
links(0).size();
189 const int num_dofs1 = dofmap1.
links(0).size();
190 const int ndim0 = bs0 * num_dofs0;
191 const int ndim1 = bs1 * num_dofs1;
192 std::vector<T> Ae(ndim0 * ndim1);
193 std::vector<double> coordinate_dofs(num_dofs_g * gdim);
194 for (std::int32_t c : active_cells)
197 auto x_dofs = x_dofmap.
links(c);
198 for (std::size_t i = 0; i < x_dofs.size(); ++i)
200 std::copy_n(x_g.
row(x_dofs[i]).begin(), gdim,
201 std::next(coordinate_dofs.begin(), i * gdim));
205 std::fill(Ae.begin(), Ae.end(), 0);
206 kernel(Ae.data(), coeffs.
row(c).data(), constants.data(),
207 coordinate_dofs.data(),
nullptr,
nullptr, cell_info[c]);
210 auto dofs0 = dofmap0.
links(c);
211 auto dofs1 = dofmap1.
links(c);
214 for (
int i = 0; i < num_dofs0; ++i)
216 for (
int k = 0; k < bs0; ++k)
218 if (bc0[bs0 * dofs0[i] + k])
221 const int row = bs0 * i + k;
222 std::fill_n(std::next(Ae.begin(), ndim1 * row), ndim1, 0.0);
230 for (
int j = 0; j < num_dofs1; ++j)
232 for (
int k = 0; k < bs1; ++k)
234 if (bc1[bs1 * dofs1[j] + k])
237 const int col = bs1 * j + k;
238 for (
int row = 0; row < ndim0; ++row)
239 Ae[row * ndim1 + col] = 0.0;
245 mat_set(dofs0.size(), dofs0.data(), dofs1.size(), dofs1.data(), Ae.data());
249 template <
typename T>
250 void assemble_exterior_facets(
251 const std::function<
int(std::int32_t,
const std::int32_t*, std::int32_t,
252 const std::int32_t*,
const T*)>& mat_set_values,
253 const mesh::Mesh& mesh,
const std::vector<std::int32_t>& active_facets,
256 const std::vector<bool>& bc0,
const std::vector<bool>& bc1,
257 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
258 const std::uint8_t*,
const std::uint32_t)>& kernel,
259 const array2d<T>& coeffs,
const std::vector<T>& constants,
260 const std::vector<std::uint32_t>& cell_info,
261 const std::vector<std::uint8_t>& perms)
270 const int num_dofs_g = x_dofmap.
num_links(0);
274 std::vector<double> coordinate_dofs(num_dofs_g * gdim);
275 const int num_dofs0 = dofmap0.
links(0).size();
276 const int num_dofs1 = dofmap1.
links(0).size();
277 const int ndim0 = bs0 * num_dofs0;
278 const int ndim1 = bs1 * num_dofs1;
279 std::vector<T> Ae(ndim0 * ndim1);
286 for (std::int32_t f : active_facets)
288 auto cells = f_to_c->links(f);
289 assert(cells.size() == 1);
292 auto facets = c_to_f->links(cells[0]);
293 auto it = std::find(facets.begin(), facets.end(), f);
294 assert(it != facets.end());
295 const int local_facet = std::distance(facets.begin(), it);
298 auto x_dofs = x_dofmap.
links(cells[0]);
299 for (std::size_t i = 0; i < x_dofs.size(); ++i)
301 std::copy_n(x_g.
row(x_dofs[i]).data(), gdim,
302 std::next(coordinate_dofs.begin(), i * gdim));
306 std::fill(Ae.begin(), Ae.end(), 0);
307 kernel(Ae.data(), coeffs.
row(cells[0]).data(), constants.data(),
308 coordinate_dofs.data(), &local_facet,
309 &perms[cells[0] * facets.size() + local_facet], cell_info[cells[0]]);
312 auto dofs0 = dofmap0.
links(cells[0]);
313 auto dofs1 = dofmap1.
links(cells[0]);
316 for (
int i = 0; i < num_dofs0; ++i)
318 for (
int k = 0; k < bs0; ++k)
320 if (bc0[bs0 * dofs0[i] + k])
323 const int row = bs0 * i + k;
324 std::fill_n(std::next(Ae.begin(), ndim1 * row), ndim1, 0.0);
331 for (
int j = 0; j < num_dofs1; ++j)
333 for (
int k = 0; k < bs1; ++k)
335 if (bc1[bs1 * dofs1[j] + k])
338 const int col = bs1 * j + k;
339 for (
int row = 0; row < ndim0; ++row)
340 Ae[row * ndim1 + col] = 0.0;
346 mat_set_values(dofs0.size(), dofs0.data(), dofs1.size(), dofs1.data(),
351 template <
typename T>
352 void assemble_interior_facets(
353 const std::function<
int(std::int32_t,
const std::int32_t*, std::int32_t,
354 const std::int32_t*,
const T*)>& mat_set_values,
355 const mesh::Mesh& mesh,
const std::vector<std::int32_t>& active_facets,
356 const DofMap& dofmap0,
int bs0,
const DofMap& dofmap1,
int bs1,
357 const std::vector<bool>& bc0,
const std::vector<bool>& bc1,
358 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
359 const std::uint8_t*,
const std::uint32_t)>& fn,
360 const array2d<T>& coeffs,
const std::vector<int>& offsets,
361 const std::vector<T>& constants,
362 const std::vector<std::uint32_t>& cell_info,
363 const std::vector<std::uint8_t>& perms)
372 const int num_dofs_g = x_dofmap.
num_links(0);
376 std::vector<double> coordinate_dofs(2 * num_dofs_g * gdim);
377 std::vector<T> Ae, be;
378 std::vector<T> coeff_array(2 * offsets.back());
379 assert(offsets.back() == coeffs.
shape[1]);
382 std::vector<std::int32_t> dmapjoint0, dmapjoint1;
389 const int offset_g = gdim * num_dofs_g;
390 for (std::int32_t facet_index : active_facets)
393 auto cells = c->links(facet_index);
394 assert(cells.size() == 2);
397 auto facets0 = c_to_f->links(cells[0]);
398 const auto* it0 = std::find(facets0.begin(), facets0.end(), facet_index);
399 assert(it0 != facets0.end());
400 const int local_facet0 = std::distance(facets0.begin(), it0);
401 auto facets1 = c_to_f->links(cells[1]);
402 const auto* it1 = std::find(facets1.begin(), facets1.end(), facet_index);
403 assert(it1 != facets1.end());
404 const int local_facet1 = std::distance(facets1.begin(), it1);
406 const std::array local_facet{local_facet0, local_facet1};
409 auto x_dofs0 = x_dofmap.
links(cells[0]);
410 auto x_dofs1 = x_dofmap.
links(cells[1]);
411 for (
int i = 0; i < num_dofs_g; ++i)
413 for (
int j = 0; j < gdim; ++j)
415 coordinate_dofs[i * gdim + j] = x_g(x_dofs0[i], j);
416 coordinate_dofs[offset_g + i * gdim + j] = x_g(x_dofs1[i], j);
421 tcb::span<const std::int32_t> dmap0_cell0 = dofmap0.cell_dofs(cells[0]);
422 tcb::span<const std::int32_t> dmap0_cell1 = dofmap0.cell_dofs(cells[1]);
423 dmapjoint0.resize(dmap0_cell0.size() + dmap0_cell1.size());
424 std::copy(dmap0_cell0.begin(), dmap0_cell0.end(), dmapjoint0.begin());
425 std::copy(dmap0_cell1.begin(), dmap0_cell1.end(),
426 std::next(dmapjoint0.begin(), dmap0_cell0.size()));
428 tcb::span<const std::int32_t> dmap1_cell0 = dofmap1.cell_dofs(cells[0]);
429 tcb::span<const std::int32_t> dmap1_cell1 = dofmap1.cell_dofs(cells[1]);
430 dmapjoint1.resize(dmap1_cell0.size() + dmap1_cell1.size());
431 std::copy(dmap1_cell0.begin(), dmap1_cell0.end(), dmapjoint1.begin());
432 std::copy(dmap1_cell1.begin(), dmap1_cell1.end(),
433 std::next(dmapjoint1.begin(), dmap1_cell0.size()));
437 auto coeff_cell0 = coeffs.
row(cells[0]);
438 auto coeff_cell1 = coeffs.
row(cells[1]);
441 for (std::size_t i = 0; i < offsets.size() - 1; ++i)
444 const int num_entries = offsets[i + 1] - offsets[i];
445 std::copy_n(coeff_cell0.data() + offsets[i], num_entries,
446 std::next(coeff_array.begin(), 2 * offsets[i]));
447 std::copy_n(coeff_cell1.data() + offsets[i], num_entries,
448 std::next(coeff_array.begin(), offsets[i + 1] + offsets[i]));
451 const int num_rows = bs0 * dmapjoint0.size();
452 const int num_cols = bs1 * dmapjoint1.size();
455 Ae.resize(num_rows * num_cols);
456 std::fill(Ae.begin(), Ae.end(), 0);
458 const int facets_per_cell = facets0.size();
459 const std::array perm{perms[cells[0] * facets_per_cell + local_facet[0]],
460 perms[cells[1] * facets_per_cell + local_facet[1]]};
461 fn(Ae.data(), coeff_array.data(), constants.data(), coordinate_dofs.data(),
462 local_facet.data(), perm.data(), cell_info[cells[0]]);
467 for (std::size_t i = 0; i < dmapjoint0.size(); ++i)
469 for (
int k = 0; k < bs0; ++k)
471 if (bc0[bs0 * dmapjoint0[i] + k])
474 std::fill_n(std::next(Ae.begin(), num_cols * (bs0 * i + k)),
482 for (std::size_t j = 0; j < dmapjoint1.size(); ++j)
484 for (
int k = 0; k < bs1; ++k)
486 if (bc1[bs1 * dmapjoint1[j] + k])
489 for (
int m = 0; m < num_rows; ++m)
490 Ae[m * num_cols + bs1 * j + k] = 0.0;
496 mat_set_values(dmapjoint0.size(), dmapjoint0.data(), dmapjoint1.size(),
497 dmapjoint1.data(), Ae.data());
This class provides a dynamic 2-dimensional row-wise array data structure.
Definition: array2d.h:21
std::array< size_type, 2 > shape
The shape of the array.
Definition: array2d.h:157
constexpr tcb::span< value_type > row(size_type i)
Access a row in the array.
Definition: array2d.h:116
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:68
tcb::span< T > links(int node)
Get the links (edges) for given node.
Definition: AdjacencyList.h:151
int num_links(int node) const
Number of connections for given node.
Definition: AdjacencyList.h:141
Geometry stores the geometry imposed on a mesh.
Definition: Geometry.h:35
const graph::AdjacencyList< std::int32_t > & dofmap() const
DOF map.
Definition: Geometry.cpp:21
array2d< double > & x()
Geometry degrees-of-freedom.
Definition: Geometry.cpp:31
int dim() const
Return Euclidean dimension of coordinate system.
Definition: Geometry.cpp:19
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:122
Topology & topology()
Get mesh topology.
Definition: Mesh.cpp:116
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:263
int dim() const
Return topological dimension.
Definition: Topology.cpp:162
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:453
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
array2d< typename U::scalar_type > pack_coefficients(const U &u)
Pack coefficients of u of generic type U ready for assembly.
Definition: utils.h:398