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