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