9 #include "DirichletBC.h"
13 #include <Eigen/Dense>
14 #include <dolfinx/common/IndexMap.h>
15 #include <dolfinx/common/types.h>
16 #include <dolfinx/fem/Constant.h>
17 #include <dolfinx/fem/FunctionSpace.h>
18 #include <dolfinx/graph/AdjacencyList.h>
19 #include <dolfinx/mesh/Geometry.h>
20 #include <dolfinx/mesh/Mesh.h>
21 #include <dolfinx/mesh/Topology.h>
26 namespace dolfinx::fem::impl
41 tcb::span<T> b,
const mesh::Geometry& geometry,
42 const std::vector<std::int32_t>& active_cells,
43 const graph::AdjacencyList<std::int32_t>& dofmap,
const int bs,
44 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
45 const std::uint8_t*,
const std::uint32_t)>& kernel,
46 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
48 const std::vector<T>& constant_values,
49 const std::vector<std::uint32_t>& cell_info);
53 void assemble_exterior_facets(
54 tcb::span<T> b,
const mesh::Mesh& mesh,
55 const std::vector<std::int32_t>& active_facets,
56 const graph::AdjacencyList<std::int32_t>& dofmap,
const int bs,
57 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
58 const std::uint8_t*,
const std::uint32_t)>& fn,
59 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
61 const std::vector<T>& constant_values,
62 const std::vector<std::uint32_t>& cell_info,
63 const std::vector<std::uint8_t>& perms);
67 void assemble_interior_facets(
68 tcb::span<T> b,
const mesh::Mesh& mesh,
69 const std::vector<std::int32_t>& active_facets,
const fem::DofMap& dofmap,
70 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
71 const std::uint8_t*,
const std::uint32_t)>& fn,
72 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
74 const std::vector<int>& offsets,
const std::vector<T>& constant_values,
75 const std::vector<std::uint32_t>& cell_info,
76 const std::vector<std::uint8_t>& perms);
97 tcb::span<T> b,
const std::vector<std::shared_ptr<
const Form<T>>> a,
98 const std::vector<std::vector<std::shared_ptr<
const DirichletBC<T>>>>& bcs1,
99 const tcb::span<const T>& x0,
double scale);
113 template <
typename T>
114 void lift_bc(tcb::span<T> b,
const Form<T>& a,
115 const tcb::span<const T>& bc_values1,
116 const std::vector<bool>& bc_markers1,
const tcb::span<const T>& x0,
120 template <
typename T>
122 tcb::span<T> b,
const mesh::Geometry& geometry,
123 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
124 const std::uint8_t*,
const std::uint32_t)>& kernel,
125 const std::vector<std::int32_t>& active_cells,
126 const graph::AdjacencyList<std::int32_t>& dofmap0,
int bs0,
127 const graph::AdjacencyList<std::int32_t>& dofmap1,
int bs1,
128 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
130 const std::vector<T>& constant_values,
131 const std::vector<std::uint32_t>& cell_info,
132 const tcb::span<const T>& bc_values1,
const std::vector<bool>& bc_markers1,
133 const tcb::span<const T>& x0,
double scale)
136 const int gdim = geometry.dim();
137 const graph::AdjacencyList<std::int32_t>& x_dofmap = geometry.dofmap();
140 const int num_dofs_g = x_dofmap.num_links(0);
141 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
145 std::vector<double> coordinate_dofs(num_dofs_g * gdim);
146 std::vector<T> Ae, be;
148 for (std::int32_t c : active_cells)
151 auto dmap1 = dofmap1.links(c);
155 for (std::size_t j = 0; j < dmap1.size(); ++j)
157 for (
int k = 0; k < bs1; ++k)
159 assert(bs1 * dmap1[j] + k < (
int)bc_markers1.size());
160 if (bc_markers1[bs1 * dmap1[j] + k])
172 auto x_dofs = x_dofmap.links(c);
173 for (std::size_t i = 0; i < x_dofs.size(); ++i)
175 std::copy_n(x_g.row(x_dofs[i]).data(), gdim,
176 std::next(coordinate_dofs.begin(), i * gdim));
180 auto dmap0 = dofmap0.links(c);
182 const int num_rows = bs0 * dmap0.size();
183 const int num_cols = bs1 * dmap1.size();
185 auto coeff_array = coeffs.row(c);
186 Ae.resize(num_rows * num_cols);
187 std::fill(Ae.begin(), Ae.end(), 0);
188 kernel(Ae.data(), coeff_array.data(), constant_values.data(),
189 coordinate_dofs.data(),
nullptr,
nullptr, cell_info[c]);
193 std::fill(be.begin(), be.end(), 0);
194 for (std::size_t j = 0; j < dmap1.size(); ++j)
196 for (
int k = 0; k < bs1; ++k)
198 const std::int32_t jj = bs1 * dmap1[j] + k;
199 assert(jj < (
int)bc_markers1.size());
202 const T bc = bc_values1[jj];
203 const T _x0 = x0.empty() ? 0.0 : x0[jj];
205 for (
int m = 0; m < num_rows; ++m)
206 be[m] -= Ae[m * num_cols + bs1 * j + k] * scale * (bc - _x0);
211 for (std::size_t i = 0; i < dmap0.size(); ++i)
212 for (
int k = 0; k < bs0; ++k)
213 b[bs0 * dmap0[i] + k] += be[bs0 * i + k];
217 template <
typename T>
218 void _lift_bc_exterior_facets(
219 tcb::span<T> b,
const mesh::Mesh& mesh,
220 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
221 const std::uint8_t*,
const std::uint32_t)>& kernel,
222 const std::vector<std::int32_t>& active_facets,
223 const graph::AdjacencyList<std::int32_t>& dofmap0,
int bs0,
224 const graph::AdjacencyList<std::int32_t>& dofmap1,
int bs1,
225 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
227 const std::vector<T>& constant_values,
228 const std::vector<std::uint32_t>& cell_info,
229 const std::vector<std::uint8_t>& perms,
230 const tcb::span<const T>& bc_values1,
const std::vector<bool>& bc_markers1,
231 const tcb::span<const T>& x0,
double scale)
233 const int gdim = mesh.geometry().dim();
234 const int tdim = mesh.topology().dim();
237 const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
240 const int num_dofs_g = x_dofmap.num_links(0);
241 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
242 = mesh.geometry().x();
245 std::vector<double> coordinate_dofs(num_dofs_g * gdim);
246 std::vector<T> Ae, be;
249 const mesh::Topology& topology = mesh.topology();
250 auto connectivity = topology.connectivity(tdim - 1, tdim);
251 assert(connectivity);
252 auto c_to_f = topology.connectivity(tdim, tdim - 1);
254 auto map = topology.index_map(tdim - 1);
257 for (std::int32_t f : active_facets)
260 assert(connectivity->num_links(f) == 1);
261 const std::int32_t cell = connectivity->links(f)[0];
264 auto facets = c_to_f->links(cell);
265 auto it = std::find(facets.begin(), facets.end(), f);
266 assert(it != facets.end());
267 const int local_facet = std::distance(facets.begin(), it);
270 auto dmap1 = dofmap1.links(cell);
274 for (std::size_t j = 0; j < dmap1.size(); ++j)
276 for (
int k = 0; k < bs1; ++k)
278 if (bc_markers1[bs1 * dmap1[j] + k])
290 auto x_dofs = x_dofmap.links(cell);
291 for (std::size_t i = 0; i < x_dofs.size(); ++i)
293 std::copy_n(x_g.row(x_dofs[i]).data(), gdim,
294 std::next(coordinate_dofs.begin(), i * gdim));
298 auto dmap0 = dofmap0.links(cell);
300 const int num_rows = bs0 * dmap0.size();
301 const int num_cols = bs1 * dmap1.size();
303 auto coeff_array = coeffs.row(cell);
304 Ae.resize(num_rows * num_cols);
305 std::fill(Ae.begin(), Ae.end(), 0);
306 kernel(Ae.data(), coeff_array.data(), constant_values.data(),
307 coordinate_dofs.data(), &local_facet,
308 &perms[cell * facets.size() + local_facet], cell_info[cell]);
312 std::fill(be.begin(), be.end(), 0);
313 for (std::size_t j = 0; j < dmap1.size(); ++j)
315 for (
int k = 0; k < bs1; ++k)
317 const std::int32_t jj = bs1 * dmap1[j] + k;
321 const T bc = bc_values1[jj];
322 const T _x0 = x0.empty() ? 0.0 : x0[jj];
324 for (
int m = 0; m < num_rows; ++m)
325 be[m] -= Ae[m * num_cols + bs1 * j + k] * scale * (bc - _x0);
330 for (std::size_t i = 0; i < dmap0.size(); ++i)
331 for (
int k = 0; k < bs0; ++k)
332 b[bs0 * dmap0[i] + k] += be[bs0 * i + k];
336 template <
typename T>
337 void _lift_bc_interior_facets(
338 tcb::span<T> b,
const mesh::Mesh& mesh,
339 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
340 const std::uint8_t*,
const std::uint32_t)>& kernel,
341 const std::vector<std::int32_t>& active_facets,
342 const graph::AdjacencyList<std::int32_t>& dofmap0,
int bs0,
343 const graph::AdjacencyList<std::int32_t>& dofmap1,
int bs1,
344 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
346 const std::vector<int>& offsets,
const std::vector<T>& constant_values,
347 const std::vector<std::uint32_t>& cell_info,
348 const std::vector<std::uint8_t>& perms,
349 const tcb::span<const T>& bc_values1,
const std::vector<bool>& bc_markers1,
350 const tcb::span<const T>& x0,
double scale)
352 const int gdim = mesh.geometry().dim();
353 const int tdim = mesh.topology().dim();
356 const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
359 const int num_dofs_g = x_dofmap.num_links(0);
360 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
361 = mesh.geometry().x();
364 std::vector<double> coordinate_dofs(2 * num_dofs_g * gdim);
365 std::vector<T> coeff_array(2 * offsets.back());
366 assert(offsets.back() == coeffs.cols());
367 std::vector<T> Ae, be;
370 std::vector<std::int32_t> dmapjoint0, dmapjoint1;
372 const mesh::Topology& topology = mesh.topology();
373 auto connectivity = topology.connectivity(tdim - 1, tdim);
374 assert(connectivity);
375 auto c_to_f = topology.connectivity(tdim, tdim - 1);
377 auto f_to_c = topology.connectivity(tdim - 1, tdim);
379 auto map = topology.index_map(tdim - 1);
382 const int offset_g = gdim * num_dofs_g;
383 for (std::int32_t f : active_facets)
386 auto cells = f_to_c->links(f);
387 assert(
cells.size() == 2);
390 auto facets0 = c_to_f->links(cells[0]);
391 auto it0 = std::find(facets0.begin(), facets0.end(), f);
392 assert(it0 != facets0.end());
393 const int local_facet0 = std::distance(facets0.begin(), it0);
394 auto facets1 = c_to_f->links(cells[1]);
395 auto it1 = std::find(facets1.begin(), facets1.end(), f);
396 assert(it1 != facets1.end());
397 const int local_facet1 = std::distance(facets1.begin(), it1);
399 const std::array local_facet{local_facet0, local_facet1};
402 auto x_dofs0 = x_dofmap.links(cells[0]);
403 auto x_dofs1 = x_dofmap.links(cells[1]);
404 for (
int i = 0; i < num_dofs_g; ++i)
406 for (
int j = 0; j < gdim; ++j)
408 coordinate_dofs[i * gdim + j] = x_g(x_dofs0[i], j);
409 coordinate_dofs[offset_g + i * gdim + j] = x_g(x_dofs1[i], j);
414 const tcb::span<const std::int32_t> dmap0_cell0 = dofmap0.links(cells[0]);
415 const tcb::span<const std::int32_t> dmap0_cell1 = dofmap0.links(cells[1]);
416 dmapjoint0.resize(dmap0_cell0.size() + dmap0_cell1.size());
417 std::copy(dmap0_cell0.begin(), dmap0_cell0.end(), dmapjoint0.begin());
418 std::copy(dmap0_cell1.begin(), dmap0_cell1.end(),
419 std::next(dmapjoint0.begin(), dmap0_cell0.size()));
421 const tcb::span<const std::int32_t> dmap1_cell0 = dofmap1.links(cells[0]);
422 const tcb::span<const std::int32_t> dmap1_cell1 = dofmap1.links(cells[1]);
423 dmapjoint1.resize(dmap1_cell0.size() + dmap1_cell1.size());
424 std::copy(dmap1_cell0.begin(), dmap1_cell0.end(), dmapjoint1.begin());
425 std::copy(dmap1_cell1.begin(), dmap1_cell1.end(),
426 std::next(dmapjoint1.begin(), dmap1_cell0.size()));
430 for (std::size_t j = 0; j < dmap1_cell0.size(); ++j)
432 for (
int k = 0; k < bs1; ++k)
434 if (bc_markers1[bs1 * dmap1_cell0[j] + k])
443 for (std::size_t j = 0; j < dmap1_cell1.size(); ++j)
445 for (
int k = 0; k < bs1; ++k)
447 if (bc_markers1[bs1 * dmap1_cell1[j] + k])
460 auto coeff_cell0 = coeffs.row(cells[0]);
461 auto coeff_cell1 = coeffs.row(cells[1]);
464 for (std::size_t i = 0; i < offsets.size() - 1; ++i)
467 const int num_entries = offsets[i + 1] - offsets[i];
468 std::copy_n(coeff_cell0.data() + offsets[i], num_entries,
469 std::next(coeff_array.begin(), 2 * offsets[i]));
470 std::copy_n(coeff_cell1.data() + offsets[i], num_entries,
471 std::next(coeff_array.begin(), offsets[i + 1] + offsets[i]));
474 const int num_rows = bs0 * dmapjoint0.size();
475 const int num_cols = bs1 * dmapjoint1.size();
478 Ae.resize(num_rows * num_cols);
479 std::fill(Ae.begin(), Ae.end(), 0);
480 const int facets_per_cell = facets0.size();
481 const std::array perm{perms[
cells[0] * facets_per_cell + local_facet[0]],
482 perms[
cells[1] * facets_per_cell + local_facet[1]]};
483 kernel(Ae.data(), coeff_array.data(), constant_values.data(),
484 coordinate_dofs.data(), local_facet.data(), perm.data(),
485 cell_info[cells[0]]);
488 std::fill(be.begin(), be.end(), 0);
491 for (std::size_t j = 0; j < dmap1_cell0.size(); ++j)
493 for (
int k = 0; k < bs1; ++k)
495 const std::int32_t jj = bs1 * dmap1_cell0[j] + k;
498 const T bc = bc_values1[jj];
499 const T _x0 = x0.empty() ? 0.0 : x0[jj];
501 for (
int m = 0; m < num_rows; ++m)
502 be[m] -= Ae[m * num_cols + bs1 * j + k] * scale * (bc - _x0);
508 const int offset = bs1 * dmap1_cell0.size();
509 for (std::size_t j = 0; j < dmap1_cell1.size(); ++j)
511 for (
int k = 0; k < bs1; ++k)
513 const std::int32_t jj = bs1 * dmap1_cell1[j] + k;
516 const T bc = bc_values1[jj];
517 const T _x0 = x0.empty() ? 0.0 : x0[jj];
519 for (
int m = 0; m < num_rows; ++m)
522 -= Ae[m * num_cols + offset + bs1 * j + k] * scale * (bc - _x0);
528 for (std::size_t i = 0; i < dmap0_cell0.size(); ++i)
529 for (
int k = 0; k < bs0; ++k)
530 b[bs0 * dmap0_cell0[i] + k] += be[bs0 * i + k];
532 for (std::size_t i = 0; i < dmap0_cell1.size(); ++i)
533 for (
int k = 0; k < bs0; ++k)
534 b[bs0 * dmap0_cell1[i] + k] += be[offset + bs0 * i + k];
538 template <
typename T>
541 std::shared_ptr<const mesh::Mesh> mesh = L.mesh();
543 const int tdim = mesh->topology().dim();
544 const std::int32_t num_cells
545 = mesh->topology().connectivity(tdim, 0)->num_nodes();
548 assert(L.function_spaces().at(0));
549 std::shared_ptr<const fem::DofMap> dofmap
550 = L.function_spaces().at(0)->dofmap();
552 const graph::AdjacencyList<std::int32_t>& dofs = dofmap->list();
553 const int bs = dofmap->bs();
559 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> coeffs
562 const bool needs_permutation_data = L.needs_permutation_data();
563 if (needs_permutation_data)
564 mesh->topology_mutable().create_entity_permutations();
565 const std::vector<std::uint32_t>& cell_info
566 = needs_permutation_data ? mesh->topology().get_cell_permutation_info()
567 : std::vector<std::uint32_t>(num_cells);
569 for (
int i : L.integral_ids(IntegralType::cell))
571 const auto& fn = L.kernel(IntegralType::cell, i);
572 const std::vector<std::int32_t>& active_cells
573 = L.domains(IntegralType::cell, i);
574 fem::impl::assemble_cells(b, mesh->geometry(), active_cells, dofs, bs, fn,
575 coeffs, constant_values, cell_info);
578 if (L.num_integrals(IntegralType::exterior_facet) > 0
579 or L.num_integrals(IntegralType::interior_facet) > 0)
582 mesh->topology_mutable().create_entities(tdim - 1);
583 mesh->topology_mutable().create_connectivity(tdim - 1, tdim);
584 mesh->topology_mutable().create_entity_permutations();
586 const std::vector<std::uint8_t>& perms
587 = mesh->topology().get_facet_permutations();
588 for (
int i : L.integral_ids(IntegralType::exterior_facet))
590 const auto& fn = L.kernel(IntegralType::exterior_facet, i);
591 const std::vector<std::int32_t>& active_facets
592 = L.domains(IntegralType::exterior_facet, i);
593 fem::impl::assemble_exterior_facets(b, *mesh, active_facets, dofs, bs, fn,
594 coeffs, constant_values, cell_info,
598 const std::vector<int> c_offsets = L.coefficient_offsets();
599 for (
int i : L.integral_ids(IntegralType::interior_facet))
601 const auto& fn = L.kernel(IntegralType::interior_facet, i);
602 const std::vector<std::int32_t>& active_facets
603 = L.domains(IntegralType::interior_facet, i);
604 fem::impl::assemble_interior_facets(b, *mesh, active_facets, *dofmap, fn,
605 coeffs, c_offsets, constant_values,
611 template <
typename T>
613 tcb::span<T> b,
const mesh::Geometry& geometry,
614 const std::vector<std::int32_t>& active_cells,
615 const graph::AdjacencyList<std::int32_t>& dofmap,
const int bs,
616 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
617 const std::uint8_t*,
const std::uint32_t)>& kernel,
618 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
620 const std::vector<T>& constant_values,
621 const std::vector<std::uint32_t>& cell_info)
623 const int gdim = geometry.dim();
626 const graph::AdjacencyList<std::int32_t>& x_dofmap = geometry.dofmap();
629 const int num_dofs_g = x_dofmap.num_links(0);
630 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
635 const int num_dofs = dofmap.links(0).size();
636 std::vector<double> coordinate_dofs(num_dofs_g * gdim);
637 std::vector<T> be(bs * num_dofs);
640 for (std::int32_t c : active_cells)
643 auto x_dofs = x_dofmap.links(c);
644 for (std::size_t i = 0; i < x_dofs.size(); ++i)
646 std::copy_n(x_g.row(x_dofs[i]).data(), gdim,
647 std::next(coordinate_dofs.begin(), i * gdim));
651 std::fill(be.begin(), be.end(), 0);
652 kernel(be.data(), coeffs.row(c).data(), constant_values.data(),
653 coordinate_dofs.data(),
nullptr,
nullptr, cell_info[c]);
656 auto dofs = dofmap.links(c);
657 for (
int i = 0; i < num_dofs; ++i)
658 for (
int k = 0; k < bs; ++k)
659 b[bs * dofs[i] + k] += be[bs * i + k];
663 template <
typename T>
664 void assemble_exterior_facets(
665 tcb::span<T> b,
const mesh::Mesh& mesh,
666 const std::vector<std::int32_t>& active_facets,
667 const graph::AdjacencyList<std::int32_t>& dofmap,
const int bs,
668 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
669 const std::uint8_t*,
const std::uint32_t)>& fn,
670 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
672 const std::vector<T>& constant_values,
673 const std::vector<std::uint32_t>& cell_info,
674 const std::vector<std::uint8_t>& perms)
676 const int gdim = mesh.geometry().dim();
677 const int tdim = mesh.topology().dim();
680 const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
683 const int num_dofs_g = x_dofmap.num_links(0);
684 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
685 = mesh.geometry().x();
689 const int num_dofs = dofmap.links(0).size();
690 std::vector<double> coordinate_dofs(num_dofs_g * gdim);
691 std::vector<T> be(bs * num_dofs);
693 auto f_to_c = mesh.topology().connectivity(tdim - 1, tdim);
695 auto c_to_f = mesh.topology().connectivity(tdim, tdim - 1);
697 for (
const auto& f : active_facets)
700 assert(f_to_c->num_links(f) > 0);
701 const std::int32_t cell = f_to_c->links(f)[0];
704 auto facets = c_to_f->links(cell);
705 auto it = std::find(facets.begin(), facets.end(), f);
706 assert(it != facets.end());
707 const int local_facet = std::distance(facets.begin(), it);
710 auto x_dofs = x_dofmap.links(cell);
711 for (std::size_t i = 0; i < x_dofs.size(); ++i)
713 std::copy_n(x_g.row(x_dofs[i]).data(), gdim,
714 std::next(coordinate_dofs.begin(), i * gdim));
718 std::fill(be.begin(), be.end(), 0);
719 fn(be.data(), coeffs.row(cell).data(), constant_values.data(),
720 coordinate_dofs.data(), &local_facet,
721 &perms[cell * facets.size() + local_facet], cell_info[cell]);
724 auto dofs = dofmap.links(cell);
725 for (
int i = 0; i < num_dofs; ++i)
726 for (
int k = 0; k < bs; ++k)
727 b[bs * dofs[i] + k] += be[bs * i + k];
731 template <
typename T>
732 void assemble_interior_facets(
733 tcb::span<T> b,
const mesh::Mesh& mesh,
734 const std::vector<std::int32_t>& active_facets,
const fem::DofMap& dofmap,
735 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
736 const std::uint8_t*,
const std::uint32_t)>& fn,
737 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
739 const std::vector<int>& offsets,
const std::vector<T>& constant_values,
740 const std::vector<std::uint32_t>& cell_info,
741 const std::vector<std::uint8_t>& perms)
743 const int gdim = mesh.geometry().dim();
744 const int tdim = mesh.topology().dim();
747 const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
750 const int num_dofs_g = x_dofmap.num_links(0);
751 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
752 = mesh.geometry().x();
755 std::vector<double> coordinate_dofs(2 * num_dofs_g * gdim);
757 std::vector<T> coeff_array(2 * offsets.back());
758 assert(offsets.back() == coeffs.cols());
760 const int bs = dofmap.bs();
761 auto f_to_c = mesh.topology().connectivity(tdim - 1, tdim);
763 auto c_to_f = mesh.topology().connectivity(tdim, tdim - 1);
765 const int offset_g = gdim * num_dofs_g;
766 for (
const auto& f : active_facets)
769 auto cells = f_to_c->links(f);
770 assert(
cells.size() == 2);
772 const int facets_per_cell = c_to_f->num_links(cells[0]);
775 std::array<int, 2> local_facet;
776 for (
int i = 0; i < 2; ++i)
778 auto facets = c_to_f->links(cells[i]);
779 auto it = std::find(facets.begin(), facets.end(), f);
780 assert(it != facets.end());
781 local_facet[i] = std::distance(facets.begin(), it);
785 auto x_dofs0 = x_dofmap.links(cells[0]);
786 auto x_dofs1 = x_dofmap.links(cells[1]);
787 for (
int i = 0; i < num_dofs_g; ++i)
789 for (
int j = 0; j < gdim; ++j)
791 coordinate_dofs[i * gdim + j] = x_g(x_dofs0[i], j);
792 coordinate_dofs[offset_g + i * gdim + j] = x_g(x_dofs1[i], j);
798 auto coeff_cell0 = coeffs.row(cells[0]);
799 auto coeff_cell1 = coeffs.row(cells[1]);
802 for (std::size_t i = 0; i < offsets.size() - 1; ++i)
805 const int num_entries = offsets[i + 1] - offsets[i];
806 std::copy_n(coeff_cell0.data() + offsets[i], num_entries,
807 std::next(coeff_array.begin(), 2 * offsets[i]));
808 std::copy_n(coeff_cell1.data() + offsets[i], num_entries,
809 std::next(coeff_array.begin(), offsets[i + 1] + offsets[i]));
813 tcb::span<const std::int32_t> dmap0 = dofmap.cell_dofs(cells[0]);
814 tcb::span<const std::int32_t> dmap1 = dofmap.cell_dofs(cells[1]);
817 be.resize(bs * (dmap0.size() + dmap1.size()));
818 std::fill(be.begin(), be.end(), 0);
819 const std::array perm{perms[
cells[0] * facets_per_cell + local_facet[0]],
820 perms[
cells[1] * facets_per_cell + local_facet[1]]};
821 fn(be.data(), coeff_array.data(), constant_values.data(),
822 coordinate_dofs.data(), local_facet.data(), perm.data(),
823 cell_info[cells[0]]);
826 for (std::size_t i = 0; i < dmap0.size(); ++i)
827 for (
int k = 0; k < bs; ++k)
828 b[bs * dmap0[i] + k] += be[bs * i + k];
829 for (std::size_t i = 0; i < dmap1.size(); ++i)
830 for (
int k = 0; k < bs; ++k)
831 b[bs * dmap1[i] + k] += be[bs * (i + dmap0.size()) + k];
835 template <
typename T>
837 tcb::span<T> b,
const std::vector<std::shared_ptr<
const Form<T>>> a,
838 const std::vector<std::vector<std::shared_ptr<
const DirichletBC<T>>>>& bcs1,
839 const std::vector<tcb::span<const T>>& x0,
double scale)
842 if (!x0.empty() and x0.size() != a.size())
844 throw std::runtime_error(
845 "Mismatch in size between x0 and bilinear form in assembler.");
848 if (a.size() != bcs1.size())
850 throw std::runtime_error(
851 "Mismatch in size between a and bcs in assembler.");
854 for (std::size_t j = 0; j < a.size(); ++j)
856 std::vector<bool> bc_markers1;
857 std::vector<T> bc_values1;
858 if (a[j] and !bcs1[j].empty())
860 auto V1 = a[j]->function_spaces()[1];
862 auto map1 = V1->dofmap()->index_map;
863 const int bs1 = V1->dofmap()->index_map_bs();
865 const int crange = bs1 * (map1->size_local() + map1->num_ghosts());
866 bc_markers1.assign(crange,
false);
867 bc_values1.assign(crange, 0.0);
868 for (
const std::shared_ptr<
const DirichletBC<T>>& bc : bcs1[j])
870 bc->mark_dofs(bc_markers1);
871 bc->dof_values(bc_values1);
875 lift_bc<T>(b, *a[j], bc_values1, bc_markers1, x0[j], scale);
878 lift_bc<T>(b, *a[j], bc_values1, bc_markers1, tcb::span<const T>(),
885 template <
typename T>
886 void lift_bc(tcb::span<T> b,
const Form<T>& a,
887 const tcb::span<const T>& bc_values1,
888 const std::vector<bool>& bc_markers1,
const tcb::span<const T>& x0,
891 std::shared_ptr<const mesh::Mesh> mesh = a.mesh();
893 const int tdim = mesh->topology().dim();
894 const std::int32_t num_cells
895 = mesh->topology().connectivity(tdim, 0)->num_nodes();
898 assert(a.function_spaces().at(0));
899 assert(a.function_spaces().at(1));
900 const graph::AdjacencyList<std::int32_t>& dofmap0
901 = a.function_spaces()[0]->dofmap()->list();
902 const int bs0 = a.function_spaces()[0]->dofmap()->bs();
903 const graph::AdjacencyList<std::int32_t>& dofmap1
904 = a.function_spaces()[1]->dofmap()->list();
905 const int bs1 = a.function_spaces()[1]->dofmap()->bs();
911 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> coeffs
914 const bool needs_permutation_data = a.needs_permutation_data();
915 if (needs_permutation_data)
916 mesh->topology_mutable().create_entity_permutations();
917 const std::vector<std::uint32_t>& cell_info
918 = needs_permutation_data ? mesh->topology().get_cell_permutation_info()
919 : std::vector<std::uint32_t>(num_cells);
921 for (
int i : a.integral_ids(IntegralType::cell))
923 const auto& kernel = a.kernel(IntegralType::cell, i);
924 const std::vector<std::int32_t>& active_cells
925 = a.domains(IntegralType::cell, i);
926 _lift_bc_cells(b, mesh->geometry(), kernel, active_cells, dofmap0, bs0,
927 dofmap1, bs1, coeffs, constant_values, cell_info, bc_values1,
928 bc_markers1, x0, scale);
931 if (a.num_integrals(IntegralType::exterior_facet) > 0
932 or a.num_integrals(IntegralType::interior_facet) > 0)
935 mesh->topology_mutable().create_entities(tdim - 1);
936 mesh->topology_mutable().create_connectivity(tdim - 1, tdim);
937 mesh->topology_mutable().create_entity_permutations();
939 const std::vector<std::uint8_t>& perms
940 = mesh->topology().get_facet_permutations();
941 for (
int i : a.integral_ids(IntegralType::exterior_facet))
943 const auto& kernel = a.kernel(IntegralType::exterior_facet, i);
944 const std::vector<std::int32_t>& active_facets
945 = a.domains(IntegralType::exterior_facet, i);
946 _lift_bc_exterior_facets(b, *mesh, kernel, active_facets, dofmap0, bs0,
947 dofmap1, bs1, coeffs, constant_values, cell_info,
948 perms, bc_values1, bc_markers1, x0, scale);
951 const std::vector<int> c_offsets = a.coefficient_offsets();
952 for (
int i : a.integral_ids(IntegralType::interior_facet))
954 const auto& kernel = a.kernel(IntegralType::interior_facet, i);
955 const std::vector<std::int32_t>& active_facets
956 = a.domains(IntegralType::interior_facet, i);
957 _lift_bc_interior_facets(b, *mesh, kernel, active_facets, dofmap0, bs0,
958 dofmap1, bs1, coeffs, c_offsets, constant_values,
959 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:455
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
Eigen::Array< typename U::scalar_type, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > pack_coefficients(const U &u)
Pack coefficients of u of generic type U ready for assembly.
Definition: utils.h:400
void assemble_vector(tcb::span< T > b, const Form< T > &L)
Assemble linear form into a vector.
Definition: assembler.h:45