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>>>&
117 template <
typename T>
119 Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
const Form<T>& a,
120 const Eigen::Ref<
const Eigen::Matrix<T, Eigen::Dynamic, 1>>& bc_values1,
121 const std::vector<bool>& bc_markers1,
double scale);
135 template <
typename T>
137 Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
const Form<T>& a,
138 const Eigen::Ref<
const Eigen::Matrix<T, Eigen::Dynamic, 1>>& bc_values1,
139 const std::vector<bool>& bc_markers1,
140 const Eigen::Ref<
const Eigen::Matrix<T, Eigen::Dynamic, 1>>& x0,
144 template <
typename T>
146 Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
const Form<T>& a,
147 const Eigen::Ref<
const Eigen::Matrix<T, Eigen::Dynamic, 1>>& bc_values1,
148 const std::vector<bool>& bc_markers1,
149 const Eigen::Ref<
const Eigen::Matrix<T, Eigen::Dynamic, 1>>& x0,
152 assert(a.rank() == 2);
155 std::shared_ptr<const mesh::Mesh> mesh = a.mesh();
158 mesh->topology_mutable().create_entity_permutations();
161 assert(a.function_space(0));
162 assert(a.function_space(1));
163 std::shared_ptr<const fem::DofMap> dofmap0 = a.function_space(0)->dofmap();
164 std::shared_ptr<const fem::DofMap> dofmap1 = a.function_space(1)->dofmap();
169 Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> coeffs
172 const std::function<void(T*,
const T*,
const T*,
const double*,
const int*,
173 const std::uint8_t*,
const std::uint32_t)>& fn
174 = a.integrals().get_tabulate_tensor(IntegralType::cell, 0);
177 const int gdim = mesh->geometry().dim();
178 const graph::AdjacencyList<std::int32_t>& x_dofmap
179 = mesh->geometry().dofmap();
182 const int num_dofs_g = x_dofmap.num_links(0);
183 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
184 = mesh->geometry().x();
187 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
188 coordinate_dofs(num_dofs_g, gdim);
189 Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> Ae;
190 Eigen::Matrix<T, Eigen::Dynamic, 1> be;
193 if (!a.all_constants_set())
194 throw std::runtime_error(
"Unset constant in Form");
195 const Eigen::Array<T, Eigen::Dynamic, 1> constant_values =
pack_constants(a);
197 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info
198 = mesh->topology().get_cell_permutation_info();
201 const int tdim = mesh->topology().dim();
202 auto map = mesh->topology().index_map(tdim);
204 const int num_cells = map->size_local();
205 for (
int c = 0; c < num_cells; ++c)
208 auto dmap1 = dofmap1->cell_dofs(c);
212 for (Eigen::Index j = 0; j < dmap1.size(); ++j)
214 if (bc_markers1[dmap1[j]])
225 auto x_dofs = x_dofmap.links(c);
226 for (
int i = 0; i < num_dofs_g; ++i)
227 for (
int j = 0; j < gdim; ++j)
228 coordinate_dofs(i, j) = x_g(x_dofs[i], j);
231 auto dmap0 = dofmap0->cell_dofs(c);
233 auto coeff_array = coeffs.row(c);
234 Ae.setZero(dmap0.size(), dmap1.size());
235 fn(Ae.data(), coeff_array.data(), constant_values.data(),
236 coordinate_dofs.data(),
nullptr,
nullptr, cell_info[c]);
239 be.setZero(dmap0.size());
240 for (Eigen::Index j = 0; j < dmap1.size(); ++j)
242 const std::int32_t jj = dmap1[j];
245 const T bc = bc_values1[jj];
247 be -= Ae.col(j) * scale * (bc - x0[jj]);
249 be -= Ae.col(j) * scale * bc;
253 for (Eigen::Index k = 0; k < dmap0.size(); ++k)
254 b[dmap0[k]] += be[k];
258 template <
typename T>
259 void _lift_bc_exterior_facets(
260 Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
const Form<T>& a,
261 const Eigen::Ref<
const Eigen::Matrix<T, Eigen::Dynamic, 1>>& bc_values1,
262 const std::vector<bool>& bc_markers1,
263 const Eigen::Ref<
const Eigen::Matrix<T, Eigen::Dynamic, 1>>& x0,
266 assert(a.rank() == 2);
269 std::shared_ptr<const mesh::Mesh> mesh = a.mesh();
272 mesh->topology_mutable().create_entity_permutations();
274 const int gdim = mesh->geometry().dim();
275 const int tdim = mesh->topology().dim();
278 mesh->topology_mutable().create_entities(tdim - 1);
279 mesh->topology_mutable().create_connectivity(tdim - 1, tdim);
281 mesh->topology_mutable().create_entity_permutations();
284 assert(a.function_space(0));
285 assert(a.function_space(1));
286 std::shared_ptr<const fem::DofMap> dofmap0 = a.function_space(0)->dofmap();
287 std::shared_ptr<const fem::DofMap> dofmap1 = a.function_space(1)->dofmap();
292 Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> coeffs
295 const std::function<void(T*,
const T*,
const T*,
const double*,
const int*,
296 const std::uint8_t*,
const std::uint32_t)>& fn
297 = a.integrals().get_tabulate_tensor(IntegralType::exterior_facet, 0);
300 const graph::AdjacencyList<std::int32_t>& x_dofmap
301 = mesh->geometry().dofmap();
303 const int num_dofs_g = x_dofmap.num_links(0);
304 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
305 = mesh->geometry().x();
308 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
309 coordinate_dofs(num_dofs_g, gdim);
310 Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> Ae;
311 Eigen::Matrix<T, Eigen::Dynamic, 1> be;
314 if (!a.all_constants_set())
315 throw std::runtime_error(
"Unset constant in Form");
316 const Eigen::Array<T, Eigen::Dynamic, 1> constant_values =
pack_constants(a);
319 const mesh::Topology& topology = mesh->topology();
320 auto connectivity = topology.connectivity(tdim - 1, tdim);
321 assert(connectivity);
322 auto c_to_f = topology.connectivity(tdim, tdim - 1);
324 auto map = topology.index_map(tdim - 1);
327 const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>& perms
328 = mesh->topology().get_facet_permutations();
329 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info
330 = mesh->topology().get_cell_permutation_info();
332 std::set<std::int32_t> fwd_shared_facets;
334 if (topology.index_map(tdim)->num_ghosts() == 0)
336 fwd_shared_facets.insert(
337 topology.index_map(tdim - 1)->shared_indices().begin(),
338 topology.index_map(tdim - 1)->shared_indices().end());
341 for (
int f = 0; f < map->size_local(); ++f)
346 if (connectivity->num_links(f) == 2
347 or fwd_shared_facets.find(f) != fwd_shared_facets.end())
351 assert(connectivity->num_links(f) == 1);
352 const std::int32_t cell = connectivity->links(f)[0];
355 auto facets = c_to_f->links(cell);
356 const auto* it = std::find(facets.data(), facets.data() + facets.rows(), f);
357 assert(it != (facets.data() + facets.rows()));
358 const int local_facet = std::distance(facets.data(), it);
360 const std::uint8_t perm = perms(local_facet, cell);
363 auto dmap1 = dofmap1->cell_dofs(cell);
367 for (Eigen::Index j = 0; j < dmap1.size(); ++j)
369 if (bc_markers1[dmap1[j]])
380 auto x_dofs = x_dofmap.links(cell);
381 for (
int i = 0; i < num_dofs_g; ++i)
382 for (
int j = 0; j < gdim; ++j)
383 coordinate_dofs(i, j) = x_g(x_dofs[i], j);
386 auto dmap0 = dofmap0->cell_dofs(cell);
391 auto coeff_array = coeffs.row(cell);
392 Ae.setZero(dmap0.size(), dmap1.size());
393 fn(Ae.data(), coeff_array.data(), constant_values.data(),
394 coordinate_dofs.data(), &local_facet, &perm, cell_info[cell]);
397 be.setZero(dmap0.size());
398 for (Eigen::Index j = 0; j < dmap1.size(); ++j)
400 const std::int32_t jj = dmap1[j];
403 const T bc = bc_values1[jj];
405 be -= Ae.col(j) * scale * (bc - x0[jj]);
407 be -= Ae.col(j) * scale * bc;
411 for (Eigen::Index k = 0; k < dmap0.size(); ++k)
412 b[dmap0[k]] += be[k];
417 template <
typename T>
418 void assemble_vector(Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
421 std::shared_ptr<const mesh::Mesh> mesh = L.mesh();
423 const int tdim = mesh->topology().dim();
424 const std::int32_t num_cells
425 = mesh->topology().connectivity(tdim, 0)->num_nodes();
428 assert(L.function_space(0));
429 std::shared_ptr<const fem::DofMap> dofmap = L.function_space(0)->dofmap();
431 const graph::AdjacencyList<std::int32_t>& dofs = dofmap->list();
434 if (!L.all_constants_set())
435 throw std::runtime_error(
"Unset constant in Form");
436 const Eigen::Array<T, Eigen::Dynamic, 1> constant_values =
pack_constants(L);
439 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> coeffs
442 const FormIntegrals<T>& integrals = L.integrals();
443 const bool needs_permutation_data = integrals.needs_permutation_data();
444 if (needs_permutation_data)
445 mesh->topology_mutable().create_entity_permutations();
446 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info
447 = needs_permutation_data
448 ? mesh->topology().get_cell_permutation_info()
449 : Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>(num_cells);
451 for (
int i = 0; i < integrals.num_integrals(IntegralType::cell); ++i)
453 const auto& fn = integrals.get_tabulate_tensor(IntegralType::cell, i);
454 const std::vector<std::int32_t>& active_cells
455 = integrals.integral_domains(IntegralType::cell, i);
456 fem::impl::assemble_cells(b, mesh->geometry(), active_cells, dofs, fn,
457 coeffs, constant_values, cell_info);
460 if (integrals.num_integrals(IntegralType::exterior_facet) > 0
461 or integrals.num_integrals(IntegralType::interior_facet) > 0)
464 mesh->topology_mutable().create_entities(tdim - 1);
465 mesh->topology_mutable().create_connectivity(tdim - 1, tdim);
467 const int facets_per_cell
469 const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>& perms
470 = needs_permutation_data
471 ? mesh->topology().get_facet_permutations()
472 : Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>(
473 facets_per_cell, num_cells);
474 for (
int i = 0; i < integrals.num_integrals(IntegralType::exterior_facet);
478 = integrals.get_tabulate_tensor(IntegralType::exterior_facet, i);
479 const std::vector<std::int32_t>& active_facets
480 = integrals.integral_domains(IntegralType::exterior_facet, i);
481 fem::impl::assemble_exterior_facets(b, *mesh, active_facets, dofs, fn,
482 coeffs, constant_values, cell_info,
486 const std::vector<int> c_offsets = L.coefficients().offsets();
487 for (
int i = 0; i < integrals.num_integrals(IntegralType::interior_facet);
491 = integrals.get_tabulate_tensor(IntegralType::interior_facet, i);
492 const std::vector<std::int32_t>& active_facets
493 = integrals.integral_domains(IntegralType::interior_facet, i);
494 fem::impl::assemble_interior_facets(b, *mesh, active_facets, *dofmap, fn,
495 coeffs, c_offsets, constant_values,
501 template <
typename T>
503 Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
504 const mesh::Geometry& geometry,
505 const std::vector<std::int32_t>& active_cells,
506 const graph::AdjacencyList<std::int32_t>& dofmap,
507 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
508 const std::uint8_t*,
const std::uint32_t)>& kernel,
509 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
511 const Eigen::Array<T, Eigen::Dynamic, 1>& constant_values,
512 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info)
514 const int gdim = geometry.dim();
517 const graph::AdjacencyList<std::int32_t>& x_dofmap = geometry.dofmap();
520 const int num_dofs_g = x_dofmap.num_links(0);
521 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
526 const int num_dofs = dofmap.links(0).size();
527 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
528 coordinate_dofs(num_dofs_g, gdim);
529 Eigen::Matrix<T, Eigen::Dynamic, 1> be(num_dofs);
532 for (std::int32_t c : active_cells)
535 auto x_dofs = x_dofmap.links(c);
536 for (
int i = 0; i < num_dofs_g; ++i)
537 for (
int j = 0; j < gdim; ++j)
538 coordinate_dofs(i, j) = x_g(x_dofs[i], j);
541 std::fill(be.data(), be.data() + num_dofs, 0);
542 kernel(be.data(), coeffs.row(c).data(), constant_values.data(),
543 coordinate_dofs.data(),
nullptr,
nullptr, cell_info[c]);
546 auto dofs = dofmap.links(c);
547 for (Eigen::Index i = 0; i < num_dofs; ++i)
552 template <
typename T>
553 void assemble_exterior_facets(
554 Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
const mesh::Mesh& mesh,
555 const std::vector<std::int32_t>& active_facets,
556 const graph::AdjacencyList<std::int32_t>& dofmap,
557 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
558 const std::uint8_t*,
const std::uint32_t)>& fn,
559 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
561 const Eigen::Array<T, Eigen::Dynamic, 1>& constant_values,
562 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info,
563 const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>& perms)
565 const int gdim = mesh.geometry().dim();
566 const int tdim = mesh.topology().dim();
569 const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
572 const int num_dofs_g = x_dofmap.num_links(0);
573 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
574 = mesh.geometry().x();
578 const int num_dofs = dofmap.links(0).size();
579 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
580 coordinate_dofs(num_dofs_g, gdim);
581 Eigen::Matrix<T, Eigen::Dynamic, 1> be(num_dofs);
583 auto f_to_c = mesh.topology().connectivity(tdim - 1, tdim);
585 auto c_to_f = mesh.topology().connectivity(tdim, tdim - 1);
587 for (
const auto& f : active_facets)
590 assert(f_to_c->num_links(f) > 0);
591 const std::int32_t cell = f_to_c->links(f)[0];
594 auto facets = c_to_f->links(cell);
595 const auto* it = std::find(facets.data(), facets.data() + facets.rows(), f);
596 assert(it != (facets.data() + facets.rows()));
597 const int local_facet = std::distance(facets.data(), it);
600 auto x_dofs = x_dofmap.links(cell);
601 for (
int i = 0; i < num_dofs_g; ++i)
602 for (
int j = 0; j < gdim; ++j)
603 coordinate_dofs(i, j) = x_g(x_dofs[i], j);
606 std::fill(be.data(), be.data() + num_dofs, 0);
607 fn(be.data(), coeffs.row(cell).data(), constant_values.data(),
608 coordinate_dofs.data(), &local_facet, &perms(local_facet, cell),
612 auto dofs = dofmap.links(cell);
613 for (Eigen::Index i = 0; i < num_dofs; ++i)
618 template <
typename T>
619 void assemble_interior_facets(
620 Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
const mesh::Mesh& mesh,
621 const std::vector<std::int32_t>& active_facets,
const fem::DofMap& dofmap,
622 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
623 const std::uint8_t*,
const std::uint32_t)>& fn,
624 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
626 const std::vector<int>& offsets,
627 const Eigen::Array<T, Eigen::Dynamic, 1>& constant_values,
628 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info,
629 const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>& perms)
631 const int gdim = mesh.geometry().dim();
632 const int tdim = mesh.topology().dim();
635 const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
638 const int num_dofs_g = x_dofmap.num_links(0);
639 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
640 = mesh.geometry().x();
643 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
644 coordinate_dofs(2 * num_dofs_g, gdim);
645 Eigen::Matrix<T, Eigen::Dynamic, 1> be;
646 Eigen::Array<T, Eigen::Dynamic, 1> coeff_array(2 * offsets.back());
647 assert(offsets.back() == coeffs.cols());
649 auto f_to_c = mesh.topology().connectivity(tdim - 1, tdim);
651 auto c_to_f = mesh.topology().connectivity(tdim, tdim - 1);
653 for (
const auto& f : active_facets)
656 auto cells = f_to_c->links(f);
657 assert(
cells.rows() == 2);
660 std::array<int, 2> local_facet;
661 for (
int i = 0; i < 2; ++i)
663 auto facets = c_to_f->links(cells[i]);
665 = std::find(facets.data(), facets.data() + facets.rows(), f);
666 assert(it != (facets.data() + facets.rows()));
667 local_facet[i] = std::distance(facets.data(), it);
671 auto x_dofs0 = x_dofmap.links(cells[0]);
672 auto x_dofs1 = x_dofmap.links(cells[1]);
673 for (
int i = 0; i < num_dofs_g; ++i)
675 for (
int j = 0; j < gdim; ++j)
677 coordinate_dofs(i, j) = x_g(x_dofs0[i], j);
678 coordinate_dofs(i + num_dofs_g, j) = x_g(x_dofs1[i], j);
683 auto dmap0 = dofmap.cell_dofs(cells[0]);
684 auto dmap1 = dofmap.cell_dofs(cells[1]);
688 auto coeff_cell0 = coeffs.row(cells[0]);
689 auto coeff_cell1 = coeffs.row(cells[1]);
692 for (std::size_t i = 0; i < offsets.size() - 1; ++i)
695 const int num_entries = offsets[i + 1] - offsets[i];
696 coeff_array.segment(2 * offsets[i], num_entries)
697 = coeff_cell0.segment(offsets[i], num_entries);
698 coeff_array.segment(offsets[i + 1] + offsets[i], num_entries)
699 = coeff_cell1.segment(offsets[i], num_entries);
703 be.setZero(dmap0.size() + dmap1.size());
705 const std::array perm{perms(local_facet[0], cells[0]),
706 perms(local_facet[1], cells[1])};
707 fn(be.data(), coeff_array.data(), constant_values.data(),
708 coordinate_dofs.data(), local_facet.data(), perm.data(),
709 cell_info[cells[0]]);
712 for (Eigen::Index i = 0; i < dmap0.size(); ++i)
713 b[dmap0[i]] += be[i];
714 for (Eigen::Index i = 0; i < dmap1.size(); ++i)
715 b[dmap1[i]] += be[i + dmap0.size()];
719 template <
typename T>
721 Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
722 const std::vector<std::shared_ptr<
const Form<T>>> a,
723 const std::vector<std::vector<std::shared_ptr<
const DirichletBC<T>>>>& bcs1,
724 const std::vector<Eigen::Ref<
const Eigen::Matrix<T, Eigen::Dynamic, 1>>>&
729 if (!x0.empty() and x0.size() != a.size())
731 throw std::runtime_error(
732 "Mismatch in size between x0 and bilinear form in assembler.");
735 if (a.size() != bcs1.size())
737 throw std::runtime_error(
738 "Mismatch in size between a and bcs in assembler.");
741 for (std::size_t j = 0; j < a.size(); ++j)
743 std::vector<bool> bc_markers1;
744 Eigen::Matrix<T, Eigen::Dynamic, 1> bc_values1;
745 if (a[j] and !bcs1[j].empty())
747 auto V1 = a[j]->function_space(1);
749 auto map1 = V1->dofmap()->index_map;
752 = map1->block_size() * (map1->size_local() + map1->num_ghosts());
753 bc_markers1.assign(crange,
false);
754 bc_values1 = Eigen::Matrix<T, Eigen::Dynamic, 1>::Zero(crange);
755 for (
const std::shared_ptr<
const DirichletBC<T>>& bc : bcs1[j])
757 bc->mark_dofs(bc_markers1);
758 bc->dof_values(bc_values1);
763 lift_bc<T>(b, *a[j], bc_values1, bc_markers1, x0[j], scale);
765 lift_bc<T>(b, *a[j], bc_values1, bc_markers1, scale);
770 template <
typename T>
772 Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
const Form<T>& a,
773 const Eigen::Ref<
const Eigen::Matrix<T, Eigen::Dynamic, 1>>& bc_values1,
774 const std::vector<bool>& bc_markers1,
double scale)
776 const Eigen::Matrix<T, Eigen::Dynamic, 1> x0(0);
777 if (a.integrals().num_integrals(fem::IntegralType::cell) > 0)
778 _lift_bc_cells<T>(b, a, bc_values1, bc_markers1, x0, scale);
779 if (a.integrals().num_integrals(fem::IntegralType::exterior_facet) > 0)
780 _lift_bc_exterior_facets<T>(b, a, bc_values1, bc_markers1, x0, scale);
783 template <
typename T>
785 Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
const Form<T>& a,
786 const Eigen::Ref<
const Eigen::Matrix<T, Eigen::Dynamic, 1>>& bc_values1,
787 const std::vector<bool>& bc_markers1,
788 const Eigen::Ref<
const Eigen::Matrix<T, Eigen::Dynamic, 1>>& x0,
791 if (a.integrals().num_integrals(fem::IntegralType::cell) > 0)
792 _lift_bc_cells(b, a, bc_values1, bc_markers1, x0, scale);
793 if (a.integrals().num_integrals(fem::IntegralType::exterior_facet) > 0)
794 _lift_bc_exterior_facets(b, a, bc_values1, bc_markers1, x0, scale);