DOLFINx
DOLFINx C++ interface
assemble_matrix_impl.h
1// Copyright (C) 2018-2019 Garth N. Wells
2//
3// This file is part of DOLFINx (https://www.fenicsproject.org)
4//
5// SPDX-License-Identifier: LGPL-3.0-or-later
6
7#pragma once
8
9#include "DofMap.h"
10#include "Form.h"
11#include "FunctionSpace.h"
12#include "utils.h"
13#include <dolfinx/common/utils.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>
19#include <functional>
20#include <iterator>
21#include <span>
22#include <vector>
23
24namespace dolfinx::fem::impl
25{
26
32template <typename T, typename U>
34 U mat_set_values, const Form<T>& a, const std::span<const T>& constants,
35 const std::map<std::pair<IntegralType, int>,
36 std::pair<std::span<const T>, int>>& coefficients,
37 const std::span<const std::int8_t>& bc0,
38 const std::span<const std::int8_t>& bc1);
39
41template <typename T, typename U>
42void assemble_cells(
43 U mat_set, const mesh::Geometry& geometry,
44 const std::span<const std::int32_t>& cells,
45 const std::function<void(const std::span<T>&,
46 const std::span<const std::uint32_t>&,
47 std::int32_t, int)>& dof_transform,
48 const graph::AdjacencyList<std::int32_t>& dofmap0, int bs0,
49 const std::function<void(const std::span<T>&,
50 const std::span<const std::uint32_t>&,
51 std::int32_t, int)>& dof_transform_to_transpose,
52 const graph::AdjacencyList<std::int32_t>& dofmap1, int bs1,
53 const std::span<const std::int8_t>& bc0,
54 const std::span<const std::int8_t>& bc1,
55 const std::function<void(T*, const T*, const T*,
56 const scalar_value_type_t<T>*, const int*,
57 const std::uint8_t*)>& kernel,
58 const std::span<const T>& coeffs, int cstride,
59 const std::span<const T>& constants,
60 const std::span<const std::uint32_t>& cell_info)
61{
62 if (cells.empty())
63 return;
64
65 // Prepare cell geometry
66 const graph::AdjacencyList<std::int32_t>& x_dofmap = geometry.dofmap();
67 const std::size_t num_dofs_g = geometry.cmap().dim();
68 std::span<const double> x_g = geometry.x();
69
70 // Iterate over active cells
71 const int num_dofs0 = dofmap0.links(0).size();
72 const int num_dofs1 = dofmap1.links(0).size();
73 const int ndim0 = bs0 * num_dofs0;
74 const int ndim1 = bs1 * num_dofs1;
75 std::vector<T> Ae(ndim0 * ndim1);
76 const std::span<T> _Ae(Ae);
77 std::vector<scalar_value_type_t<T>> coordinate_dofs(3 * num_dofs_g);
78
79 // Iterate over active cells
80 for (std::size_t index = 0; index < cells.size(); ++index)
81 {
82 std::int32_t c = cells[index];
83
84 // Get cell coordinates/geometry
85 auto x_dofs = x_dofmap.links(c);
86 for (std::size_t i = 0; i < x_dofs.size(); ++i)
87 {
88 common::impl::copy_N<3>(std::next(x_g.begin(), 3 * x_dofs[i]),
89 std::next(coordinate_dofs.begin(), 3 * i));
90 }
91
92 // Tabulate tensor
93 std::fill(Ae.begin(), Ae.end(), 0);
94 kernel(Ae.data(), coeffs.data() + index * cstride, constants.data(),
95 coordinate_dofs.data(), nullptr, nullptr);
96
97 dof_transform(_Ae, cell_info, c, ndim1);
98 dof_transform_to_transpose(_Ae, cell_info, c, ndim0);
99
100 // Zero rows/columns for essential bcs
101 auto dofs0 = dofmap0.links(c);
102 auto dofs1 = dofmap1.links(c);
103 if (!bc0.empty())
104 {
105 for (int i = 0; i < num_dofs0; ++i)
106 {
107 for (int k = 0; k < bs0; ++k)
108 {
109 if (bc0[bs0 * dofs0[i] + k])
110 {
111 // Zero row bs0 * i + k
112 const int row = bs0 * i + k;
113 std::fill_n(std::next(Ae.begin(), ndim1 * row), ndim1, 0.0);
114 }
115 }
116 }
117 }
118
119 if (!bc1.empty())
120 {
121 for (int j = 0; j < num_dofs1; ++j)
122 {
123 for (int k = 0; k < bs1; ++k)
124 {
125 if (bc1[bs1 * dofs1[j] + k])
126 {
127 // Zero column bs1 * j + k
128 const int col = bs1 * j + k;
129 for (int row = 0; row < ndim0; ++row)
130 Ae[row * ndim1 + col] = 0.0;
131 }
132 }
133 }
134 }
135
136 mat_set(dofs0, dofs1, Ae);
137 }
138}
139
141template <typename T, typename U>
142void assemble_exterior_facets(
143 U mat_set, const mesh::Mesh& mesh,
144 const std::span<const std::int32_t>& facets,
145 const std::function<void(const std::span<T>&,
146 const std::span<const std::uint32_t>&,
147 std::int32_t, int)>& dof_transform,
148 const graph::AdjacencyList<std::int32_t>& dofmap0, int bs0,
149 const std::function<void(const std::span<T>&,
150 const std::span<const std::uint32_t>&,
151 std::int32_t, int)>& dof_transform_to_transpose,
152 const graph::AdjacencyList<std::int32_t>& dofmap1, int bs1,
153 const std::span<const std::int8_t>& bc0,
154 const std::span<const std::int8_t>& bc1,
155 const std::function<void(T*, const T*, const T*,
156 const scalar_value_type_t<T>*, const int*,
157 const std::uint8_t*)>& kernel,
158 const std::span<const T>& coeffs, int cstride,
159 const std::span<const T>& constants,
160 const std::span<const std::uint32_t>& cell_info)
161{
162 if (facets.empty())
163 return;
164
165 // Prepare cell geometry
166 const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
167 const std::size_t num_dofs_g = mesh.geometry().cmap().dim();
168 std::span<const double> x_g = mesh.geometry().x();
169
170 // Data structures used in assembly
171 std::vector<scalar_value_type_t<T>> coordinate_dofs(3 * num_dofs_g);
172 const int num_dofs0 = dofmap0.links(0).size();
173 const int num_dofs1 = dofmap1.links(0).size();
174 const int ndim0 = bs0 * num_dofs0;
175 const int ndim1 = bs1 * num_dofs1;
176 std::vector<T> Ae(ndim0 * ndim1);
177 const std::span<T> _Ae(Ae);
178 assert(facets.size() % 2 == 0);
179 for (std::size_t index = 0; index < facets.size(); index += 2)
180 {
181 std::int32_t cell = facets[index];
182 std::int32_t local_facet = facets[index + 1];
183
184 // Get cell coordinates/geometry
185 auto x_dofs = x_dofmap.links(cell);
186 for (std::size_t i = 0; i < x_dofs.size(); ++i)
187 {
188 common::impl::copy_N<3>(std::next(x_g.begin(), 3 * x_dofs[i]),
189 std::next(coordinate_dofs.begin(), 3 * i));
190 }
191
192 // Tabulate tensor
193 std::fill(Ae.begin(), Ae.end(), 0);
194 kernel(Ae.data(), coeffs.data() + index / 2 * cstride, constants.data(),
195 coordinate_dofs.data(), &local_facet, nullptr);
196
197 dof_transform(_Ae, cell_info, cell, ndim1);
198 dof_transform_to_transpose(_Ae, cell_info, cell, ndim0);
199
200 // Zero rows/columns for essential bcs
201 auto dofs0 = dofmap0.links(cell);
202 auto dofs1 = dofmap1.links(cell);
203 if (!bc0.empty())
204 {
205 for (int i = 0; i < num_dofs0; ++i)
206 {
207 for (int k = 0; k < bs0; ++k)
208 {
209 if (bc0[bs0 * dofs0[i] + k])
210 {
211 // Zero row bs0 * i + k
212 const int row = bs0 * i + k;
213 std::fill_n(std::next(Ae.begin(), ndim1 * row), ndim1, 0.0);
214 }
215 }
216 }
217 }
218 if (!bc1.empty())
219 {
220 for (int j = 0; j < num_dofs1; ++j)
221 {
222 for (int k = 0; k < bs1; ++k)
223 {
224 if (bc1[bs1 * dofs1[j] + k])
225 {
226 // Zero column bs1 * j + k
227 const int col = bs1 * j + k;
228 for (int row = 0; row < ndim0; ++row)
229 Ae[row * ndim1 + col] = 0.0;
230 }
231 }
232 }
233 }
234
235 mat_set(dofs0, dofs1, Ae);
236 }
237}
238
240template <typename T, typename U>
241void assemble_interior_facets(
242 U mat_set, const mesh::Mesh& mesh,
243 const std::span<const std::int32_t>& facets,
244 const std::function<void(const std::span<T>&,
245 const std::span<const std::uint32_t>&,
246 std::int32_t, int)>& dof_transform,
247 const DofMap& dofmap0, int bs0,
248 const std::function<void(const std::span<T>&,
249 const std::span<const std::uint32_t>&,
250 std::int32_t, int)>& dof_transform_to_transpose,
251 const DofMap& dofmap1, int bs1, const std::span<const std::int8_t>& bc0,
252 const std::span<const std::int8_t>& bc1,
253 const std::function<void(T*, const T*, const T*,
254 const scalar_value_type_t<T>*, const int*,
255 const std::uint8_t*)>& kernel,
256 const std::span<const T>& coeffs, int cstride,
257 const std::span<const int>& offsets, const std::span<const T>& constants,
258 const std::span<const std::uint32_t>& cell_info,
259 const std::function<std::uint8_t(std::size_t)>& get_perm)
260{
261 if (facets.empty())
262 return;
263
264 const int tdim = mesh.topology().dim();
265
266 // Prepare cell geometry
267 const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
268 const std::size_t num_dofs_g = mesh.geometry().cmap().dim();
269 std::span<const double> x_g = mesh.geometry().x();
270
271 // Data structures used in assembly
272 using X = scalar_value_type_t<T>;
273 std::vector<X> coordinate_dofs(2 * num_dofs_g * 3);
274 std::span<X> cdofs0(coordinate_dofs.data(), num_dofs_g * 3);
275 std::span<X> cdofs1(coordinate_dofs.data() + num_dofs_g * 3, num_dofs_g * 3);
276
277 std::vector<T> Ae, be;
278 std::vector<T> coeff_array(2 * offsets.back());
279 assert(offsets.back() == cstride);
280
281 const int num_cell_facets
282 = mesh::cell_num_entities(mesh.topology().cell_type(), tdim - 1);
283
284 // Temporaries for joint dofmaps
285 std::vector<std::int32_t> dmapjoint0, dmapjoint1;
286 assert(facets.size() % 4 == 0);
287 for (std::size_t index = 0; index < facets.size(); index += 4)
288 {
289 std::array<std::int32_t, 2> cells = {facets[index], facets[index + 2]};
290 std::array<std::int32_t, 2> local_facet
291 = {facets[index + 1], facets[index + 3]};
292
293 // Get cell geometry
294 auto x_dofs0 = x_dofmap.links(cells[0]);
295 for (std::size_t i = 0; i < x_dofs0.size(); ++i)
296 {
297 common::impl::copy_N<3>(std::next(x_g.begin(), 3 * x_dofs0[i]),
298 std::next(cdofs0.begin(), 3 * i));
299 }
300 auto x_dofs1 = x_dofmap.links(cells[1]);
301 for (std::size_t i = 0; i < x_dofs1.size(); ++i)
302 {
303 common::impl::copy_N<3>(std::next(x_g.begin(), 3 * x_dofs1[i]),
304 std::next(cdofs1.begin(), 3 * i));
305 }
306
307 // Get dof maps for cells and pack
308 std::span<const std::int32_t> dmap0_cell0 = dofmap0.cell_dofs(cells[0]);
309 std::span<const std::int32_t> dmap0_cell1 = dofmap0.cell_dofs(cells[1]);
310 dmapjoint0.resize(dmap0_cell0.size() + dmap0_cell1.size());
311 std::copy(dmap0_cell0.begin(), dmap0_cell0.end(), dmapjoint0.begin());
312 std::copy(dmap0_cell1.begin(), dmap0_cell1.end(),
313 std::next(dmapjoint0.begin(), dmap0_cell0.size()));
314
315 std::span<const std::int32_t> dmap1_cell0 = dofmap1.cell_dofs(cells[0]);
316 std::span<const std::int32_t> dmap1_cell1 = dofmap1.cell_dofs(cells[1]);
317 dmapjoint1.resize(dmap1_cell0.size() + dmap1_cell1.size());
318 std::copy(dmap1_cell0.begin(), dmap1_cell0.end(), dmapjoint1.begin());
319 std::copy(dmap1_cell1.begin(), dmap1_cell1.end(),
320 std::next(dmapjoint1.begin(), dmap1_cell0.size()));
321
322 const int num_rows = bs0 * dmapjoint0.size();
323 const int num_cols = bs1 * dmapjoint1.size();
324
325 // Tabulate tensor
326 Ae.resize(num_rows * num_cols);
327 std::fill(Ae.begin(), Ae.end(), 0);
328
329 const std::array perm{
330 get_perm(cells[0] * num_cell_facets + local_facet[0]),
331 get_perm(cells[1] * num_cell_facets + local_facet[1])};
332 kernel(Ae.data(), coeffs.data() + index / 2 * cstride, constants.data(),
333 coordinate_dofs.data(), local_facet.data(), perm.data());
334
335 const std::span<T> _Ae(Ae);
336
337 const std::span<T> sub_Ae0
338 = _Ae.subspan(bs0 * dmap0_cell0.size() * num_cols,
339 bs0 * dmap0_cell1.size() * num_cols);
340 const std::span<T> sub_Ae1
341 = _Ae.subspan(bs1 * dmap1_cell0.size(),
342 num_rows * num_cols - bs1 * dmap1_cell0.size());
343
344 // Need to apply DOF transformations for parts of the matrix due to cell 0
345 // and cell 1. For example, if the space has 3 DOFs, then Ae will be 6 by 6
346 // (3 rows/columns for each cell). Subspans are used to offset to the right
347 // blocks of the matrix
348
349 dof_transform(_Ae, cell_info, cells[0], num_cols);
350 dof_transform(sub_Ae0, cell_info, cells[1], num_cols);
351 dof_transform_to_transpose(_Ae, cell_info, cells[0], num_rows);
352 dof_transform_to_transpose(sub_Ae1, cell_info, cells[1], num_rows);
353
354 // Zero rows/columns for essential bcs
355 if (!bc0.empty())
356 {
357 for (std::size_t i = 0; i < dmapjoint0.size(); ++i)
358 {
359 for (int k = 0; k < bs0; ++k)
360 {
361 if (bc0[bs0 * dmapjoint0[i] + k])
362 {
363 // Zero row bs0 * i + k
364 std::fill_n(std::next(Ae.begin(), num_cols * (bs0 * i + k)),
365 num_cols, 0.0);
366 }
367 }
368 }
369 }
370 if (!bc1.empty())
371 {
372 for (std::size_t j = 0; j < dmapjoint1.size(); ++j)
373 {
374 for (int k = 0; k < bs1; ++k)
375 {
376 if (bc1[bs1 * dmapjoint1[j] + k])
377 {
378 // Zero column bs1 * j + k
379 for (int m = 0; m < num_rows; ++m)
380 Ae[m * num_cols + bs1 * j + k] = 0.0;
381 }
382 }
383 }
384 }
385
386 mat_set(dmapjoint0, dmapjoint1, Ae);
387 }
388}
389
390template <typename T, typename U>
391void assemble_matrix(
392 U mat_set, const Form<T>& a, const std::span<const T>& constants,
393 const std::map<std::pair<IntegralType, int>,
394 std::pair<std::span<const T>, int>>& coefficients,
395 const std::span<const std::int8_t>& bc0,
396 const std::span<const std::int8_t>& bc1)
397{
398 std::shared_ptr<const mesh::Mesh> mesh = a.mesh();
399 assert(mesh);
400
401 // Get dofmap data
402 std::shared_ptr<const fem::DofMap> dofmap0
403 = a.function_spaces().at(0)->dofmap();
404 std::shared_ptr<const fem::DofMap> dofmap1
405 = a.function_spaces().at(1)->dofmap();
406 assert(dofmap0);
407 assert(dofmap1);
408 const graph::AdjacencyList<std::int32_t>& dofs0 = dofmap0->list();
409 const int bs0 = dofmap0->bs();
410 const graph::AdjacencyList<std::int32_t>& dofs1 = dofmap1->list();
411 const int bs1 = dofmap1->bs();
412
413 std::shared_ptr<const fem::FiniteElement> element0
414 = a.function_spaces().at(0)->element();
415 std::shared_ptr<const fem::FiniteElement> element1
416 = a.function_spaces().at(1)->element();
417 const std::function<void(const std::span<T>&,
418 const std::span<const std::uint32_t>&, std::int32_t,
419 int)>& dof_transform
420 = element0->get_dof_transformation_function<T>();
421 const std::function<void(const std::span<T>&,
422 const std::span<const std::uint32_t>&, std::int32_t,
423 int)>& dof_transform_to_transpose
424 = element1->get_dof_transformation_to_transpose_function<T>();
425
426 const bool needs_transformation_data
427 = element0->needs_dof_transformations()
428 or element1->needs_dof_transformations()
430 std::span<const std::uint32_t> cell_info;
431 if (needs_transformation_data)
432 {
433 mesh->topology_mutable().create_entity_permutations();
434 cell_info = std::span(mesh->topology().get_cell_permutation_info());
435 }
436
437 for (int i : a.integral_ids(IntegralType::cell))
438 {
439 const auto& fn = a.kernel(IntegralType::cell, i);
440 const auto& [coeffs, cstride] = coefficients.at({IntegralType::cell, i});
441 const std::vector<std::int32_t>& cells = a.cell_domains(i);
442 impl::assemble_cells(mat_set, mesh->geometry(), cells, dof_transform, dofs0,
443 bs0, dof_transform_to_transpose, dofs1, bs1, bc0, bc1,
444 fn, coeffs, cstride, constants, cell_info);
445 }
446
448 {
449 const auto& fn = a.kernel(IntegralType::exterior_facet, i);
450 const auto& [coeffs, cstride]
451 = coefficients.at({IntegralType::exterior_facet, i});
452 const std::vector<std::int32_t>& facets = a.exterior_facet_domains(i);
453 impl::assemble_exterior_facets(mat_set, *mesh, facets, dof_transform, dofs0,
454 bs0, dof_transform_to_transpose, dofs1, bs1,
455 bc0, bc1, fn, coeffs, cstride, constants,
456 cell_info);
457 }
458
460 {
461 std::function<std::uint8_t(std::size_t)> get_perm;
463 {
464 mesh->topology_mutable().create_entity_permutations();
465 const std::vector<std::uint8_t>& perms
466 = mesh->topology().get_facet_permutations();
467 get_perm = [&perms](std::size_t i) { return perms[i]; };
468 }
469 else
470 get_perm = [](std::size_t) { return 0; };
471
472 const std::vector<int> c_offsets = a.coefficient_offsets();
474 {
475 const auto& fn = a.kernel(IntegralType::interior_facet, i);
476 const auto& [coeffs, cstride]
477 = coefficients.at({IntegralType::interior_facet, i});
478 const std::vector<std::int32_t>& facets = a.interior_facet_domains(i);
479 impl::assemble_interior_facets(
480 mat_set, *mesh, facets, dof_transform, *dofmap0, bs0,
481 dof_transform_to_transpose, *dofmap1, bs1, bc0, bc1, fn, coeffs,
482 cstride, c_offsets, constants, cell_info, get_perm);
483 }
484 }
485}
486
487} // namespace dolfinx::fem::impl
Degree-of-freedeom map representations ans tools.
int dim() const
The dimension of the geometry element space.
Definition: CoordinateElement.cpp:183
Degree-of-freedom map.
Definition: DofMap.h:71
A representation of finite element variational forms.
Definition: Form.h:63
bool needs_facet_permutations() const
Get bool indicating whether permutation data needs to be passed into these integrals.
Definition: Form.h:325
const std::vector< std::int32_t > & cell_domains(int i) const
Get the list of cell indices for the ith integral (kernel) for the cell domain type.
Definition: Form.h:280
const std::vector< std::int32_t > & exterior_facet_domains(int i) const
Get the list of (cell_index, local_facet_index) pairs for the ith integral (kernel) for the exterior ...
Definition: Form.h:293
std::vector< int > coefficient_offsets() const
Offset for each coefficient expansion array on a cell. Used to pack data for multiple coefficients in...
Definition: Form.h:330
const std::function< void(T *, const T *, const T *, const scalar_value_type_t *, const int *, const std::uint8_t *)> & kernel(IntegralType type, int i) const
Get the function for 'kernel' for integral i of given type.
Definition: Form.h:193
const std::vector< std::shared_ptr< const FunctionSpace > > & function_spaces() const
Return function spaces for all arguments.
Definition: Form.h:181
std::vector< int > integral_ids(IntegralType type) const
Get the IDs for integrals (kernels) for given integral type. The IDs correspond to the domain IDs whi...
Definition: Form.h:248
std::shared_ptr< const mesh::Mesh > mesh() const
Extract common mesh for the form.
Definition: Form.h:176
const std::vector< std::int32_t > & interior_facet_domains(int i) const
Get the list of (cell_index_0, local_facet_index_0, cell_index_1, local_facet_index_1) quadruplets fo...
Definition: Form.h:308
int num_integrals(IntegralType type) const
Number of integrals of given type.
Definition: Form.h:227
This class provides a static adjacency list data structure. It is commonly used to store directed gra...
Definition: AdjacencyList.h:26
std::span< T > links(int node)
Get the links (edges) for given node.
Definition: AdjacencyList.h:111
Geometry stores the geometry imposed on a mesh.
Definition: Geometry.h:28
const graph::AdjacencyList< std::int32_t > & dofmap() const
DOF map.
Definition: Geometry.cpp:21
std::span< const double > x() const
Access geometry degrees-of-freedom data (const version).
Definition: Geometry.cpp:33
const fem::CoordinateElement & cmap() const
The element that describes the geometry map.
Definition: Geometry.cpp:35
A Mesh consists of a set of connected and numbered mesh topological entities, and geometry data.
Definition: Mesh.h:33
Geometry & geometry()
Get mesh geometry.
Definition: Mesh.cpp:450
Topology & topology()
Get mesh topology.
Definition: Mesh.cpp:444
int dim() const noexcept
Return the topological dimension of the mesh.
Definition: Topology.cpp:776
CellType cell_type() const noexcept
Cell type.
Definition: Topology.cpp:910
void assemble_matrix(U mat_add, const Form< T > &a, const std::span< const T > &constants, const std::map< std::pair< IntegralType, int >, std::pair< std::span< const T >, int > > &coefficients, const std::vector< std::shared_ptr< const DirichletBC< T > > > &bcs)
Assemble bilinear form into a matrix.
Definition: assembler.h:200
@ interior_facet
Interior facet.
@ exterior_facet
Exterior facet.
int cell_num_entities(CellType type, int dim)
Number of entities of dimension dim.
Definition: cell_types.cpp:182