9 #include "DirichletBC.h"
13 #include <dolfinx/common/IndexMap.h>
14 #include <dolfinx/common/types.h>
15 #include <dolfinx/fem/Constant.h>
16 #include <dolfinx/fem/FunctionSpace.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>
25 namespace dolfinx::fem::impl
40 tcb::span<T> b,
const mesh::Geometry& geometry,
41 const std::vector<std::int32_t>& active_cells,
42 const graph::AdjacencyList<std::int32_t>& dofmap,
const int bs,
43 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
44 const std::uint8_t*,
const std::uint32_t)>& kernel,
45 const array2d<T>& coeffs,
const std::vector<T>& constant_values,
46 const std::vector<std::uint32_t>& cell_info);
50 void assemble_exterior_facets(
51 tcb::span<T> b,
const mesh::Mesh& mesh,
52 const std::vector<std::int32_t>& active_facets,
53 const graph::AdjacencyList<std::int32_t>& dofmap,
const int bs,
54 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
55 const std::uint8_t*,
const std::uint32_t)>& fn,
56 const array2d<T>& coeffs,
const std::vector<T>& constant_values,
57 const std::vector<std::uint32_t>& cell_info,
58 const std::vector<std::uint8_t>& perms);
62 void assemble_interior_facets(
63 tcb::span<T> b,
const mesh::Mesh& mesh,
64 const std::vector<std::int32_t>& active_facets,
const fem::DofMap& dofmap,
65 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
66 const std::uint8_t*,
const std::uint32_t)>& fn,
67 const array2d<T>& coeffs,
const std::vector<int>& offsets,
68 const std::vector<T>& constant_values,
69 const std::vector<std::uint32_t>& cell_info,
70 const std::vector<std::uint8_t>& perms);
91 tcb::span<T> b,
const std::vector<std::shared_ptr<
const Form<T>>> a,
92 const std::vector<std::vector<std::shared_ptr<
const DirichletBC<T>>>>& bcs1,
93 const tcb::span<const T>& x0,
double scale);
107 template <
typename T>
108 void lift_bc(tcb::span<T> b,
const Form<T>& a,
109 const tcb::span<const T>& bc_values1,
110 const std::vector<bool>& bc_markers1,
const tcb::span<const T>& x0,
114 template <
typename T>
116 tcb::span<T> b,
const mesh::Geometry& geometry,
117 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
118 const std::uint8_t*,
const std::uint32_t)>& kernel,
119 const std::vector<std::int32_t>& active_cells,
120 const graph::AdjacencyList<std::int32_t>& dofmap0,
int bs0,
121 const graph::AdjacencyList<std::int32_t>& dofmap1,
int bs1,
122 const array2d<T>& coeffs,
const std::vector<T>& constant_values,
123 const std::vector<std::uint32_t>& cell_info,
124 const tcb::span<const T>& bc_values1,
const std::vector<bool>& bc_markers1,
125 const tcb::span<const T>& x0,
double scale)
128 const int gdim = geometry.dim();
129 const graph::AdjacencyList<std::int32_t>& x_dofmap = geometry.dofmap();
132 const int num_dofs_g = x_dofmap.num_links(0);
133 const array2d<double>& x_g = geometry.x();
136 std::vector<double> coordinate_dofs(num_dofs_g * gdim);
137 std::vector<T> Ae, be;
139 for (std::int32_t c : active_cells)
142 auto dmap1 = dofmap1.links(c);
146 for (std::size_t j = 0; j < dmap1.size(); ++j)
148 for (
int k = 0; k < bs1; ++k)
150 assert(bs1 * dmap1[j] + k < (
int)bc_markers1.size());
151 if (bc_markers1[bs1 * dmap1[j] + k])
163 auto x_dofs = x_dofmap.links(c);
164 for (std::size_t i = 0; i < x_dofs.size(); ++i)
166 std::copy_n(x_g.row(x_dofs[i]).data(), gdim,
167 std::next(coordinate_dofs.begin(), i * gdim));
171 auto dmap0 = dofmap0.links(c);
173 const int num_rows = bs0 * dmap0.size();
174 const int num_cols = bs1 * dmap1.size();
176 auto coeff_array = coeffs.row(c);
177 Ae.resize(num_rows * num_cols);
178 std::fill(Ae.begin(), Ae.end(), 0);
179 kernel(Ae.data(), coeff_array.data(), constant_values.data(),
180 coordinate_dofs.data(),
nullptr,
nullptr, cell_info[c]);
184 std::fill(be.begin(), be.end(), 0);
185 for (std::size_t j = 0; j < dmap1.size(); ++j)
187 for (
int k = 0; k < bs1; ++k)
189 const std::int32_t jj = bs1 * dmap1[j] + k;
190 assert(jj < (
int)bc_markers1.size());
193 const T bc = bc_values1[jj];
194 const T _x0 = x0.empty() ? 0.0 : x0[jj];
196 for (
int m = 0; m < num_rows; ++m)
197 be[m] -= Ae[m * num_cols + bs1 * j + k] * scale * (bc - _x0);
202 for (std::size_t i = 0; i < dmap0.size(); ++i)
203 for (
int k = 0; k < bs0; ++k)
204 b[bs0 * dmap0[i] + k] += be[bs0 * i + k];
208 template <
typename T>
209 void _lift_bc_exterior_facets(
210 tcb::span<T> b,
const mesh::Mesh& mesh,
211 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
212 const std::uint8_t*,
const std::uint32_t)>& kernel,
213 const std::vector<std::int32_t>& active_facets,
214 const graph::AdjacencyList<std::int32_t>& dofmap0,
int bs0,
215 const graph::AdjacencyList<std::int32_t>& dofmap1,
int bs1,
216 const array2d<T>& coeffs,
const std::vector<T>& constant_values,
217 const std::vector<std::uint32_t>& cell_info,
218 const std::vector<std::uint8_t>& perms,
219 const tcb::span<const T>& bc_values1,
const std::vector<bool>& bc_markers1,
220 const tcb::span<const T>& x0,
double scale)
222 const int gdim = mesh.geometry().dim();
223 const int tdim = mesh.topology().dim();
226 const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
229 const int num_dofs_g = x_dofmap.num_links(0);
230 const array2d<double>& x_g = mesh.geometry().x();
233 std::vector<double> coordinate_dofs(num_dofs_g * gdim);
234 std::vector<T> Ae, be;
237 const mesh::Topology& topology = mesh.topology();
238 auto connectivity = topology.connectivity(tdim - 1, tdim);
239 assert(connectivity);
240 auto c_to_f = topology.connectivity(tdim, tdim - 1);
242 auto map = topology.index_map(tdim - 1);
245 for (std::int32_t f : active_facets)
248 assert(connectivity->num_links(f) == 1);
249 const std::int32_t cell = connectivity->links(f)[0];
252 auto facets = c_to_f->links(cell);
253 auto it = std::find(facets.begin(), facets.end(), f);
254 assert(it != facets.end());
255 const int local_facet = std::distance(facets.begin(), it);
258 auto dmap1 = dofmap1.links(cell);
262 for (std::size_t j = 0; j < dmap1.size(); ++j)
264 for (
int k = 0; k < bs1; ++k)
266 if (bc_markers1[bs1 * dmap1[j] + k])
278 auto x_dofs = x_dofmap.links(cell);
279 for (std::size_t i = 0; i < x_dofs.size(); ++i)
281 std::copy_n(x_g.row(x_dofs[i]).data(), gdim,
282 std::next(coordinate_dofs.begin(), i * gdim));
286 auto dmap0 = dofmap0.links(cell);
288 const int num_rows = bs0 * dmap0.size();
289 const int num_cols = bs1 * dmap1.size();
291 auto coeff_array = coeffs.row(cell);
292 Ae.resize(num_rows * num_cols);
293 std::fill(Ae.begin(), Ae.end(), 0);
294 kernel(Ae.data(), coeff_array.data(), constant_values.data(),
295 coordinate_dofs.data(), &local_facet,
296 &perms[cell * facets.size() + local_facet], cell_info[cell]);
300 std::fill(be.begin(), be.end(), 0);
301 for (std::size_t j = 0; j < dmap1.size(); ++j)
303 for (
int k = 0; k < bs1; ++k)
305 const std::int32_t jj = bs1 * dmap1[j] + k;
309 const T bc = bc_values1[jj];
310 const T _x0 = x0.empty() ? 0.0 : x0[jj];
312 for (
int m = 0; m < num_rows; ++m)
313 be[m] -= Ae[m * num_cols + bs1 * j + k] * scale * (bc - _x0);
318 for (std::size_t i = 0; i < dmap0.size(); ++i)
319 for (
int k = 0; k < bs0; ++k)
320 b[bs0 * dmap0[i] + k] += be[bs0 * i + k];
324 template <
typename T>
325 void _lift_bc_interior_facets(
326 tcb::span<T> b,
const mesh::Mesh& mesh,
327 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
328 const std::uint8_t*,
const std::uint32_t)>& kernel,
329 const std::vector<std::int32_t>& active_facets,
330 const graph::AdjacencyList<std::int32_t>& dofmap0,
int bs0,
331 const graph::AdjacencyList<std::int32_t>& dofmap1,
int bs1,
332 const array2d<T>& coeffs,
const std::vector<int>& offsets,
333 const std::vector<T>& constant_values,
334 const std::vector<std::uint32_t>& cell_info,
335 const std::vector<std::uint8_t>& perms,
336 const tcb::span<const T>& bc_values1,
const std::vector<bool>& bc_markers1,
337 const tcb::span<const T>& x0,
double scale)
339 const int gdim = mesh.geometry().dim();
340 const int tdim = mesh.topology().dim();
343 const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
346 const int num_dofs_g = x_dofmap.num_links(0);
347 const array2d<double>& x_g = mesh.geometry().x();
350 std::vector<double> coordinate_dofs(2 * num_dofs_g * gdim);
351 std::vector<T> coeff_array(2 * offsets.back());
352 assert(offsets.back() ==
int(coeffs.shape[1]));
353 std::vector<T> Ae, be;
356 std::vector<std::int32_t> dmapjoint0, dmapjoint1;
358 const mesh::Topology& topology = mesh.topology();
359 auto connectivity = topology.connectivity(tdim - 1, tdim);
360 assert(connectivity);
361 auto c_to_f = topology.connectivity(tdim, tdim - 1);
363 auto f_to_c = topology.connectivity(tdim - 1, tdim);
365 auto map = topology.index_map(tdim - 1);
368 const int offset_g = gdim * num_dofs_g;
369 for (std::int32_t f : active_facets)
372 auto cells = f_to_c->links(f);
373 assert(
cells.size() == 2);
376 auto facets0 = c_to_f->links(cells[0]);
377 auto it0 = std::find(facets0.begin(), facets0.end(), f);
378 assert(it0 != facets0.end());
379 const int local_facet0 = std::distance(facets0.begin(), it0);
380 auto facets1 = c_to_f->links(cells[1]);
381 auto it1 = std::find(facets1.begin(), facets1.end(), f);
382 assert(it1 != facets1.end());
383 const int local_facet1 = std::distance(facets1.begin(), it1);
385 const std::array local_facet{local_facet0, local_facet1};
388 auto x_dofs0 = x_dofmap.links(cells[0]);
389 auto x_dofs1 = x_dofmap.links(cells[1]);
390 for (
int i = 0; i < num_dofs_g; ++i)
392 for (
int j = 0; j < gdim; ++j)
394 coordinate_dofs[i * gdim + j] = x_g(x_dofs0[i], j);
395 coordinate_dofs[offset_g + i * gdim + j] = x_g(x_dofs1[i], j);
400 const tcb::span<const std::int32_t> dmap0_cell0 = dofmap0.links(cells[0]);
401 const tcb::span<const std::int32_t> dmap0_cell1 = dofmap0.links(cells[1]);
402 dmapjoint0.resize(dmap0_cell0.size() + dmap0_cell1.size());
403 std::copy(dmap0_cell0.begin(), dmap0_cell0.end(), dmapjoint0.begin());
404 std::copy(dmap0_cell1.begin(), dmap0_cell1.end(),
405 std::next(dmapjoint0.begin(), dmap0_cell0.size()));
407 const tcb::span<const std::int32_t> dmap1_cell0 = dofmap1.links(cells[0]);
408 const tcb::span<const std::int32_t> dmap1_cell1 = dofmap1.links(cells[1]);
409 dmapjoint1.resize(dmap1_cell0.size() + dmap1_cell1.size());
410 std::copy(dmap1_cell0.begin(), dmap1_cell0.end(), dmapjoint1.begin());
411 std::copy(dmap1_cell1.begin(), dmap1_cell1.end(),
412 std::next(dmapjoint1.begin(), dmap1_cell0.size()));
416 for (std::size_t j = 0; j < dmap1_cell0.size(); ++j)
418 for (
int k = 0; k < bs1; ++k)
420 if (bc_markers1[bs1 * dmap1_cell0[j] + k])
429 for (std::size_t j = 0; j < dmap1_cell1.size(); ++j)
431 for (
int k = 0; k < bs1; ++k)
433 if (bc_markers1[bs1 * dmap1_cell1[j] + k])
446 const auto coeff_cell0 = coeffs.row(cells[0]);
447 const auto coeff_cell1 = coeffs.row(cells[1]);
450 for (std::size_t i = 0; i < offsets.size() - 1; ++i)
453 const int num_entries = offsets[i + 1] - offsets[i];
454 std::copy_n(coeff_cell0.data() + offsets[i], num_entries,
455 std::next(coeff_array.begin(), 2 * offsets[i]));
456 std::copy_n(coeff_cell1.data() + offsets[i], num_entries,
457 std::next(coeff_array.begin(), offsets[i + 1] + offsets[i]));
460 const int num_rows = bs0 * dmapjoint0.size();
461 const int num_cols = bs1 * dmapjoint1.size();
464 Ae.resize(num_rows * num_cols);
465 std::fill(Ae.begin(), Ae.end(), 0);
466 const int facets_per_cell = facets0.size();
467 const std::array perm{perms[
cells[0] * facets_per_cell + local_facet[0]],
468 perms[
cells[1] * facets_per_cell + local_facet[1]]};
469 kernel(Ae.data(), coeff_array.data(), constant_values.data(),
470 coordinate_dofs.data(), local_facet.data(), perm.data(),
471 cell_info[cells[0]]);
474 std::fill(be.begin(), be.end(), 0);
477 for (std::size_t j = 0; j < dmap1_cell0.size(); ++j)
479 for (
int k = 0; k < bs1; ++k)
481 const std::int32_t jj = bs1 * dmap1_cell0[j] + k;
484 const T bc = bc_values1[jj];
485 const T _x0 = x0.empty() ? 0.0 : x0[jj];
487 for (
int m = 0; m < num_rows; ++m)
488 be[m] -= Ae[m * num_cols + bs1 * j + k] * scale * (bc - _x0);
494 const int offset = bs1 * dmap1_cell0.size();
495 for (std::size_t j = 0; j < dmap1_cell1.size(); ++j)
497 for (
int k = 0; k < bs1; ++k)
499 const std::int32_t jj = bs1 * dmap1_cell1[j] + k;
502 const T bc = bc_values1[jj];
503 const T _x0 = x0.empty() ? 0.0 : x0[jj];
505 for (
int m = 0; m < num_rows; ++m)
508 -= Ae[m * num_cols + offset + bs1 * j + k] * scale * (bc - _x0);
514 for (std::size_t i = 0; i < dmap0_cell0.size(); ++i)
515 for (
int k = 0; k < bs0; ++k)
516 b[bs0 * dmap0_cell0[i] + k] += be[bs0 * i + k];
518 for (std::size_t i = 0; i < dmap0_cell1.size(); ++i)
519 for (
int k = 0; k < bs0; ++k)
520 b[bs0 * dmap0_cell1[i] + k] += be[offset + bs0 * i + k];
524 template <
typename T>
527 std::shared_ptr<const mesh::Mesh> mesh = L.mesh();
529 const int tdim = mesh->topology().dim();
530 const std::int32_t num_cells
531 = mesh->topology().connectivity(tdim, 0)->num_nodes();
534 assert(L.function_spaces().at(0));
535 std::shared_ptr<const fem::DofMap> dofmap
536 = L.function_spaces().at(0)->dofmap();
538 const graph::AdjacencyList<std::int32_t>& dofs = dofmap->list();
539 const int bs = dofmap->bs();
547 const bool needs_permutation_data = L.needs_permutation_data();
548 if (needs_permutation_data)
549 mesh->topology_mutable().create_entity_permutations();
550 const std::vector<std::uint32_t>& cell_info
551 = needs_permutation_data ? mesh->topology().get_cell_permutation_info()
552 : std::vector<std::uint32_t>(num_cells);
554 for (
int i : L.integral_ids(IntegralType::cell))
556 const auto& fn = L.kernel(IntegralType::cell, i);
557 const std::vector<std::int32_t>& active_cells
558 = L.domains(IntegralType::cell, i);
559 fem::impl::assemble_cells(b, mesh->geometry(), active_cells, dofs, bs, fn,
560 coeffs, constant_values, cell_info);
563 if (L.num_integrals(IntegralType::exterior_facet) > 0
564 or L.num_integrals(IntegralType::interior_facet) > 0)
567 mesh->topology_mutable().create_entities(tdim - 1);
568 mesh->topology_mutable().create_connectivity(tdim - 1, tdim);
569 mesh->topology_mutable().create_entity_permutations();
571 const std::vector<std::uint8_t>& perms
572 = mesh->topology().get_facet_permutations();
573 for (
int i : L.integral_ids(IntegralType::exterior_facet))
575 const auto& fn = L.kernel(IntegralType::exterior_facet, i);
576 const std::vector<std::int32_t>& active_facets
577 = L.domains(IntegralType::exterior_facet, i);
578 fem::impl::assemble_exterior_facets(b, *mesh, active_facets, dofs, bs, fn,
579 coeffs, constant_values, cell_info,
583 const std::vector<int> c_offsets = L.coefficient_offsets();
584 for (
int i : L.integral_ids(IntegralType::interior_facet))
586 const auto& fn = L.kernel(IntegralType::interior_facet, i);
587 const std::vector<std::int32_t>& active_facets
588 = L.domains(IntegralType::interior_facet, i);
589 fem::impl::assemble_interior_facets(b, *mesh, active_facets, *dofmap, fn,
590 coeffs, c_offsets, constant_values,
596 template <
typename T>
598 tcb::span<T> b,
const mesh::Geometry& geometry,
599 const std::vector<std::int32_t>& active_cells,
600 const graph::AdjacencyList<std::int32_t>& dofmap,
const int bs,
601 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
602 const std::uint8_t*,
const std::uint32_t)>& kernel,
603 const array2d<T>& coeffs,
const std::vector<T>& constant_values,
604 const std::vector<std::uint32_t>& cell_info)
606 const int gdim = geometry.dim();
609 const graph::AdjacencyList<std::int32_t>& x_dofmap = geometry.dofmap();
612 const int num_dofs_g = x_dofmap.num_links(0);
613 const array2d<double>& x_g = geometry.x();
617 const int num_dofs = dofmap.links(0).size();
618 std::vector<double> coordinate_dofs(num_dofs_g * gdim);
619 std::vector<T> be(bs * num_dofs);
622 for (std::int32_t c : active_cells)
625 auto x_dofs = x_dofmap.links(c);
626 for (std::size_t i = 0; i < x_dofs.size(); ++i)
628 std::copy_n(x_g.row(x_dofs[i]).begin(), gdim,
629 std::next(coordinate_dofs.begin(), i * gdim));
633 std::fill(be.begin(), be.end(), 0);
634 kernel(be.data(), coeffs.row(c).data(), constant_values.data(),
635 coordinate_dofs.data(),
nullptr,
nullptr, cell_info[c]);
638 auto dofs = dofmap.links(c);
639 for (
int i = 0; i < num_dofs; ++i)
640 for (
int k = 0; k < bs; ++k)
641 b[bs * dofs[i] + k] += be[bs * i + k];
645 template <
typename T>
646 void assemble_exterior_facets(
647 tcb::span<T> b,
const mesh::Mesh& mesh,
648 const std::vector<std::int32_t>& active_facets,
649 const graph::AdjacencyList<std::int32_t>& dofmap,
const int bs,
650 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
651 const std::uint8_t*,
const std::uint32_t)>& fn,
652 const array2d<T>& coeffs,
const std::vector<T>& constant_values,
653 const std::vector<std::uint32_t>& cell_info,
654 const std::vector<std::uint8_t>& perms)
656 const int gdim = mesh.geometry().dim();
657 const int tdim = mesh.topology().dim();
660 const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
663 const int num_dofs_g = x_dofmap.num_links(0);
664 const array2d<double>& x_g = mesh.geometry().x();
668 const int num_dofs = dofmap.links(0).size();
669 std::vector<double> coordinate_dofs(num_dofs_g * gdim);
670 std::vector<T> be(bs * num_dofs);
672 auto f_to_c = mesh.topology().connectivity(tdim - 1, tdim);
674 auto c_to_f = mesh.topology().connectivity(tdim, tdim - 1);
676 for (std::int32_t f : active_facets)
679 assert(f_to_c->num_links(f) > 0);
680 const std::int32_t cell = f_to_c->links(f)[0];
683 auto facets = c_to_f->links(cell);
684 auto it = std::find(facets.begin(), facets.end(), f);
685 assert(it != facets.end());
686 const int local_facet = std::distance(facets.begin(), it);
689 auto x_dofs = x_dofmap.links(cell);
690 for (std::size_t i = 0; i < x_dofs.size(); ++i)
692 std::copy_n(x_g.row(x_dofs[i]).begin(), gdim,
693 std::next(coordinate_dofs.begin(), i * gdim));
697 std::fill(be.begin(), be.end(), 0);
698 fn(be.data(), coeffs.row(cell).data(), constant_values.data(),
699 coordinate_dofs.data(), &local_facet,
700 &perms[cell * facets.size() + local_facet], cell_info[cell]);
703 auto dofs = dofmap.links(cell);
704 for (
int i = 0; i < num_dofs; ++i)
705 for (
int k = 0; k < bs; ++k)
706 b[bs * dofs[i] + k] += be[bs * i + k];
710 template <
typename T>
711 void assemble_interior_facets(
712 tcb::span<T> b,
const mesh::Mesh& mesh,
713 const std::vector<std::int32_t>& active_facets,
const fem::DofMap& dofmap,
714 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
715 const std::uint8_t*,
const std::uint32_t)>& fn,
716 const array2d<T>& coeffs,
const std::vector<int>& offsets,
717 const std::vector<T>& constant_values,
718 const std::vector<std::uint32_t>& cell_info,
719 const std::vector<std::uint8_t>& perms)
721 const int gdim = mesh.geometry().dim();
722 const int tdim = mesh.topology().dim();
725 const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
728 const int num_dofs_g = x_dofmap.num_links(0);
729 const array2d<double>& x_g = mesh.geometry().x();
732 std::vector<double> coordinate_dofs(2 * num_dofs_g * gdim);
734 std::vector<T> coeff_array(2 * offsets.back());
735 assert(offsets.back() ==
int(coeffs.shape[1]));
737 const int bs = dofmap.bs();
738 auto f_to_c = mesh.topology().connectivity(tdim - 1, tdim);
740 auto c_to_f = mesh.topology().connectivity(tdim, tdim - 1);
742 const int offset_g = gdim * num_dofs_g;
743 for (std::int32_t f : active_facets)
746 auto cells = f_to_c->links(f);
747 assert(
cells.size() == 2);
749 const int facets_per_cell = c_to_f->num_links(cells[0]);
752 std::array<int, 2> local_facet;
753 for (
int i = 0; i < 2; ++i)
755 auto facets = c_to_f->links(cells[i]);
756 auto it = std::find(facets.begin(), facets.end(), f);
757 assert(it != facets.end());
758 local_facet[i] = std::distance(facets.begin(), it);
762 auto x_dofs0 = x_dofmap.links(cells[0]);
763 auto x_dofs1 = x_dofmap.links(cells[1]);
764 for (
int i = 0; i < num_dofs_g; ++i)
766 for (
int j = 0; j < gdim; ++j)
768 coordinate_dofs[i * gdim + j] = x_g(x_dofs0[i], j);
769 coordinate_dofs[offset_g + i * gdim + j] = x_g(x_dofs1[i], j);
775 auto coeff_cell0 = coeffs.row(cells[0]);
776 auto coeff_cell1 = coeffs.row(cells[1]);
779 for (std::size_t i = 0; i < offsets.size() - 1; ++i)
782 const int num_entries = offsets[i + 1] - offsets[i];
783 std::copy_n(coeff_cell0.data() + offsets[i], num_entries,
784 std::next(coeff_array.begin(), 2 * offsets[i]));
785 std::copy_n(coeff_cell1.data() + offsets[i], num_entries,
786 std::next(coeff_array.begin(), offsets[i + 1] + offsets[i]));
790 tcb::span<const std::int32_t> dmap0 = dofmap.cell_dofs(cells[0]);
791 tcb::span<const std::int32_t> dmap1 = dofmap.cell_dofs(cells[1]);
794 be.resize(bs * (dmap0.size() + dmap1.size()));
795 std::fill(be.begin(), be.end(), 0);
796 const std::array perm{perms[
cells[0] * facets_per_cell + local_facet[0]],
797 perms[
cells[1] * facets_per_cell + local_facet[1]]};
798 fn(be.data(), coeff_array.data(), constant_values.data(),
799 coordinate_dofs.data(), local_facet.data(), perm.data(),
800 cell_info[cells[0]]);
803 for (std::size_t i = 0; i < dmap0.size(); ++i)
804 for (
int k = 0; k < bs; ++k)
805 b[bs * dmap0[i] + k] += be[bs * i + k];
806 for (std::size_t i = 0; i < dmap1.size(); ++i)
807 for (
int k = 0; k < bs; ++k)
808 b[bs * dmap1[i] + k] += be[bs * (i + dmap0.size()) + k];
812 template <
typename T>
814 tcb::span<T> b,
const std::vector<std::shared_ptr<
const Form<T>>> a,
815 const std::vector<std::vector<std::shared_ptr<
const DirichletBC<T>>>>& bcs1,
816 const std::vector<tcb::span<const T>>& x0,
double scale)
819 if (!x0.empty() and x0.size() != a.size())
821 throw std::runtime_error(
822 "Mismatch in size between x0 and bilinear form in assembler.");
825 if (a.size() != bcs1.size())
827 throw std::runtime_error(
828 "Mismatch in size between a and bcs in assembler.");
831 for (std::size_t j = 0; j < a.size(); ++j)
833 std::vector<bool> bc_markers1;
834 std::vector<T> bc_values1;
835 if (a[j] and !bcs1[j].empty())
837 auto V1 = a[j]->function_spaces()[1];
839 auto map1 = V1->dofmap()->index_map;
840 const int bs1 = V1->dofmap()->index_map_bs();
842 const int crange = bs1 * (map1->size_local() + map1->num_ghosts());
843 bc_markers1.assign(crange,
false);
844 bc_values1.assign(crange, 0.0);
845 for (
const std::shared_ptr<
const DirichletBC<T>>& bc : bcs1[j])
847 bc->mark_dofs(bc_markers1);
848 bc->dof_values(bc_values1);
852 lift_bc<T>(b, *a[j], bc_values1, bc_markers1, x0[j], scale);
855 lift_bc<T>(b, *a[j], bc_values1, bc_markers1, tcb::span<const T>(),
862 template <
typename T>
863 void lift_bc(tcb::span<T> b,
const Form<T>& a,
864 const tcb::span<const T>& bc_values1,
865 const std::vector<bool>& bc_markers1,
const tcb::span<const T>& x0,
868 std::shared_ptr<const mesh::Mesh> mesh = a.mesh();
870 const int tdim = mesh->topology().dim();
871 const std::int32_t num_cells
872 = mesh->topology().connectivity(tdim, 0)->num_nodes();
875 assert(a.function_spaces().at(0));
876 assert(a.function_spaces().at(1));
877 const graph::AdjacencyList<std::int32_t>& dofmap0
878 = a.function_spaces()[0]->dofmap()->list();
879 const int bs0 = a.function_spaces()[0]->dofmap()->bs();
880 const graph::AdjacencyList<std::int32_t>& dofmap1
881 = a.function_spaces()[1]->dofmap()->list();
882 const int bs1 = a.function_spaces()[1]->dofmap()->bs();
890 const bool needs_permutation_data = a.needs_permutation_data();
891 if (needs_permutation_data)
892 mesh->topology_mutable().create_entity_permutations();
893 const std::vector<std::uint32_t>& cell_info
894 = needs_permutation_data ? mesh->topology().get_cell_permutation_info()
895 : std::vector<std::uint32_t>(num_cells);
897 for (
int i : a.integral_ids(IntegralType::cell))
899 const auto& kernel = a.kernel(IntegralType::cell, i);
900 const std::vector<std::int32_t>& active_cells
901 = a.domains(IntegralType::cell, i);
902 _lift_bc_cells(b, mesh->geometry(), kernel, active_cells, dofmap0, bs0,
903 dofmap1, bs1, coeffs, constant_values, cell_info, bc_values1,
904 bc_markers1, x0, scale);
907 if (a.num_integrals(IntegralType::exterior_facet) > 0
908 or a.num_integrals(IntegralType::interior_facet) > 0)
911 mesh->topology_mutable().create_entities(tdim - 1);
912 mesh->topology_mutable().create_connectivity(tdim - 1, tdim);
913 mesh->topology_mutable().create_entity_permutations();
915 const std::vector<std::uint8_t>& perms
916 = mesh->topology().get_facet_permutations();
917 for (
int i : a.integral_ids(IntegralType::exterior_facet))
919 const auto& kernel = a.kernel(IntegralType::exterior_facet, i);
920 const std::vector<std::int32_t>& active_facets
921 = a.domains(IntegralType::exterior_facet, i);
922 _lift_bc_exterior_facets(b, *mesh, kernel, active_facets, dofmap0, bs0,
923 dofmap1, bs1, coeffs, constant_values, cell_info,
924 perms, bc_values1, bc_markers1, x0, scale);
927 const std::vector<int> c_offsets = a.coefficient_offsets();
928 for (
int i : a.integral_ids(IntegralType::interior_facet))
930 const auto& kernel = a.kernel(IntegralType::interior_facet, i);
931 const std::vector<std::int32_t>& active_facets
932 = a.domains(IntegralType::interior_facet, i);
933 _lift_bc_interior_facets(b, *mesh, kernel, active_facets, dofmap0, bs0,
934 dofmap1, bs1, coeffs, c_offsets, constant_values,
935 cell_info, perms, bc_values1, bc_markers1, x0,
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
std::vector< typename U::scalar_type > pack_constants(const U &u)
Pack constants of u of generic type U ready for assembly.
Definition: utils.h:453
void apply_lifting(tcb::span< T > b, const std::vector< std::shared_ptr< const Form< T >>> &a, const std::vector< std::vector< std::shared_ptr< const DirichletBC< T >>>> &bcs1, const std::vector< tcb::span< const T >> &x0, double scale)
Modify b such that:
Definition: assembler.h:69
array2d< typename U::scalar_type > pack_coefficients(const U &u)
Pack coefficients of u of generic type U ready for assembly.
Definition: utils.h:398
void assemble_vector(tcb::span< T > b, const Form< T > &L)
Assemble linear form into a vector.
Definition: assembler.h:45