DOLFINx
DOLFINx C++ interface
assemble_vector_impl.h
1// Copyright (C) 2018-2021 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 "Constant.h"
10#include "DirichletBC.h"
11#include "DofMap.h"
12#include "Form.h"
13#include "FunctionSpace.h"
14#include "utils.h"
15#include <dolfinx/common/IndexMap.h>
16#include <dolfinx/common/utils.h>
17#include <dolfinx/graph/AdjacencyList.h>
18#include <dolfinx/mesh/Geometry.h>
19#include <dolfinx/mesh/Mesh.h>
20#include <dolfinx/mesh/Topology.h>
21#include <functional>
22#include <memory>
23#include <span>
24#include <vector>
25
26namespace dolfinx::fem::impl
27{
28
30
38template <typename T, int _bs0 = -1, int _bs1 = -1>
39void _lift_bc_cells(
40 std::span<T> b, const mesh::Geometry& geometry,
41 const std::function<void(T*, const T*, const T*,
42 const scalar_value_type_t<T>*, const int*,
43 const std::uint8_t*)>& kernel,
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 T>& constants, const std::span<const T>& coeffs,
54 int cstride, const std::span<const std::uint32_t>& cell_info,
55 const std::span<const T>& bc_values1,
56 const std::span<const std::int8_t>& bc_markers1,
57 const std::span<const T>& x0, double scale)
58{
59 assert(_bs0 < 0 or _bs0 == bs0);
60 assert(_bs1 < 0 or _bs1 == bs1);
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 // Data structures used in bc application
71 std::vector<scalar_value_type_t<T>> coordinate_dofs(3 * num_dofs_g);
72 std::vector<T> Ae, be;
73 const scalar_value_type_t<T> _scale
74 = static_cast<scalar_value_type_t<T>>(scale);
75 for (std::size_t index = 0; index < cells.size(); ++index)
76 {
77 std::int32_t c = cells[index];
78
79 // Get dof maps for cell
80 auto dmap1 = dofmap1.links(c);
81
82 // Check if bc is applied to cell
83 bool has_bc = false;
84 for (std::size_t j = 0; j < dmap1.size(); ++j)
85 {
86 if constexpr (_bs1 > 0)
87 {
88 for (int k = 0; k < _bs1; ++k)
89 {
90 assert(_bs1 * dmap1[j] + k < (int)bc_markers1.size());
91 if (bc_markers1[_bs1 * dmap1[j] + k])
92 {
93 has_bc = true;
94 break;
95 }
96 }
97 }
98 else
99 {
100 for (int k = 0; k < bs1; ++k)
101 {
102 assert(bs1 * dmap1[j] + k < (int)bc_markers1.size());
103 if (bc_markers1[bs1 * dmap1[j] + k])
104 {
105 has_bc = true;
106 break;
107 }
108 }
109 }
110 }
111
112 if (!has_bc)
113 continue;
114
115 // Get cell coordinates/geometry
116 auto x_dofs = x_dofmap.links(c);
117 for (std::size_t i = 0; i < x_dofs.size(); ++i)
118 {
119 common::impl::copy_N<3>(std::next(x_g.begin(), 3 * x_dofs[i]),
120 std::next(coordinate_dofs.begin(), 3 * i));
121 }
122
123 // Size data structure for assembly
124 auto dmap0 = dofmap0.links(c);
125
126 const int num_rows = bs0 * dmap0.size();
127 const int num_cols = bs1 * dmap1.size();
128
129 const T* coeff_array = coeffs.data() + index * cstride;
130 Ae.resize(num_rows * num_cols);
131 std::fill(Ae.begin(), Ae.end(), 0);
132 kernel(Ae.data(), coeff_array, constants.data(), coordinate_dofs.data(),
133 nullptr, nullptr);
134 dof_transform(Ae, cell_info, c, num_cols);
135 dof_transform_to_transpose(Ae, cell_info, c, num_rows);
136
137 // Size data structure for assembly
138 be.resize(num_rows);
139 std::fill(be.begin(), be.end(), 0);
140 for (std::size_t j = 0; j < dmap1.size(); ++j)
141 {
142 if constexpr (_bs1 > 0)
143 {
144 for (int k = 0; k < _bs1; ++k)
145 {
146 const std::int32_t jj = _bs1 * dmap1[j] + k;
147 assert(jj < (int)bc_markers1.size());
148 if (bc_markers1[jj])
149 {
150 const T bc = bc_values1[jj];
151 const T _x0 = x0.empty() ? 0.0 : x0[jj];
152 // const T _x0 = 0.0;
153 // be -= Ae.col(bs1 * j + k) * scale * (bc - _x0);
154 for (int m = 0; m < num_rows; ++m)
155 be[m] -= Ae[m * num_cols + _bs1 * j + k] * _scale * (bc - _x0);
156 }
157 }
158 }
159 else
160 {
161 for (int k = 0; k < bs1; ++k)
162 {
163 const std::int32_t jj = bs1 * dmap1[j] + k;
164 assert(jj < (int)bc_markers1.size());
165 if (bc_markers1[jj])
166 {
167 const T bc = bc_values1[jj];
168 const T _x0 = x0.empty() ? 0.0 : x0[jj];
169 // be -= Ae.col(bs1 * j + k) * scale * (bc - _x0);
170 for (int m = 0; m < num_rows; ++m)
171 be[m] -= Ae[m * num_cols + bs1 * j + k] * _scale * (bc - _x0);
172 }
173 }
174 }
175 }
176
177 for (std::size_t i = 0; i < dmap0.size(); ++i)
178 {
179 if constexpr (_bs0 > 0)
180 {
181 for (int k = 0; k < _bs0; ++k)
182 b[_bs0 * dmap0[i] + k] += be[_bs0 * i + k];
183 }
184 else
185 {
186 for (int k = 0; k < bs0; ++k)
187 b[bs0 * dmap0[i] + k] += be[bs0 * i + k];
188 }
189 }
190 }
191}
192
199template <typename T, int _bs = -1>
200void _lift_bc_exterior_facets(
201 std::span<T> b, const mesh::Mesh& mesh,
202 const std::function<void(T*, const T*, const T*,
203 const scalar_value_type_t<T>*, const int*,
204 const std::uint8_t*)>& kernel,
205 const std::span<const std::int32_t>& facets,
206 const std::function<void(const std::span<T>&,
207 const std::span<const std::uint32_t>&,
208 std::int32_t, int)>& dof_transform,
209 const graph::AdjacencyList<std::int32_t>& dofmap0, int bs0,
210 const std::function<void(const std::span<T>&,
211 const std::span<const std::uint32_t>&,
212 std::int32_t, int)>& dof_transform_to_transpose,
213 const graph::AdjacencyList<std::int32_t>& dofmap1, int bs1,
214 const std::span<const T>& constants, const std::span<const T>& coeffs,
215 int cstride, const std::span<const std::uint32_t>& cell_info,
216 const std::span<const T>& bc_values1,
217 const std::span<const std::int8_t>& bc_markers1,
218 const std::span<const T>& x0, double scale)
219{
220 if (facets.empty())
221 return;
222
223 // Prepare cell geometry
224 const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
225 const std::size_t num_dofs_g = mesh.geometry().cmap().dim();
226 std::span<const double> x_g = mesh.geometry().x();
227
228 // Data structures used in bc application
229 std::vector<scalar_value_type_t<T>> coordinate_dofs(3 * num_dofs_g);
230 std::vector<T> Ae, be;
231 assert(facets.size() % 2 == 0);
232 const scalar_value_type_t<T> _scale
233 = static_cast<scalar_value_type_t<T>>(scale);
234 for (std::size_t index = 0; index < facets.size(); index += 2)
235 {
236 std::int32_t cell = facets[index];
237 std::int32_t local_facet = facets[index + 1];
238
239 // Get dof maps for cell
240 auto dmap1 = dofmap1.links(cell);
241
242 // Check if bc is applied to cell
243 bool has_bc = false;
244 for (std::size_t j = 0; j < dmap1.size(); ++j)
245 {
246 for (int k = 0; k < bs1; ++k)
247 {
248 if (bc_markers1[bs1 * dmap1[j] + k])
249 {
250 has_bc = true;
251 break;
252 }
253 }
254 }
255
256 if (!has_bc)
257 continue;
258
259 // Get cell coordinates/geometry
260 auto x_dofs = x_dofmap.links(cell);
261 for (std::size_t i = 0; i < x_dofs.size(); ++i)
262 {
263 common::impl::copy_N<3>(std::next(x_g.begin(), 3 * x_dofs[i]),
264 std::next(coordinate_dofs.begin(), 3 * i));
265 }
266
267 // Size data structure for assembly
268 auto dmap0 = dofmap0.links(cell);
269 const int num_rows = bs0 * dmap0.size();
270 const int num_cols = bs1 * dmap1.size();
271
272 const T* coeff_array = coeffs.data() + index / 2 * cstride;
273 Ae.resize(num_rows * num_cols);
274 std::fill(Ae.begin(), Ae.end(), 0);
275 kernel(Ae.data(), coeff_array, constants.data(), coordinate_dofs.data(),
276 &local_facet, nullptr);
277 dof_transform(Ae, cell_info, cell, num_cols);
278 dof_transform_to_transpose(Ae, cell_info, cell, num_rows);
279
280 // Size data structure for assembly
281 be.resize(num_rows);
282 std::fill(be.begin(), be.end(), 0);
283 for (std::size_t j = 0; j < dmap1.size(); ++j)
284 {
285 for (int k = 0; k < bs1; ++k)
286 {
287 const std::int32_t jj = bs1 * dmap1[j] + k;
288 if (bc_markers1[jj])
289 {
290 const T bc = bc_values1[jj];
291 const T _x0 = x0.empty() ? 0.0 : x0[jj];
292 // be -= Ae.col(bs1 * j + k) * scale * (bc - _x0);
293 for (int m = 0; m < num_rows; ++m)
294 be[m] -= Ae[m * num_cols + bs1 * j + k] * _scale * (bc - _x0);
295 }
296 }
297 }
298
299 for (std::size_t i = 0; i < dmap0.size(); ++i)
300 for (int k = 0; k < bs0; ++k)
301 b[bs0 * dmap0[i] + k] += be[bs0 * i + k];
302 }
303}
304
311template <typename T, int _bs = -1>
312void _lift_bc_interior_facets(
313 std::span<T> b, const mesh::Mesh& mesh,
314 const std::function<void(T*, const T*, const T*,
315 const scalar_value_type_t<T>*, const int*,
316 const std::uint8_t*)>& kernel,
317 const std::span<const std::int32_t>& facets,
318 const std::function<void(const std::span<T>&,
319 const std::span<const std::uint32_t>&,
320 std::int32_t, int)>& dof_transform,
321 const graph::AdjacencyList<std::int32_t>& dofmap0, int bs0,
322 const std::function<void(const std::span<T>&,
323 const std::span<const std::uint32_t>&,
324 std::int32_t, int)>& dof_transform_to_transpose,
325 const graph::AdjacencyList<std::int32_t>& dofmap1, int bs1,
326 const std::span<const T>& constants, const std::span<const T>& coeffs,
327 int cstride, const std::span<const std::uint32_t>& cell_info,
328 const std::function<std::uint8_t(std::size_t)>& get_perm,
329 const std::span<const T>& bc_values1,
330 const std::span<const std::int8_t>& bc_markers1,
331 const std::span<const T>& x0, double scale)
332{
333 if (facets.empty())
334 return;
335
336 const int tdim = mesh.topology().dim();
337
338 // Prepare cell geometry
339 const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
340 const std::size_t num_dofs_g = mesh.geometry().cmap().dim();
341 std::span<const double> x_g = mesh.geometry().x();
342
343 const int num_cell_facets
344 = mesh::cell_num_entities(mesh.topology().cell_type(), tdim - 1);
345
346 // Data structures used in assembly
347 using X = scalar_value_type_t<T>;
348 std::vector<X> coordinate_dofs(2 * num_dofs_g * 3);
349 std::span<X> cdofs0(coordinate_dofs.data(), num_dofs_g * 3);
350 std::span<X> cdofs1(coordinate_dofs.data() + num_dofs_g * 3, num_dofs_g * 3);
351 std::vector<T> Ae, be;
352
353 // Temporaries for joint dofmaps
354 std::vector<std::int32_t> dmapjoint0, dmapjoint1;
355 assert(facets.size() % 4 == 0);
356
357 const scalar_value_type_t<T> _scale
358 = static_cast<scalar_value_type_t<T>>(scale);
359 for (std::size_t index = 0; index < facets.size(); index += 4)
360 {
361 std::array<std::int32_t, 2> cells = {facets[index], facets[index + 2]};
362 std::array<std::int32_t, 2> local_facet
363 = {facets[index + 1], facets[index + 3]};
364
365 // Get cell geometry
366 auto x_dofs0 = x_dofmap.links(cells[0]);
367 for (std::size_t i = 0; i < x_dofs0.size(); ++i)
368 {
369 common::impl::copy_N<3>(std::next(x_g.begin(), 3 * x_dofs0[i]),
370 std::next(cdofs0.begin(), 3 * i));
371 }
372 auto x_dofs1 = x_dofmap.links(cells[1]);
373 for (std::size_t i = 0; i < x_dofs1.size(); ++i)
374 {
375 common::impl::copy_N<3>(std::next(x_g.begin(), 3 * x_dofs1[i]),
376 std::next(cdofs1.begin(), 3 * i));
377 }
378
379 // Get dof maps for cells and pack
380 const std::span<const std::int32_t> dmap0_cell0 = dofmap0.links(cells[0]);
381 const std::span<const std::int32_t> dmap0_cell1 = dofmap0.links(cells[1]);
382 dmapjoint0.resize(dmap0_cell0.size() + dmap0_cell1.size());
383 std::copy(dmap0_cell0.begin(), dmap0_cell0.end(), dmapjoint0.begin());
384 std::copy(dmap0_cell1.begin(), dmap0_cell1.end(),
385 std::next(dmapjoint0.begin(), dmap0_cell0.size()));
386
387 const std::span<const std::int32_t> dmap1_cell0 = dofmap1.links(cells[0]);
388 const std::span<const std::int32_t> dmap1_cell1 = dofmap1.links(cells[1]);
389 dmapjoint1.resize(dmap1_cell0.size() + dmap1_cell1.size());
390 std::copy(dmap1_cell0.begin(), dmap1_cell0.end(), dmapjoint1.begin());
391 std::copy(dmap1_cell1.begin(), dmap1_cell1.end(),
392 std::next(dmapjoint1.begin(), dmap1_cell0.size()));
393
394 // Check if bc is applied to cell0
395 bool has_bc = false;
396 for (std::size_t j = 0; j < dmap1_cell0.size(); ++j)
397 {
398 for (int k = 0; k < bs1; ++k)
399 {
400 if (bc_markers1[bs1 * dmap1_cell0[j] + k])
401 {
402 has_bc = true;
403 break;
404 }
405 }
406 }
407
408 // Check if bc is applied to cell1
409 for (std::size_t j = 0; j < dmap1_cell1.size(); ++j)
410 {
411 for (int k = 0; k < bs1; ++k)
412 {
413 if (bc_markers1[bs1 * dmap1_cell1[j] + k])
414 {
415 has_bc = true;
416 break;
417 }
418 }
419 }
420
421 if (!has_bc)
422 continue;
423
424 const int num_rows = bs0 * dmapjoint0.size();
425 const int num_cols = bs1 * dmapjoint1.size();
426
427 // Tabulate tensor
428 Ae.resize(num_rows * num_cols);
429 std::fill(Ae.begin(), Ae.end(), 0);
430 const std::array perm{
431 get_perm(cells[0] * num_cell_facets + local_facet[0]),
432 get_perm(cells[1] * num_cell_facets + local_facet[1])};
433 kernel(Ae.data(), coeffs.data() + index / 2 * cstride, constants.data(),
434 coordinate_dofs.data(), local_facet.data(), perm.data());
435
436 const std::span<T> _Ae(Ae);
437 const std::span<T> sub_Ae0
438 = _Ae.subspan(bs0 * dmap0_cell0.size() * num_cols,
439 bs0 * dmap0_cell1.size() * num_cols);
440 const std::span<T> sub_Ae1
441 = _Ae.subspan(bs1 * dmap1_cell0.size(),
442 num_rows * num_cols - bs1 * dmap1_cell0.size());
443
444 dof_transform(_Ae, cell_info, cells[0], num_cols);
445 dof_transform(sub_Ae0, cell_info, cells[1], num_cols);
446 dof_transform_to_transpose(_Ae, cell_info, cells[0], num_rows);
447 dof_transform_to_transpose(sub_Ae1, cell_info, cells[1], num_rows);
448
449 be.resize(num_rows);
450 std::fill(be.begin(), be.end(), 0);
451
452 // Compute b = b - A*b for cell0
453 for (std::size_t j = 0; j < dmap1_cell0.size(); ++j)
454 {
455 for (int k = 0; k < bs1; ++k)
456 {
457 const std::int32_t jj = bs1 * dmap1_cell0[j] + k;
458 if (bc_markers1[jj])
459 {
460 const T bc = bc_values1[jj];
461 const T _x0 = x0.empty() ? 0.0 : x0[jj];
462 // be -= Ae.col(bs1 * j + k) * scale * (bc - _x0);
463 for (int m = 0; m < num_rows; ++m)
464 be[m] -= Ae[m * num_cols + bs1 * j + k] * _scale * (bc - _x0);
465 }
466 }
467 }
468
469 // Compute b = b - A*b for cell1
470 const int offset = bs1 * dmap1_cell0.size();
471 for (std::size_t j = 0; j < dmap1_cell1.size(); ++j)
472 {
473 for (int k = 0; k < bs1; ++k)
474 {
475 const std::int32_t jj = bs1 * dmap1_cell1[j] + k;
476 if (bc_markers1[jj])
477 {
478 const T bc = bc_values1[jj];
479 const T _x0 = x0.empty() ? 0.0 : x0[jj];
480 // be -= Ae.col(offset + bs1 * j + k) * scale * (bc - x0[jj]);
481 for (int m = 0; m < num_rows; ++m)
482 {
483 be[m] -= Ae[m * num_cols + offset + bs1 * j + k] * _scale
484 * (bc - _x0);
485 }
486 }
487 }
488 }
489
490 for (std::size_t i = 0; i < dmap0_cell0.size(); ++i)
491 for (int k = 0; k < bs0; ++k)
492 b[bs0 * dmap0_cell0[i] + k] += be[bs0 * i + k];
493
494 const int offset_be = bs0 * dmap0_cell0.size();
495 for (std::size_t i = 0; i < dmap0_cell1.size(); ++i)
496 for (int k = 0; k < bs0; ++k)
497 b[bs0 * dmap0_cell1[i] + k] += be[offset_be + bs0 * i + k];
498 }
499}
506template <typename T, int _bs = -1>
507void assemble_cells(
508 const std::function<void(const std::span<T>&,
509 const std::span<const std::uint32_t>&,
510 std::int32_t, int)>& dof_transform,
511 std::span<T> b, const mesh::Geometry& geometry,
512 const std::span<const std::int32_t>& cells,
513 const graph::AdjacencyList<std::int32_t>& dofmap, int bs,
514 const std::function<void(T*, const T*, const T*,
515 const scalar_value_type_t<T>*, const int*,
516 const std::uint8_t*)>& kernel,
517 const std::span<const T>& constants, const std::span<const T>& coeffs,
518 int cstride, const std::span<const std::uint32_t>& cell_info)
519{
520 assert(_bs < 0 or _bs == bs);
521
522 if (cells.empty())
523 return;
524
525 // Prepare cell geometry
526 const graph::AdjacencyList<std::int32_t>& x_dofmap = geometry.dofmap();
527 const std::size_t num_dofs_g = geometry.cmap().dim();
528 std::span<const double> x_g = geometry.x();
529
530 // FIXME: Add proper interface for num_dofs
531 // Create data structures used in assembly
532 const int num_dofs = dofmap.links(0).size();
533 std::vector<scalar_value_type_t<T>> coordinate_dofs(3 * num_dofs_g);
534 std::vector<T> be(bs * num_dofs);
535 const std::span<T> _be(be);
536
537 // Iterate over active cells
538 for (std::size_t index = 0; index < cells.size(); ++index)
539 {
540 std::int32_t c = cells[index];
541
542 // Get cell coordinates/geometry
543 auto x_dofs = x_dofmap.links(c);
544 for (std::size_t i = 0; i < x_dofs.size(); ++i)
545 {
546 common::impl::copy_N<3>(std::next(x_g.begin(), 3 * x_dofs[i]),
547 std::next(coordinate_dofs.begin(), 3 * i));
548 }
549
550 // Tabulate vector for cell
551 std::fill(be.begin(), be.end(), 0);
552 kernel(be.data(), coeffs.data() + index * cstride, constants.data(),
553 coordinate_dofs.data(), nullptr, nullptr);
554 dof_transform(_be, cell_info, c, 1);
555
556 // Scatter cell vector to 'global' vector array
557 auto dofs = dofmap.links(c);
558 if constexpr (_bs > 0)
559 {
560 for (int i = 0; i < num_dofs; ++i)
561 for (int k = 0; k < _bs; ++k)
562 b[_bs * dofs[i] + k] += be[_bs * i + k];
563 }
564 else
565 {
566 for (int i = 0; i < num_dofs; ++i)
567 for (int k = 0; k < bs; ++k)
568 b[bs * dofs[i] + k] += be[bs * i + k];
569 }
570 }
571}
572
579template <typename T, int _bs = -1>
580void assemble_exterior_facets(
581 const std::function<void(const std::span<T>&,
582 const std::span<const std::uint32_t>&,
583 std::int32_t, int)>& dof_transform,
584 std::span<T> b, const mesh::Mesh& mesh,
585 const std::span<const std::int32_t>& facets,
586 const graph::AdjacencyList<std::int32_t>& dofmap, int bs,
587 const std::function<void(T*, const T*, const T*,
588 const scalar_value_type_t<T>*, const int*,
589 const std::uint8_t*)>& fn,
590 const std::span<const T>& constants, const std::span<const T>& coeffs,
591 int cstride, const std::span<const std::uint32_t>& cell_info)
592{
593 assert(_bs < 0 or _bs == bs);
594
595 if (facets.empty())
596 return;
597
598 // Prepare cell geometry
599 const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
600 const std::size_t num_dofs_g = mesh.geometry().cmap().dim();
601 std::span<const double> x_g = mesh.geometry().x();
602
603 // FIXME: Add proper interface for num_dofs
604 // Create data structures used in assembly
605 const int num_dofs = dofmap.links(0).size();
606 std::vector<scalar_value_type_t<T>> coordinate_dofs(3 * num_dofs_g);
607 std::vector<T> be(bs * num_dofs);
608 const std::span<T> _be(be);
609 assert(facets.size() % 2 == 0);
610 for (std::size_t index = 0; index < facets.size(); index += 2)
611 {
612 std::int32_t cell = facets[index];
613 std::int32_t local_facet = facets[index + 1];
614
615 // Get cell coordinates/geometry
616 auto x_dofs = x_dofmap.links(cell);
617 for (std::size_t i = 0; i < x_dofs.size(); ++i)
618 {
619 common::impl::copy_N<3>(std::next(x_g.begin(), 3 * x_dofs[i]),
620 std::next(coordinate_dofs.begin(), 3 * i));
621 }
622
623 // Tabulate element vector
624 std::fill(be.begin(), be.end(), 0);
625 fn(be.data(), coeffs.data() + index / 2 * cstride, constants.data(),
626 coordinate_dofs.data(), &local_facet, nullptr);
627
628 dof_transform(_be, cell_info, cell, 1);
629
630 // Add element vector to global vector
631 auto dofs = dofmap.links(cell);
632 if constexpr (_bs > 0)
633 {
634 for (int i = 0; i < num_dofs; ++i)
635 for (int k = 0; k < _bs; ++k)
636 b[_bs * dofs[i] + k] += be[_bs * i + k];
637 }
638 else
639 {
640 for (int i = 0; i < num_dofs; ++i)
641 for (int k = 0; k < bs; ++k)
642 b[bs * dofs[i] + k] += be[bs * i + k];
643 }
644 }
645}
646
653template <typename T, int _bs = -1>
654void assemble_interior_facets(
655 const std::function<void(const std::span<T>&,
656 const std::span<const std::uint32_t>&,
657 std::int32_t, int)>& dof_transform,
658 std::span<T> b, const mesh::Mesh& mesh,
659 const std::span<const std::int32_t>& facets, const fem::DofMap& dofmap,
660 const std::function<void(T*, const T*, const T*,
661 const scalar_value_type_t<T>*, const int*,
662 const std::uint8_t*)>& fn,
663 const std::span<const T>& constants, const std::span<const T>& coeffs,
664 int cstride, const std::span<const std::uint32_t>& cell_info,
665 const std::function<std::uint8_t(std::size_t)>& get_perm)
666{
667 const int tdim = mesh.topology().dim();
668
669 // Prepare cell geometry
670 const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
671 const std::size_t num_dofs_g = mesh.geometry().cmap().dim();
672 std::span<const double> x_g = mesh.geometry().x();
673
674 // Create data structures used in assembly
675 using X = scalar_value_type_t<T>;
676 std::vector<X> coordinate_dofs(2 * num_dofs_g * 3);
677 std::span<X> cdofs0(coordinate_dofs.data(), num_dofs_g * 3);
678 std::span<X> cdofs1(coordinate_dofs.data() + num_dofs_g * 3, num_dofs_g * 3);
679 std::vector<T> be;
680
681 const int num_cell_facets
682 = mesh::cell_num_entities(mesh.topology().cell_type(), tdim - 1);
683
684 const int bs = dofmap.bs();
685 assert(_bs < 0 or _bs == bs);
686 assert(facets.size() % 4 == 0);
687 for (std::size_t index = 0; index < facets.size(); index += 4)
688 {
689 std::array<std::int32_t, 2> cells = {facets[index], facets[index + 2]};
690 std::array<std::int32_t, 2> local_facet
691 = {facets[index + 1], facets[index + 3]};
692
693 // Get cell geometry
694 auto x_dofs0 = x_dofmap.links(cells[0]);
695 for (std::size_t i = 0; i < x_dofs0.size(); ++i)
696 {
697 common::impl::copy_N<3>(std::next(x_g.begin(), 3 * x_dofs0[i]),
698 std::next(cdofs0.begin(), 3 * i));
699 }
700 auto x_dofs1 = x_dofmap.links(cells[1]);
701 for (std::size_t i = 0; i < x_dofs1.size(); ++i)
702 {
703 common::impl::copy_N<3>(std::next(x_g.begin(), 3 * x_dofs1[i]),
704 std::next(cdofs1.begin(), 3 * i));
705 }
706
707 // Get dofmaps for cells
708 std::span<const std::int32_t> dmap0 = dofmap.cell_dofs(cells[0]);
709 std::span<const std::int32_t> dmap1 = dofmap.cell_dofs(cells[1]);
710
711 // Tabulate element vector
712 be.resize(bs * (dmap0.size() + dmap1.size()));
713 std::fill(be.begin(), be.end(), 0);
714 const std::array perm{
715 get_perm(cells[0] * num_cell_facets + local_facet[0]),
716 get_perm(cells[1] * num_cell_facets + local_facet[1])};
717 fn(be.data(), coeffs.data() + index / 2 * cstride, constants.data(),
718 coordinate_dofs.data(), local_facet.data(), perm.data());
719
720 const std::span<T> _be(be);
721 const std::span<T> sub_be
722 = _be.subspan(bs * dmap0.size(), bs * dmap1.size());
723
724 dof_transform(be, cell_info, cells[0], 1);
725 dof_transform(sub_be, cell_info, cells[1], 1);
726
727 // Add element vector to global vector
728 if constexpr (_bs > 0)
729 {
730 for (std::size_t i = 0; i < dmap0.size(); ++i)
731 for (int k = 0; k < _bs; ++k)
732 b[_bs * dmap0[i] + k] += be[_bs * i + k];
733 for (std::size_t i = 0; i < dmap1.size(); ++i)
734 for (int k = 0; k < _bs; ++k)
735 b[_bs * dmap1[i] + k] += be[_bs * (i + dmap0.size()) + k];
736 }
737 else
738 {
739 for (std::size_t i = 0; i < dmap0.size(); ++i)
740 for (int k = 0; k < bs; ++k)
741 b[bs * dmap0[i] + k] += be[bs * i + k];
742 for (std::size_t i = 0; i < dmap1.size(); ++i)
743 for (int k = 0; k < bs; ++k)
744 b[bs * dmap1[i] + k] += be[bs * (i + dmap0.size()) + k];
745 }
746 }
747}
748
763template <typename T>
764void lift_bc(std::span<T> b, const Form<T>& a,
765 const std::span<const T>& constants,
766 const std::map<std::pair<IntegralType, int>,
767 std::pair<std::span<const T>, int>>& coefficients,
768 const std::span<const T>& bc_values1,
769 const std::span<const std::int8_t>& bc_markers1,
770 const std::span<const T>& x0, double scale)
771{
772 std::shared_ptr<const mesh::Mesh> mesh = a.mesh();
773 assert(mesh);
774
775 // Get dofmap for columns and rows of a
776 assert(a.function_spaces().at(0));
777 assert(a.function_spaces().at(1));
778 const graph::AdjacencyList<std::int32_t>& dofmap0
779 = a.function_spaces()[0]->dofmap()->list();
780 const int bs0 = a.function_spaces()[0]->dofmap()->bs();
781 std::shared_ptr<const fem::FiniteElement> element0
782 = a.function_spaces()[0]->element();
783 const graph::AdjacencyList<std::int32_t>& dofmap1
784 = a.function_spaces()[1]->dofmap()->list();
785 const int bs1 = a.function_spaces()[1]->dofmap()->bs();
786 std::shared_ptr<const fem::FiniteElement> element1
787 = a.function_spaces()[1]->element();
788
789 const bool needs_transformation_data
790 = element0->needs_dof_transformations()
791 or element1->needs_dof_transformations()
792 or a.needs_facet_permutations();
793
794 std::span<const std::uint32_t> cell_info;
795 if (needs_transformation_data)
796 {
797 mesh->topology_mutable().create_entity_permutations();
798 cell_info = std::span(mesh->topology().get_cell_permutation_info());
799 }
800
801 const std::function<void(const std::span<T>&,
802 const std::span<const std::uint32_t>&, std::int32_t,
803 int)>
804 dof_transform = element0->get_dof_transformation_function<T>();
805 const std::function<void(const std::span<T>&,
806 const std::span<const std::uint32_t>&, std::int32_t,
807 int)>
808 dof_transform_to_transpose
809 = element1->get_dof_transformation_to_transpose_function<T>();
810
811 for (int i : a.integral_ids(IntegralType::cell))
812 {
813 const auto& kernel = a.kernel(IntegralType::cell, i);
814 const auto& [coeffs, cstride] = coefficients.at({IntegralType::cell, i});
815 const std::vector<std::int32_t>& cells = a.cell_domains(i);
816 if (bs0 == 1 and bs1 == 1)
817 {
818 _lift_bc_cells<T, 1, 1>(b, mesh->geometry(), kernel, cells, dof_transform,
819 dofmap0, bs0, dof_transform_to_transpose, dofmap1,
820 bs1, constants, coeffs, cstride, cell_info,
821 bc_values1, bc_markers1, x0, scale);
822 }
823 else if (bs0 == 3 and bs1 == 3)
824 {
825 _lift_bc_cells<T, 3, 3>(b, mesh->geometry(), kernel, cells, dof_transform,
826 dofmap0, bs0, dof_transform_to_transpose, dofmap1,
827 bs1, constants, coeffs, cstride, cell_info,
828 bc_values1, bc_markers1, x0, scale);
829 }
830 else
831 {
832 _lift_bc_cells(b, mesh->geometry(), kernel, cells, dof_transform, dofmap0,
833 bs0, dof_transform_to_transpose, dofmap1, bs1, constants,
834 coeffs, cstride, cell_info, bc_values1, bc_markers1, x0,
835 scale);
836 }
837 }
838
839 for (int i : a.integral_ids(IntegralType::exterior_facet))
840 {
841 const auto& kernel = a.kernel(IntegralType::exterior_facet, i);
842 const auto& [coeffs, cstride]
843 = coefficients.at({IntegralType::exterior_facet, i});
844 const std::vector<std::int32_t>& facets = a.exterior_facet_domains(i);
845 _lift_bc_exterior_facets(b, *mesh, kernel, facets, dof_transform, dofmap0,
846 bs0, dof_transform_to_transpose, dofmap1, bs1,
847 constants, coeffs, cstride, cell_info, bc_values1,
848 bc_markers1, x0, scale);
849 }
850
851 if (a.num_integrals(IntegralType::interior_facet) > 0)
852 {
853 std::function<std::uint8_t(std::size_t)> get_perm;
854 if (a.needs_facet_permutations())
855 {
856 mesh->topology_mutable().create_entity_permutations();
857 const std::vector<std::uint8_t>& perms
858 = mesh->topology().get_facet_permutations();
859 get_perm = [&perms](std::size_t i) { return perms[i]; };
860 }
861 else
862 get_perm = [](std::size_t) { return 0; };
863
864 for (int i : a.integral_ids(IntegralType::interior_facet))
865 {
866 const auto& kernel = a.kernel(IntegralType::interior_facet, i);
867 const auto& [coeffs, cstride]
868 = coefficients.at({IntegralType::interior_facet, i});
869 const std::vector<std::int32_t>& facets = a.interior_facet_domains(i);
870 _lift_bc_interior_facets(b, *mesh, kernel, facets, dof_transform, dofmap0,
871 bs0, dof_transform_to_transpose, dofmap1, bs1,
872 constants, coeffs, cstride, cell_info, get_perm,
873 bc_values1, bc_markers1, x0, scale);
874 }
875 }
876}
877
897template <typename T>
898void apply_lifting(
899 std::span<T> b, const std::vector<std::shared_ptr<const Form<T>>> a,
900 const std::vector<std::span<const T>>& constants,
901 const std::vector<std::map<std::pair<IntegralType, int>,
902 std::pair<std::span<const T>, int>>>& coeffs,
903 const std::vector<std::vector<std::shared_ptr<const DirichletBC<T>>>>& bcs1,
904 const std::vector<std::span<const T>>& x0, double scale)
905{
906 // FIXME: make changes to reactivate this check
907 if (!x0.empty() and x0.size() != a.size())
908 {
909 throw std::runtime_error(
910 "Mismatch in size between x0 and bilinear form in assembler.");
911 }
912
913 if (a.size() != bcs1.size())
914 {
915 throw std::runtime_error(
916 "Mismatch in size between a and bcs in assembler.");
917 }
918
919 for (std::size_t j = 0; j < a.size(); ++j)
920 {
921 std::vector<std::int8_t> bc_markers1;
922 std::vector<T> bc_values1;
923 if (a[j] and !bcs1[j].empty())
924 {
925 assert(a[j]->function_spaces().at(0));
926
927 auto V1 = a[j]->function_spaces()[1];
928 assert(V1);
929 auto map1 = V1->dofmap()->index_map;
930 const int bs1 = V1->dofmap()->index_map_bs();
931 assert(map1);
932 const int crange = bs1 * (map1->size_local() + map1->num_ghosts());
933 bc_markers1.assign(crange, false);
934 bc_values1.assign(crange, 0.0);
935 for (const std::shared_ptr<const DirichletBC<T>>& bc : bcs1[j])
936 {
937 bc->mark_dofs(bc_markers1);
938 bc->dof_values(bc_values1);
939 }
940
941 if (!x0.empty())
942 {
943 lift_bc<T>(b, *a[j], constants[j], coeffs[j], bc_values1, bc_markers1,
944 x0[j], scale);
945 }
946 else
947 {
948 lift_bc<T>(b, *a[j], constants[j], coeffs[j], bc_values1, bc_markers1,
949 std::span<const T>(), scale);
950 }
951 }
952 }
953}
954
961template <typename T>
962void assemble_vector(
963 std::span<T> b, const Form<T>& L, const std::span<const T>& constants,
964 const std::map<std::pair<IntegralType, int>,
965 std::pair<std::span<const T>, int>>& coefficients)
966{
967 std::shared_ptr<const mesh::Mesh> mesh = L.mesh();
968 assert(mesh);
969
970 // Get dofmap data
971 assert(L.function_spaces().at(0));
972 std::shared_ptr<const fem::FiniteElement> element
973 = L.function_spaces().at(0)->element();
974 std::shared_ptr<const fem::DofMap> dofmap
975 = L.function_spaces().at(0)->dofmap();
976 assert(dofmap);
977 const graph::AdjacencyList<std::int32_t>& dofs = dofmap->list();
978 const int bs = dofmap->bs();
979
980 const std::function<void(const std::span<T>&,
981 const std::span<const std::uint32_t>&, std::int32_t,
982 int)>
983 dof_transform = element->get_dof_transformation_function<T>();
984
985 const bool needs_transformation_data
986 = element->needs_dof_transformations() or L.needs_facet_permutations();
987 std::span<const std::uint32_t> cell_info;
988 if (needs_transformation_data)
989 {
990 mesh->topology_mutable().create_entity_permutations();
991 cell_info = std::span(mesh->topology().get_cell_permutation_info());
992 }
993
994 for (int i : L.integral_ids(IntegralType::cell))
995 {
996 const auto& fn = L.kernel(IntegralType::cell, i);
997 const auto& [coeffs, cstride] = coefficients.at({IntegralType::cell, i});
998 const std::vector<std::int32_t>& cells = L.cell_domains(i);
999 if (bs == 1)
1000 {
1001 impl::assemble_cells<T, 1>(dof_transform, b, mesh->geometry(), cells,
1002 dofs, bs, fn, constants, coeffs, cstride,
1003 cell_info);
1004 }
1005 else if (bs == 3)
1006 {
1007 impl::assemble_cells<T, 3>(dof_transform, b, mesh->geometry(), cells,
1008 dofs, bs, fn, constants, coeffs, cstride,
1009 cell_info);
1010 }
1011 else
1012 {
1013 impl::assemble_cells(dof_transform, b, mesh->geometry(), cells, dofs, bs,
1014 fn, constants, coeffs, cstride, cell_info);
1015 }
1016 }
1017
1018 for (int i : L.integral_ids(IntegralType::exterior_facet))
1019 {
1020 const auto& fn = L.kernel(IntegralType::exterior_facet, i);
1021 const auto& [coeffs, cstride]
1022 = coefficients.at({IntegralType::exterior_facet, i});
1023 const std::vector<std::int32_t>& facets = L.exterior_facet_domains(i);
1024 if (bs == 1)
1025 {
1026 impl::assemble_exterior_facets<T, 1>(dof_transform, b, *mesh, facets,
1027 dofs, bs, fn, constants, coeffs,
1028 cstride, cell_info);
1029 }
1030 else if (bs == 3)
1031 {
1032 impl::assemble_exterior_facets<T, 3>(dof_transform, b, *mesh, facets,
1033 dofs, bs, fn, constants, coeffs,
1034 cstride, cell_info);
1035 }
1036 else
1037 {
1038 impl::assemble_exterior_facets(dof_transform, b, *mesh, facets, dofs, bs,
1039 fn, constants, coeffs, cstride, cell_info);
1040 }
1041 }
1042
1043 if (L.num_integrals(IntegralType::interior_facet) > 0)
1044 {
1045 std::function<std::uint8_t(std::size_t)> get_perm;
1046 if (L.needs_facet_permutations())
1047 {
1048 mesh->topology_mutable().create_entity_permutations();
1049 const std::vector<std::uint8_t>& perms
1050 = mesh->topology().get_facet_permutations();
1051 get_perm = [&perms](std::size_t i) { return perms[i]; };
1052 }
1053 else
1054 get_perm = [](std::size_t) { return 0; };
1055
1056 for (int i : L.integral_ids(IntegralType::interior_facet))
1057 {
1058 const auto& fn = L.kernel(IntegralType::interior_facet, i);
1059 const auto& [coeffs, cstride]
1060 = coefficients.at({IntegralType::interior_facet, i});
1061 const std::vector<std::int32_t>& facets = L.interior_facet_domains(i);
1062 if (bs == 1)
1063 {
1064 impl::assemble_interior_facets<T, 1>(dof_transform, b, *mesh, facets,
1065 *dofmap, fn, constants, coeffs,
1066 cstride, cell_info, get_perm);
1067 }
1068 else if (bs == 3)
1069 {
1070 impl::assemble_interior_facets<T, 3>(dof_transform, b, *mesh, facets,
1071 *dofmap, fn, constants, coeffs,
1072 cstride, cell_info, get_perm);
1073 }
1074 else
1075 {
1076 impl::assemble_interior_facets(dof_transform, b, *mesh, facets, *dofmap,
1077 fn, constants, coeffs, cstride,
1078 cell_info, get_perm);
1079 }
1080 }
1081 }
1082}
1083} // namespace dolfinx::fem::impl
Degree-of-freedeom map representations ans tools.
void cells(la::SparsityPattern &pattern, const mesh::Topology &topology, const std::array< const std::reference_wrapper< const DofMap >, 2 > &dofmaps)
Iterate over cells and insert entries into sparsity pattern.
Definition: sparsitybuild.cpp:18
void apply_lifting(std::span< T > b, const std::vector< std::shared_ptr< const Form< T > > > &a, const std::vector< std::span< const T > > &constants, const std::vector< std::map< std::pair< IntegralType, int >, std::pair< std::span< const T >, int > > > &coeffs, const std::vector< std::vector< std::shared_ptr< const DirichletBC< T > > > > &bcs1, const std::vector< std::span< const T > > &x0, double scale)
Modify b such that:
Definition: assembler.h:130
void assemble_vector(std::span< T > b, const Form< T > &L, const std::span< const T > &constants, const std::map< std::pair< IntegralType, int >, std::pair< std::span< const T >, int > > &coefficients)
Assemble linear form into a vector, The caller supplies the form constants and coefficients for this ...
Definition: assembler.h:89
@ 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