9 #include "DirichletBC.h"
13 #include <Eigen/Dense>
14 #include <dolfinx/common/IndexMap.h>
15 #include <dolfinx/common/types.h>
16 #include <dolfinx/function/Constant.h>
17 #include <dolfinx/function/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
42 Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
43 const mesh::Geometry& geometry,
44 const std::vector<std::int32_t>& active_cells,
45 const graph::AdjacencyList<std::int32_t>& dofmap,
46 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
47 const std::uint8_t*,
const std::uint32_t)>& kernel,
48 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
50 const Eigen::Array<T, Eigen::Dynamic, 1>& constant_values,
51 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info);
55 void assemble_exterior_facets(
56 Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
const mesh::Mesh& mesh,
57 const std::vector<std::int32_t>& active_facets,
58 const graph::AdjacencyList<std::int32_t>& dofmap,
59 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
60 const std::uint8_t*,
const std::uint32_t)>& fn,
61 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
63 const Eigen::Array<T, Eigen::Dynamic, 1>& constant_values,
64 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info,
65 const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>& perms);
69 void assemble_interior_facets(
70 Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
const mesh::Mesh& mesh,
71 const std::vector<std::int32_t>& active_facets,
const fem::DofMap& dofmap,
72 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
73 const std::uint8_t*,
const std::uint32_t)>& fn,
74 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
76 const std::vector<int>& offsets,
77 const Eigen::Array<T, Eigen::Dynamic, 1>& constant_values,
78 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info,
79 const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>& perms);
100 Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
101 const std::vector<std::shared_ptr<
const Form<T>>> a,
102 const std::vector<std::vector<std::shared_ptr<
const DirichletBC<T>>>>& bcs1,
103 const std::vector<Eigen::Ref<
const Eigen::Matrix<T, Eigen::Dynamic, 1>>>&
119 template <
typename T>
121 Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
const Form<T>& a,
122 const Eigen::Ref<
const Eigen::Matrix<T, Eigen::Dynamic, 1>>& bc_values1,
123 const std::vector<bool>& bc_markers1,
124 const Eigen::Ref<
const Eigen::Matrix<T, Eigen::Dynamic, 1>>& x0,
128 template <
typename T>
130 Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
131 const mesh::Geometry& geometry,
132 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
133 const std::uint8_t*,
const std::uint32_t)>& kernel,
134 const std::vector<std::int32_t>& active_cells,
135 const graph::AdjacencyList<std::int32_t>& dofmap0,
136 const graph::AdjacencyList<std::int32_t>& dofmap1,
137 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
139 const Eigen::Array<T, Eigen::Dynamic, 1>& constant_values,
140 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info,
141 const Eigen::Ref<
const Eigen::Matrix<T, Eigen::Dynamic, 1>>& bc_values1,
142 const std::vector<bool>& bc_markers1,
143 const Eigen::Ref<
const Eigen::Matrix<T, Eigen::Dynamic, 1>>& x0,
147 const int gdim = geometry.dim();
148 const graph::AdjacencyList<std::int32_t>& x_dofmap = geometry.dofmap();
151 const int num_dofs_g = x_dofmap.num_links(0);
152 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
156 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
157 coordinate_dofs(num_dofs_g, gdim);
158 Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> Ae;
159 Eigen::Matrix<T, Eigen::Dynamic, 1> be;
161 for (std::int32_t c : active_cells)
164 auto dmap1 = dofmap1.links(c);
168 for (Eigen::Index j = 0; j < dmap1.size(); ++j)
170 if (bc_markers1[dmap1[j]])
181 auto x_dofs = x_dofmap.links(c);
182 for (
int i = 0; i < num_dofs_g; ++i)
183 for (
int j = 0; j < gdim; ++j)
184 coordinate_dofs(i, j) = x_g(x_dofs[i], j);
187 auto dmap0 = dofmap0.links(c);
189 auto coeff_array = coeffs.row(c);
190 Ae.setZero(dmap0.size(), dmap1.size());
191 kernel(Ae.data(), coeff_array.data(), constant_values.data(),
192 coordinate_dofs.data(),
nullptr,
nullptr, cell_info[c]);
195 be.setZero(dmap0.size());
196 for (Eigen::Index j = 0; j < dmap1.size(); ++j)
198 const std::int32_t jj = dmap1[j];
201 const T bc = bc_values1[jj];
203 be -= Ae.col(j) * scale * (bc - x0[jj]);
205 be -= Ae.col(j) * scale * bc;
209 for (Eigen::Index k = 0; k < dmap0.size(); ++k)
210 b[dmap0[k]] += be[k];
214 template <
typename T>
215 void _lift_bc_exterior_facets(
216 Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
const mesh::Mesh& mesh,
217 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
218 const std::uint8_t*,
const std::uint32_t)>& kernel,
219 const std::vector<std::int32_t>& active_facets,
220 const graph::AdjacencyList<std::int32_t>& dofmap0,
221 const graph::AdjacencyList<std::int32_t>& dofmap1,
222 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
224 const Eigen::Array<T, Eigen::Dynamic, 1>& constant_values,
225 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info,
226 const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>& perms,
227 const Eigen::Ref<
const Eigen::Matrix<T, Eigen::Dynamic, 1>>& bc_values1,
228 const std::vector<bool>& bc_markers1,
229 const Eigen::Ref<
const Eigen::Matrix<T, Eigen::Dynamic, 1>>& x0,
232 const int gdim = mesh.geometry().dim();
233 const int tdim = mesh.topology().dim();
236 const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
238 const int num_dofs_g = x_dofmap.num_links(0);
239 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
240 = mesh.geometry().x();
243 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
244 coordinate_dofs(num_dofs_g, gdim);
245 Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> Ae;
246 Eigen::Matrix<T, Eigen::Dynamic, 1> 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 const auto* it = std::find(facets.data(), facets.data() + facets.rows(), f);
266 assert(it != (facets.data() + facets.rows()));
267 const int local_facet = std::distance(facets.data(), it);
270 auto dmap1 = dofmap1.links(cell);
274 for (Eigen::Index j = 0; j < dmap1.size(); ++j)
276 if (bc_markers1[dmap1[j]])
287 auto x_dofs = x_dofmap.links(cell);
288 for (
int i = 0; i < num_dofs_g; ++i)
289 for (
int j = 0; j < gdim; ++j)
290 coordinate_dofs(i, j) = x_g(x_dofs[i], j);
293 auto dmap0 = dofmap0.links(cell);
298 auto coeff_array = coeffs.row(cell);
299 Ae.setZero(dmap0.size(), dmap1.size());
300 kernel(Ae.data(), coeff_array.data(), constant_values.data(),
301 coordinate_dofs.data(), &local_facet, &perms(local_facet, cell),
305 be.setZero(dmap0.size());
306 for (Eigen::Index j = 0; j < dmap1.size(); ++j)
308 const std::int32_t jj = dmap1[j];
311 const T bc = bc_values1[jj];
313 be -= Ae.col(j) * scale * (bc - x0[jj]);
315 be -= Ae.col(j) * scale * bc;
319 for (Eigen::Index k = 0; k < dmap0.size(); ++k)
320 b[dmap0[k]] += be[k];
324 template <
typename T>
325 void _lift_bc_interior_facets(
326 Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> 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,
331 const graph::AdjacencyList<std::int32_t>& dofmap1,
332 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
334 const std::vector<int>& offsets,
335 const Eigen::Array<T, Eigen::Dynamic, 1>& constant_values,
336 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info,
337 const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>& perms,
338 const Eigen::Ref<
const Eigen::Matrix<T, Eigen::Dynamic, 1>>& bc_values1,
339 const std::vector<bool>& bc_markers1,
340 const Eigen::Ref<
const Eigen::Matrix<T, Eigen::Dynamic, 1>>& x0,
343 const int gdim = mesh.geometry().dim();
344 const int tdim = mesh.topology().dim();
347 const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
349 const int num_dofs_g = x_dofmap.num_links(0);
350 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
351 = mesh.geometry().x();
354 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
355 coordinate_dofs(2 * num_dofs_g, gdim);
356 Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> Ae;
357 Eigen::Array<T, Eigen::Dynamic, 1> coeff_array(2 * offsets.back());
358 assert(offsets.back() == coeffs.cols());
359 Eigen::Matrix<T, Eigen::Dynamic, 1> be;
362 Eigen::Array<std::int32_t, Eigen::Dynamic, 1> dmapjoint0, dmapjoint1;
364 const mesh::Topology& topology = mesh.topology();
365 auto connectivity = topology.connectivity(tdim - 1, tdim);
366 assert(connectivity);
367 auto c_to_f = topology.connectivity(tdim, tdim - 1);
369 auto f_to_c = topology.connectivity(tdim - 1, tdim);
371 auto map = topology.index_map(tdim - 1);
374 for (std::int32_t f : active_facets)
377 auto cells = f_to_c->links(f);
378 assert(
cells.rows() == 2);
381 auto facets0 = c_to_f->links(cells[0]);
383 = std::find(facets0.data(), facets0.data() + facets0.rows(), f);
384 assert(it0 != (facets0.data() + facets0.rows()));
385 const int local_facet0 = std::distance(facets0.data(), it0);
386 auto facets1 = c_to_f->links(cells[1]);
388 = std::find(facets1.data(), facets1.data() + facets1.rows(), f);
389 assert(it1 != (facets1.data() + facets1.rows()));
390 const int local_facet1 = std::distance(facets1.data(), it1);
392 const std::array local_facet{local_facet0, local_facet1};
395 auto x_dofs0 = x_dofmap.links(cells[0]);
396 auto x_dofs1 = x_dofmap.links(cells[1]);
397 for (
int i = 0; i < num_dofs_g; ++i)
399 for (
int j = 0; j < gdim; ++j)
401 coordinate_dofs(i, j) = x_g(x_dofs0[i], j);
402 coordinate_dofs(i + num_dofs_g, j) = x_g(x_dofs1[i], j);
407 auto dmap0_cell0 = dofmap0.links(cells[0]);
408 auto dmap0_cell1 = dofmap0.links(cells[1]);
409 dmapjoint0.resize(dmap0_cell0.size() + dmap0_cell1.size());
410 dmapjoint0.head(dmap0_cell0.size()) = dmap0_cell0;
411 dmapjoint0.tail(dmap0_cell1.size()) = dmap0_cell1;
413 auto dmap1_cell0 = dofmap1.links(cells[0]);
414 auto dmap1_cell1 = dofmap1.links(cells[1]);
415 dmapjoint1.resize(dmap1_cell0.size() + dmap1_cell1.size());
416 dmapjoint1.head(dmap1_cell0.size()) = dmap1_cell0;
417 dmapjoint1.tail(dmap1_cell1.size()) = dmap1_cell1;
421 for (Eigen::Index j = 0; j < dmap1_cell0.size(); ++j)
423 if (bc_markers1[dmap1_cell0[j]])
430 for (Eigen::Index j = 0; j < dmap1_cell1.size(); ++j)
432 if (bc_markers1[dmap1_cell1[j]])
444 auto coeff_cell0 = coeffs.row(cells[0]);
445 auto coeff_cell1 = coeffs.row(cells[1]);
448 for (std::size_t i = 0; i < offsets.size() - 1; ++i)
451 const int num_entries = offsets[i + 1] - offsets[i];
452 coeff_array.segment(2 * offsets[i], num_entries)
453 = coeff_cell0.segment(offsets[i], num_entries);
454 coeff_array.segment(offsets[i + 1] + offsets[i], num_entries)
455 = coeff_cell1.segment(offsets[i], num_entries);
459 Ae.setZero(dmapjoint0.size(), dmapjoint1.size());
460 const std::array perm{perms(local_facet[0], cells[0]),
461 perms(local_facet[1], cells[1])};
462 kernel(Ae.data(), coeff_array.data(), constant_values.data(),
463 coordinate_dofs.data(), local_facet.data(), perm.data(),
464 cell_info[cells[0]]);
467 be.setZero(dmap0_cell0.size() + dmap1_cell0.size());
468 for (Eigen::Index j = 0; j < dmap1_cell0.size(); ++j)
470 const std::int32_t jj = dmap1_cell0[j];
473 const T bc = bc_values1[jj];
475 be -= Ae.col(j) * scale * (bc - x0[jj]);
477 be -= Ae.col(j) * scale * bc;
481 for (Eigen::Index k = 0; k < dmap0_cell0.size(); ++k)
482 b[dmap0_cell0[k]] += be[k];
485 be.setZero(dmap0_cell1.size() + dmap1_cell1.size());
486 for (Eigen::Index j = 0; j < dmap1_cell1.size(); ++j)
488 const std::int32_t jj = dmap1_cell1[j];
491 const T bc = bc_values1[jj];
493 be -= Ae.col(j) * scale * (bc - x0[jj]);
495 be -= Ae.col(j) * scale * bc;
499 for (Eigen::Index k = 0; k < dmap0_cell1.size(); ++k)
500 b[dmap0_cell1[k]] += be[k];
504 template <
typename T>
505 void assemble_vector(Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
508 std::shared_ptr<const mesh::Mesh> mesh = L.mesh();
510 const int tdim = mesh->topology().dim();
511 const std::int32_t num_cells
512 = mesh->topology().connectivity(tdim, 0)->num_nodes();
515 assert(L.function_spaces().at(0));
516 std::shared_ptr<const fem::DofMap> dofmap
517 = L.function_spaces().at(0)->dofmap();
519 const graph::AdjacencyList<std::int32_t>& dofs = dofmap->list();
522 const Eigen::Array<T, Eigen::Dynamic, 1> constant_values =
pack_constants(L);
525 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> coeffs
528 const bool needs_permutation_data = L.needs_permutation_data();
529 if (needs_permutation_data)
530 mesh->topology_mutable().create_entity_permutations();
531 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info
532 = needs_permutation_data
533 ? mesh->topology().get_cell_permutation_info()
534 : Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>(num_cells);
536 for (
int i : L.integral_ids(IntegralType::cell))
538 const auto& fn = L.kernel(IntegralType::cell, i);
539 const std::vector<std::int32_t>& active_cells
540 = L.domains(IntegralType::cell, i);
541 fem::impl::assemble_cells(b, mesh->geometry(), active_cells, dofs, fn,
542 coeffs, constant_values, cell_info);
545 if (L.num_integrals(IntegralType::exterior_facet) > 0
546 or L.num_integrals(IntegralType::interior_facet) > 0)
549 mesh->topology_mutable().create_entities(tdim - 1);
550 mesh->topology_mutable().create_connectivity(tdim - 1, tdim);
552 const int facets_per_cell
554 const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>& perms
555 = needs_permutation_data
556 ? mesh->topology().get_facet_permutations()
557 : Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>(
558 facets_per_cell, num_cells);
559 for (
int i : L.integral_ids(IntegralType::exterior_facet))
561 const auto& fn = L.kernel(IntegralType::exterior_facet, i);
562 const std::vector<std::int32_t>& active_facets
563 = L.domains(IntegralType::exterior_facet, i);
564 fem::impl::assemble_exterior_facets(b, *mesh, active_facets, dofs, fn,
565 coeffs, constant_values, cell_info,
569 const std::vector<int> c_offsets = L.coefficient_offsets();
570 for (
int i : L.integral_ids(IntegralType::interior_facet))
572 const auto& fn = L.kernel(IntegralType::interior_facet, i);
573 const std::vector<std::int32_t>& active_facets
574 = L.domains(IntegralType::interior_facet, i);
575 fem::impl::assemble_interior_facets(b, *mesh, active_facets, *dofmap, fn,
576 coeffs, c_offsets, constant_values,
582 template <
typename T>
584 Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
585 const mesh::Geometry& geometry,
586 const std::vector<std::int32_t>& active_cells,
587 const graph::AdjacencyList<std::int32_t>& dofmap,
588 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
589 const std::uint8_t*,
const std::uint32_t)>& kernel,
590 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
592 const Eigen::Array<T, Eigen::Dynamic, 1>& constant_values,
593 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info)
595 const int gdim = geometry.dim();
598 const graph::AdjacencyList<std::int32_t>& x_dofmap = geometry.dofmap();
601 const int num_dofs_g = x_dofmap.num_links(0);
602 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
607 const int num_dofs = dofmap.links(0).size();
608 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
609 coordinate_dofs(num_dofs_g, gdim);
610 Eigen::Matrix<T, Eigen::Dynamic, 1> be(num_dofs);
613 for (std::int32_t c : active_cells)
616 auto x_dofs = x_dofmap.links(c);
617 for (
int i = 0; i < num_dofs_g; ++i)
618 for (
int j = 0; j < gdim; ++j)
619 coordinate_dofs(i, j) = x_g(x_dofs[i], j);
622 std::fill(be.data(), be.data() + num_dofs, 0);
623 kernel(be.data(), coeffs.row(c).data(), constant_values.data(),
624 coordinate_dofs.data(),
nullptr,
nullptr, cell_info[c]);
627 auto dofs = dofmap.links(c);
628 for (Eigen::Index i = 0; i < num_dofs; ++i)
633 template <
typename T>
634 void assemble_exterior_facets(
635 Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
const mesh::Mesh& mesh,
636 const std::vector<std::int32_t>& active_facets,
637 const graph::AdjacencyList<std::int32_t>& dofmap,
638 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
639 const std::uint8_t*,
const std::uint32_t)>& fn,
640 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
642 const Eigen::Array<T, Eigen::Dynamic, 1>& constant_values,
643 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info,
644 const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>& perms)
646 const int gdim = mesh.geometry().dim();
647 const int tdim = mesh.topology().dim();
650 const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
653 const int num_dofs_g = x_dofmap.num_links(0);
654 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
655 = mesh.geometry().x();
659 const int num_dofs = dofmap.links(0).size();
660 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
661 coordinate_dofs(num_dofs_g, gdim);
662 Eigen::Matrix<T, Eigen::Dynamic, 1> be(num_dofs);
664 auto f_to_c = mesh.topology().connectivity(tdim - 1, tdim);
666 auto c_to_f = mesh.topology().connectivity(tdim, tdim - 1);
668 for (
const auto& f : active_facets)
671 assert(f_to_c->num_links(f) > 0);
672 const std::int32_t cell = f_to_c->links(f)[0];
675 auto facets = c_to_f->links(cell);
676 const auto* it = std::find(facets.data(), facets.data() + facets.rows(), f);
677 assert(it != (facets.data() + facets.rows()));
678 const int local_facet = std::distance(facets.data(), it);
681 auto x_dofs = x_dofmap.links(cell);
682 for (
int i = 0; i < num_dofs_g; ++i)
683 for (
int j = 0; j < gdim; ++j)
684 coordinate_dofs(i, j) = x_g(x_dofs[i], j);
687 std::fill(be.data(), be.data() + num_dofs, 0);
688 fn(be.data(), coeffs.row(cell).data(), constant_values.data(),
689 coordinate_dofs.data(), &local_facet, &perms(local_facet, cell),
693 auto dofs = dofmap.links(cell);
694 for (Eigen::Index i = 0; i < num_dofs; ++i)
699 template <
typename T>
700 void assemble_interior_facets(
701 Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
const mesh::Mesh& mesh,
702 const std::vector<std::int32_t>& active_facets,
const fem::DofMap& dofmap,
703 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
704 const std::uint8_t*,
const std::uint32_t)>& fn,
705 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
707 const std::vector<int>& offsets,
708 const Eigen::Array<T, Eigen::Dynamic, 1>& constant_values,
709 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info,
710 const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>& perms)
712 const int gdim = mesh.geometry().dim();
713 const int tdim = mesh.topology().dim();
716 const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
719 const int num_dofs_g = x_dofmap.num_links(0);
720 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
721 = mesh.geometry().x();
724 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
725 coordinate_dofs(2 * num_dofs_g, gdim);
726 Eigen::Matrix<T, Eigen::Dynamic, 1> be;
727 Eigen::Array<T, Eigen::Dynamic, 1> coeff_array(2 * offsets.back());
728 assert(offsets.back() == coeffs.cols());
730 auto f_to_c = mesh.topology().connectivity(tdim - 1, tdim);
732 auto c_to_f = mesh.topology().connectivity(tdim, tdim - 1);
734 for (
const auto& f : active_facets)
737 auto cells = f_to_c->links(f);
738 assert(
cells.rows() == 2);
741 std::array<int, 2> local_facet;
742 for (
int i = 0; i < 2; ++i)
744 auto facets = c_to_f->links(cells[i]);
746 = std::find(facets.data(), facets.data() + facets.rows(), f);
747 assert(it != (facets.data() + facets.rows()));
748 local_facet[i] = std::distance(facets.data(), it);
752 auto x_dofs0 = x_dofmap.links(cells[0]);
753 auto x_dofs1 = x_dofmap.links(cells[1]);
754 for (
int i = 0; i < num_dofs_g; ++i)
756 for (
int j = 0; j < gdim; ++j)
758 coordinate_dofs(i, j) = x_g(x_dofs0[i], j);
759 coordinate_dofs(i + num_dofs_g, j) = x_g(x_dofs1[i], j);
764 auto dmap0 = dofmap.cell_dofs(cells[0]);
765 auto dmap1 = dofmap.cell_dofs(cells[1]);
769 auto coeff_cell0 = coeffs.row(cells[0]);
770 auto coeff_cell1 = coeffs.row(cells[1]);
773 for (std::size_t i = 0; i < offsets.size() - 1; ++i)
776 const int num_entries = offsets[i + 1] - offsets[i];
777 coeff_array.segment(2 * offsets[i], num_entries)
778 = coeff_cell0.segment(offsets[i], num_entries);
779 coeff_array.segment(offsets[i + 1] + offsets[i], num_entries)
780 = coeff_cell1.segment(offsets[i], num_entries);
784 be.setZero(dmap0.size() + dmap1.size());
786 const std::array perm{perms(local_facet[0], cells[0]),
787 perms(local_facet[1], cells[1])};
788 fn(be.data(), coeff_array.data(), constant_values.data(),
789 coordinate_dofs.data(), local_facet.data(), perm.data(),
790 cell_info[cells[0]]);
793 for (Eigen::Index i = 0; i < dmap0.size(); ++i)
794 b[dmap0[i]] += be[i];
795 for (Eigen::Index i = 0; i < dmap1.size(); ++i)
796 b[dmap1[i]] += be[i + dmap0.size()];
800 template <
typename T>
802 Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
803 const std::vector<std::shared_ptr<
const Form<T>>> a,
804 const std::vector<std::vector<std::shared_ptr<
const DirichletBC<T>>>>& bcs1,
805 const std::vector<Eigen::Ref<
const Eigen::Matrix<T, Eigen::Dynamic, 1>>>&
810 if (!x0.empty() and x0.size() != a.size())
812 throw std::runtime_error(
813 "Mismatch in size between x0 and bilinear form in assembler.");
816 if (a.size() != bcs1.size())
818 throw std::runtime_error(
819 "Mismatch in size between a and bcs in assembler.");
822 for (std::size_t j = 0; j < a.size(); ++j)
824 std::vector<bool> bc_markers1;
825 Eigen::Matrix<T, Eigen::Dynamic, 1> bc_values1;
826 if (a[j] and !bcs1[j].empty())
828 auto V1 = a[j]->function_spaces()[1];
830 auto map1 = V1->dofmap()->index_map;
833 = map1->block_size() * (map1->size_local() + map1->num_ghosts());
834 bc_markers1.assign(crange,
false);
835 bc_values1 = Eigen::Matrix<T, Eigen::Dynamic, 1>::Zero(crange);
836 for (
const std::shared_ptr<
const DirichletBC<T>>& bc : bcs1[j])
838 bc->mark_dofs(bc_markers1);
839 bc->dof_values(bc_values1);
843 lift_bc<T>(b, *a[j], bc_values1, bc_markers1, x0[j], scale);
846 const Eigen::Matrix<T, Eigen::Dynamic, 1> x0(0);
847 lift_bc<T>(b, *a[j], bc_values1, bc_markers1, x0, scale);
853 template <
typename T>
855 Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
const Form<T>& a,
856 const Eigen::Ref<
const Eigen::Matrix<T, Eigen::Dynamic, 1>>& bc_values1,
857 const std::vector<bool>& bc_markers1,
858 const Eigen::Ref<
const Eigen::Matrix<T, Eigen::Dynamic, 1>>& x0,
861 std::shared_ptr<const mesh::Mesh> mesh = a.mesh();
863 const int tdim = mesh->topology().dim();
864 const std::int32_t num_cells
865 = mesh->topology().connectivity(tdim, 0)->num_nodes();
868 assert(a.function_spaces().at(0));
869 assert(a.function_spaces().at(1));
870 const graph::AdjacencyList<std::int32_t>& dofmap0
871 = a.function_spaces()[0]->dofmap()->list();
872 const graph::AdjacencyList<std::int32_t>& dofmap1
873 = a.function_spaces()[1]->dofmap()->list();
876 const Eigen::Array<T, Eigen::Dynamic, 1> constant_values =
pack_constants(a);
879 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> coeffs
882 const bool needs_permutation_data = a.needs_permutation_data();
883 if (needs_permutation_data)
884 mesh->topology_mutable().create_entity_permutations();
885 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info
886 = needs_permutation_data
887 ? mesh->topology().get_cell_permutation_info()
888 : Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>(num_cells);
890 for (
int i : a.integral_ids(IntegralType::cell))
892 const auto& kernel = a.kernel(IntegralType::cell, i);
893 const std::vector<std::int32_t>& active_cells
894 = a.domains(IntegralType::cell, i);
895 _lift_bc_cells(b, mesh->geometry(), kernel, active_cells, dofmap0, dofmap1,
896 coeffs, constant_values, cell_info, bc_values1, bc_markers1,
900 if (a.num_integrals(IntegralType::exterior_facet) > 0
901 or a.num_integrals(IntegralType::interior_facet) > 0)
904 mesh->topology_mutable().create_entities(tdim - 1);
905 mesh->topology_mutable().create_connectivity(tdim - 1, tdim);
907 const int facets_per_cell
909 const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>& perms
910 = needs_permutation_data
911 ? mesh->topology().get_facet_permutations()
912 : Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>(
913 facets_per_cell, num_cells);
914 for (
int i : a.integral_ids(IntegralType::exterior_facet))
916 const auto& kernel = a.kernel(IntegralType::exterior_facet, i);
917 const std::vector<std::int32_t>& active_facets
918 = a.domains(IntegralType::exterior_facet, i);
919 _lift_bc_exterior_facets(b, *mesh, kernel, active_facets, dofmap0,
920 dofmap1, coeffs, constant_values, cell_info,
921 perms, bc_values1, bc_markers1, x0, scale);
924 const std::vector<int> c_offsets = a.coefficient_offsets();
925 for (
int i : a.integral_ids(IntegralType::interior_facet))
927 const auto& kernel = a.kernel(IntegralType::interior_facet, i);
928 const std::vector<std::int32_t>& active_facets
929 = a.domains(IntegralType::interior_facet, i);
930 _lift_bc_interior_facets(b, *mesh, kernel, active_facets, dofmap0,
931 dofmap1, coeffs, c_offsets, constant_values,
932 cell_info, perms, bc_values1, bc_markers1, x0,
void cells(la::SparsityPattern &pattern, const mesh::Topology &topology, const std::array< const fem::DofMap *, 2 > dofmaps)
Iterate over cells and insert entries into sparsity pattern.
Definition: SparsityPatternBuilder.cpp:18
Eigen::Array< typename U::scalar_type, Eigen::Dynamic, 1 > pack_constants(const U &u)
Pack constants of u of generic type U ready for assembly.
Definition: utils.h:456
void apply_lifting(Eigen::Ref< Eigen::Matrix< T, Eigen::Dynamic, 1 >> 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< Eigen::Ref< const Eigen::Matrix< T, Eigen::Dynamic, 1 >>> &x0, double scale)
Modify b such that:
Definition: assembler.h:75
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:407
void assemble_vector(Eigen::Ref< Eigen::Matrix< T, Eigen::Dynamic, 1 >> b, const Form< T > &L)
Assemble linear form into an Eigen vector.
Definition: assembler.h:50
int cell_num_entities(mesh::CellType type, int dim)
Number of entities of dimension dim.
Definition: cell_types.cpp:381