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,
const mesh::Mesh& mesh,
43 const std::vector<std::int32_t>& active_cells,
44 const graph::AdjacencyList<std::int32_t>& dofmap,
45 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
46 const std::uint8_t*,
const std::uint32_t)>& kernel,
47 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
49 const Eigen::Array<T, Eigen::Dynamic, 1>& constant_values);
53 void assemble_exterior_facets(
54 Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
const mesh::Mesh& mesh,
55 const std::vector<std::int32_t>& active_facets,
56 const graph::AdjacencyList<std::int32_t>& dofmap,
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 Eigen::Array<T, Eigen::Dynamic, 1>& constant_values);
65 void assemble_interior_facets(
66 Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
const mesh::Mesh& mesh,
67 const std::vector<std::int32_t>& active_facets,
const fem::DofMap& dofmap,
68 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
69 const std::uint8_t*,
const std::uint32_t)>& fn,
70 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
72 const std::vector<int>& offsets,
73 const Eigen::Array<T, Eigen::Dynamic, 1>& constant_values);
94 Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
95 const std::vector<std::shared_ptr<
const Form<T>>> a,
96 const std::vector<std::vector<std::shared_ptr<
const DirichletBC<T>>>>& bcs1,
97 const std::vector<Eigen::Ref<
const Eigen::Matrix<T, Eigen::Dynamic, 1>>>&
111 template <
typename T>
113 Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
const Form<T>& a,
114 const Eigen::Ref<
const Eigen::Matrix<T, Eigen::Dynamic, 1>>& bc_values1,
115 const std::vector<bool>& bc_markers1,
double scale);
129 template <
typename T>
131 Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
const Form<T>& a,
132 const Eigen::Ref<
const Eigen::Matrix<T, Eigen::Dynamic, 1>>& bc_values1,
133 const std::vector<bool>& bc_markers1,
134 const Eigen::Ref<
const Eigen::Matrix<T, Eigen::Dynamic, 1>>& x0,
138 template <
typename T>
140 Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
const Form<T>& a,
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,
146 assert(a.rank() == 2);
149 std::shared_ptr<const mesh::Mesh> mesh = a.mesh();
152 mesh->topology_mutable().create_entity_permutations();
155 assert(a.function_space(0));
156 assert(a.function_space(1));
157 std::shared_ptr<const fem::DofMap> dofmap0 = a.function_space(0)->dofmap();
158 std::shared_ptr<const fem::DofMap> dofmap1 = a.function_space(1)->dofmap();
163 Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> coeffs
166 const std::function<void(T*,
const T*,
const T*,
const double*,
const int*,
167 const std::uint8_t*,
const std::uint32_t)>& fn
168 = a.integrals().get_tabulate_tensor(IntegralType::cell, 0);
171 const int gdim = mesh->geometry().dim();
172 const graph::AdjacencyList<std::int32_t>& x_dofmap
173 = mesh->geometry().dofmap();
176 const int num_dofs_g = x_dofmap.num_links(0);
177 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
178 = mesh->geometry().x();
181 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
182 coordinate_dofs(num_dofs_g, gdim);
183 Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> Ae;
184 Eigen::Matrix<T, Eigen::Dynamic, 1> be;
187 if (!a.all_constants_set())
188 throw std::runtime_error(
"Unset constant in Form");
189 const Eigen::Array<T, Eigen::Dynamic, 1> constant_values =
pack_constants(a);
191 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info
192 = mesh->topology().get_cell_permutation_info();
195 const int tdim = mesh->topology().dim();
196 auto map = mesh->topology().index_map(tdim);
198 const int num_cells = map->size_local();
199 for (
int c = 0; c < num_cells; ++c)
202 auto dmap1 = dofmap1->cell_dofs(c);
206 for (Eigen::Index j = 0; j < dmap1.size(); ++j)
208 if (bc_markers1[dmap1[j]])
219 auto x_dofs = x_dofmap.links(c);
220 for (
int i = 0; i < num_dofs_g; ++i)
221 for (
int j = 0; j < gdim; ++j)
222 coordinate_dofs(i, j) = x_g(x_dofs[i], j);
225 auto dmap0 = dofmap0->cell_dofs(c);
227 auto coeff_array = coeffs.row(c);
228 Ae.setZero(dmap0.size(), dmap1.size());
229 fn(Ae.data(), coeff_array.data(), constant_values.data(),
230 coordinate_dofs.data(),
nullptr,
nullptr, cell_info[c]);
233 be.setZero(dmap0.size());
234 for (Eigen::Index j = 0; j < dmap1.size(); ++j)
236 const std::int32_t jj = dmap1[j];
239 const T bc = bc_values1[jj];
241 be -= Ae.col(j) * scale * (bc - x0[jj]);
243 be -= Ae.col(j) * scale * bc;
247 for (Eigen::Index k = 0; k < dmap0.size(); ++k)
248 b[dmap0[k]] += be[k];
252 template <
typename T>
253 void _lift_bc_exterior_facets(
254 Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
const Form<T>& a,
255 const Eigen::Ref<
const Eigen::Matrix<T, Eigen::Dynamic, 1>>& bc_values1,
256 const std::vector<bool>& bc_markers1,
257 const Eigen::Ref<
const Eigen::Matrix<T, Eigen::Dynamic, 1>>& x0,
260 assert(a.rank() == 2);
263 std::shared_ptr<const mesh::Mesh> mesh = a.mesh();
266 mesh->topology_mutable().create_entity_permutations();
268 const int gdim = mesh->geometry().dim();
269 const int tdim = mesh->topology().dim();
272 mesh->topology_mutable().create_entities(tdim - 1);
273 mesh->topology_mutable().create_connectivity(tdim - 1, tdim);
275 mesh->topology_mutable().create_entity_permutations();
278 assert(a.function_space(0));
279 assert(a.function_space(1));
280 std::shared_ptr<const fem::DofMap> dofmap0 = a.function_space(0)->dofmap();
281 std::shared_ptr<const fem::DofMap> dofmap1 = a.function_space(1)->dofmap();
286 Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> coeffs
289 const std::function<void(T*,
const T*,
const T*,
const double*,
const int*,
290 const std::uint8_t*,
const std::uint32_t)>& fn
291 = a.integrals().get_tabulate_tensor(IntegralType::exterior_facet, 0);
294 const graph::AdjacencyList<std::int32_t>& x_dofmap
295 = mesh->geometry().dofmap();
297 const int num_dofs_g = x_dofmap.num_links(0);
298 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
299 = mesh->geometry().x();
302 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
303 coordinate_dofs(num_dofs_g, gdim);
304 Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> Ae;
305 Eigen::Matrix<T, Eigen::Dynamic, 1> be;
308 if (!a.all_constants_set())
309 throw std::runtime_error(
"Unset constant in Form");
310 const Eigen::Array<T, Eigen::Dynamic, 1> constant_values =
pack_constants(a);
313 const mesh::Topology& topology = mesh->topology();
314 auto connectivity = topology.connectivity(tdim - 1, tdim);
315 assert(connectivity);
316 auto c_to_f = topology.connectivity(tdim, tdim - 1);
318 auto map = topology.index_map(tdim - 1);
321 const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>& perms
322 = mesh->topology().get_facet_permutations();
323 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info
324 = mesh->topology().get_cell_permutation_info();
326 std::set<std::int32_t> fwd_shared_facets;
328 if (topology.index_map(tdim)->num_ghosts() == 0)
330 fwd_shared_facets.insert(
331 topology.index_map(tdim - 1)->shared_indices().begin(),
332 topology.index_map(tdim - 1)->shared_indices().end());
335 for (
int f = 0; f < map->size_local(); ++f)
340 if (connectivity->num_links(f) == 2
341 or fwd_shared_facets.find(f) != fwd_shared_facets.end())
345 assert(connectivity->num_links(f) == 1);
346 const std::int32_t cell = connectivity->links(f)[0];
349 auto facets = c_to_f->links(cell);
350 const auto* it = std::find(facets.data(), facets.data() + facets.rows(), f);
351 assert(it != (facets.data() + facets.rows()));
352 const int local_facet = std::distance(facets.data(), it);
354 const std::uint8_t perm = perms(local_facet, cell);
357 auto dmap1 = dofmap1->cell_dofs(cell);
361 for (Eigen::Index j = 0; j < dmap1.size(); ++j)
363 if (bc_markers1[dmap1[j]])
374 auto x_dofs = x_dofmap.links(cell);
375 for (
int i = 0; i < num_dofs_g; ++i)
376 for (
int j = 0; j < gdim; ++j)
377 coordinate_dofs(i, j) = x_g(x_dofs[i], j);
380 auto dmap0 = dofmap0->cell_dofs(cell);
385 auto coeff_array = coeffs.row(cell);
386 Ae.setZero(dmap0.size(), dmap1.size());
387 fn(Ae.data(), coeff_array.data(), constant_values.data(),
388 coordinate_dofs.data(), &local_facet, &perm, cell_info[cell]);
391 be.setZero(dmap0.size());
392 for (Eigen::Index j = 0; j < dmap1.size(); ++j)
394 const std::int32_t jj = dmap1[j];
397 const T bc = bc_values1[jj];
399 be -= Ae.col(j) * scale * (bc - x0[jj]);
401 be -= Ae.col(j) * scale * bc;
405 for (Eigen::Index k = 0; k < dmap0.size(); ++k)
406 b[dmap0[k]] += be[k];
411 template <
typename T>
412 void assemble_vector(Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
415 std::shared_ptr<const mesh::Mesh> mesh = L.mesh();
419 assert(L.function_space(0));
420 std::shared_ptr<const fem::DofMap> dofmap = L.function_space(0)->dofmap();
422 const graph::AdjacencyList<std::int32_t>& dofs = dofmap->list();
425 if (!L.all_constants_set())
426 throw std::runtime_error(
"Unset constant in Form");
427 const Eigen::Array<T, Eigen::Dynamic, 1> constant_values =
pack_constants(L);
430 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> coeffs
433 const FormIntegrals<T>& integrals = L.integrals();
434 for (
int i = 0; i < integrals.num_integrals(IntegralType::cell); ++i)
436 const auto& fn = integrals.get_tabulate_tensor(IntegralType::cell, i);
437 const std::vector<std::int32_t>& active_cells
438 = integrals.integral_domains(IntegralType::cell, i);
439 fem::impl::assemble_cells(b, *mesh, active_cells, dofs, fn, coeffs,
443 for (
int i = 0; i < integrals.num_integrals(IntegralType::exterior_facet);
447 = integrals.get_tabulate_tensor(IntegralType::exterior_facet, i);
448 const std::vector<std::int32_t>& active_facets
449 = integrals.integral_domains(IntegralType::exterior_facet, i);
450 fem::impl::assemble_exterior_facets(b, *mesh, active_facets, dofs, fn,
451 coeffs, constant_values);
454 const std::vector<int> c_offsets = L.coefficients().offsets();
455 for (
int i = 0; i < integrals.num_integrals(IntegralType::interior_facet);
459 = integrals.get_tabulate_tensor(IntegralType::interior_facet, i);
460 const std::vector<std::int32_t>& active_facets
461 = integrals.integral_domains(IntegralType::interior_facet, i);
462 fem::impl::assemble_interior_facets(b, *mesh, active_facets, *dofmap, fn,
463 coeffs, c_offsets, constant_values);
467 template <
typename T>
469 Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
const mesh::Mesh& mesh,
470 const std::vector<std::int32_t>& active_cells,
471 const graph::AdjacencyList<std::int32_t>& dofmap,
472 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
473 const std::uint8_t*,
const std::uint32_t)>& kernel,
474 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
476 const Eigen::Array<T, Eigen::Dynamic, 1>& constant_values)
478 const int gdim = mesh.geometry().dim();
479 mesh.topology_mutable().create_entity_permutations();
482 const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
485 const int num_dofs_g = x_dofmap.num_links(0);
486 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
487 = mesh.geometry().x();
491 const int num_dofs = dofmap.links(0).size();
492 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
493 coordinate_dofs(num_dofs_g, gdim);
494 Eigen::Matrix<T, Eigen::Dynamic, 1> be(num_dofs);
496 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info
497 = mesh.topology().get_cell_permutation_info();
500 for (std::int32_t c : active_cells)
503 auto x_dofs = x_dofmap.links(c);
504 for (
int i = 0; i < num_dofs_g; ++i)
505 for (
int j = 0; j < gdim; ++j)
506 coordinate_dofs(i, j) = x_g(x_dofs[i], j);
509 std::fill(be.data(), be.data() + num_dofs, 0);
510 kernel(be.data(), coeffs.row(c).data(), constant_values.data(),
511 coordinate_dofs.data(),
nullptr,
nullptr, cell_info[c]);
514 auto dofs = dofmap.links(c);
515 for (Eigen::Index i = 0; i < num_dofs; ++i)
520 template <
typename T>
521 void assemble_exterior_facets(
522 Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
const mesh::Mesh& mesh,
523 const std::vector<std::int32_t>& active_facets,
524 const graph::AdjacencyList<std::int32_t>& dofmap,
525 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
526 const std::uint8_t*,
const std::uint32_t)>& fn,
527 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
529 const Eigen::Array<T, Eigen::Dynamic, 1>& constant_values)
531 const int gdim = mesh.geometry().dim();
532 const int tdim = mesh.topology().dim();
535 mesh.topology_mutable().create_entities(tdim - 1);
536 mesh.topology_mutable().create_connectivity(tdim - 1, tdim);
537 mesh.topology_mutable().create_entity_permutations();
540 const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
543 const int num_dofs_g = x_dofmap.num_links(0);
544 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
545 = mesh.geometry().x();
549 const int num_dofs = dofmap.links(0).size();
550 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
551 coordinate_dofs(num_dofs_g, gdim);
552 Eigen::Matrix<T, Eigen::Dynamic, 1> be(num_dofs);
554 const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>& perms
555 = mesh.topology().get_facet_permutations();
556 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info
557 = mesh.topology().get_cell_permutation_info();
559 auto f_to_c = mesh.topology().connectivity(tdim - 1, tdim);
561 auto c_to_f = mesh.topology().connectivity(tdim, tdim - 1);
563 for (
const auto& f : active_facets)
566 assert(f_to_c->num_links(f) > 0);
567 const std::int32_t cell = f_to_c->links(f)[0];
570 auto facets = c_to_f->links(cell);
571 const auto* it = std::find(facets.data(), facets.data() + facets.rows(), f);
572 assert(it != (facets.data() + facets.rows()));
573 const int local_facet = std::distance(facets.data(), it);
574 const std::uint8_t perm = perms(local_facet, cell);
577 auto x_dofs = x_dofmap.links(cell);
578 for (
int i = 0; i < num_dofs_g; ++i)
579 for (
int j = 0; j < gdim; ++j)
580 coordinate_dofs(i, j) = x_g(x_dofs[i], j);
583 std::fill(be.data(), be.data() + num_dofs, 0);
584 fn(be.data(), coeffs.row(cell).data(), constant_values.data(),
585 coordinate_dofs.data(), &local_facet, &perm, cell_info[cell]);
588 auto dofs = dofmap.links(cell);
589 for (Eigen::Index i = 0; i < num_dofs; ++i)
594 template <
typename T>
595 void assemble_interior_facets(
596 Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
const mesh::Mesh& mesh,
597 const std::vector<std::int32_t>& active_facets,
const fem::DofMap& dofmap,
598 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
599 const std::uint8_t*,
const std::uint32_t)>& fn,
600 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
602 const std::vector<int>& offsets,
603 const Eigen::Array<T, Eigen::Dynamic, 1>& constant_values)
605 const int gdim = mesh.geometry().dim();
606 const int tdim = mesh.topology().dim();
609 mesh.topology_mutable().create_entities(tdim - 1);
610 mesh.topology_mutable().create_connectivity(tdim - 1, tdim);
611 mesh.topology_mutable().create_entity_permutations();
614 const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
617 const int num_dofs_g = x_dofmap.num_links(0);
618 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
619 = mesh.geometry().x();
622 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
623 coordinate_dofs(2 * num_dofs_g, gdim);
624 Eigen::Matrix<T, Eigen::Dynamic, 1> be;
625 Eigen::Array<T, Eigen::Dynamic, 1> coeff_array(2 * offsets.back());
626 assert(offsets.back() == coeffs.cols());
628 const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>& perms
629 = mesh.topology().get_facet_permutations();
630 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info
631 = mesh.topology().get_cell_permutation_info();
633 auto f_to_c = mesh.topology().connectivity(tdim - 1, tdim);
635 auto c_to_f = mesh.topology().connectivity(tdim, tdim - 1);
637 for (
const auto& f : active_facets)
640 auto cells = f_to_c->links(f);
641 assert(cells.rows() == 2);
644 std::array<int, 2> local_facet;
645 for (
int i = 0; i < 2; ++i)
647 auto facets = c_to_f->links(cells[i]);
649 = std::find(facets.data(), facets.data() + facets.rows(), f);
650 assert(it != (facets.data() + facets.rows()));
651 local_facet[i] = std::distance(facets.data(), it);
655 const std::array perm{perms(local_facet[0], cells[0]),
656 perms(local_facet[1], cells[1])};
659 auto x_dofs0 = x_dofmap.links(cells[0]);
660 auto x_dofs1 = x_dofmap.links(cells[1]);
661 for (
int i = 0; i < num_dofs_g; ++i)
663 for (
int j = 0; j < gdim; ++j)
665 coordinate_dofs(i, j) = x_g(x_dofs0[i], j);
666 coordinate_dofs(i + num_dofs_g, j) = x_g(x_dofs1[i], j);
671 auto dmap0 = dofmap.cell_dofs(cells[0]);
672 auto dmap1 = dofmap.cell_dofs(cells[1]);
676 auto coeff_cell0 = coeffs.row(cells[0]);
677 auto coeff_cell1 = coeffs.row(cells[1]);
680 for (std::size_t i = 0; i < offsets.size() - 1; ++i)
683 const int num_entries = offsets[i + 1] - offsets[i];
684 coeff_array.segment(2 * offsets[i], num_entries)
685 = coeff_cell0.segment(offsets[i], num_entries);
686 coeff_array.segment(offsets[i + 1] + offsets[i], num_entries)
687 = coeff_cell1.segment(offsets[i], num_entries);
691 be.setZero(dmap0.size() + dmap1.size());
692 fn(be.data(), coeff_array.data(), constant_values.data(),
693 coordinate_dofs.data(), local_facet.data(), perm.data(),
694 cell_info[cells[0]]);
697 for (Eigen::Index i = 0; i < dmap0.size(); ++i)
698 b[dmap0[i]] += be[i];
699 for (Eigen::Index i = 0; i < dmap1.size(); ++i)
700 b[dmap1[i]] += be[i + dmap0.size()];
704 template <
typename T>
706 Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
707 const std::vector<std::shared_ptr<
const Form<T>>> a,
708 const std::vector<std::vector<std::shared_ptr<
const DirichletBC<T>>>>& bcs1,
709 const std::vector<Eigen::Ref<
const Eigen::Matrix<T, Eigen::Dynamic, 1>>>&
714 if (!x0.empty() and x0.size() != a.size())
716 throw std::runtime_error(
717 "Mismatch in size between x0 and bilinear form in assembler.");
720 if (a.size() != bcs1.size())
722 throw std::runtime_error(
723 "Mismatch in size between a and bcs in assembler.");
726 for (std::size_t j = 0; j < a.size(); ++j)
728 std::vector<bool> bc_markers1;
729 Eigen::Matrix<T, Eigen::Dynamic, 1> bc_values1;
730 if (a[j] and !bcs1[j].empty())
732 auto V1 = a[j]->function_space(1);
734 auto map1 = V1->dofmap()->index_map;
737 = map1->block_size() * (map1->size_local() + map1->num_ghosts());
738 bc_markers1.assign(crange,
false);
739 bc_values1 = Eigen::Matrix<T, Eigen::Dynamic, 1>::Zero(crange);
740 for (
const std::shared_ptr<
const DirichletBC<T>>& bc : bcs1[j])
742 bc->mark_dofs(bc_markers1);
743 bc->dof_values(bc_values1);
748 lift_bc<T>(b, *a[j], bc_values1, bc_markers1, x0[j], scale);
750 lift_bc<T>(b, *a[j], bc_values1, bc_markers1, scale);
755 template <
typename T>
757 Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
const Form<T>& a,
758 const Eigen::Ref<
const Eigen::Matrix<T, Eigen::Dynamic, 1>>& bc_values1,
759 const std::vector<bool>& bc_markers1,
double scale)
761 const Eigen::Matrix<T, Eigen::Dynamic, 1> x0(0);
762 if (a.integrals().num_integrals(fem::IntegralType::cell) > 0)
763 _lift_bc_cells<T>(b, a, bc_values1, bc_markers1, x0, scale);
764 if (a.integrals().num_integrals(fem::IntegralType::exterior_facet) > 0)
765 _lift_bc_exterior_facets<T>(b, a, bc_values1, bc_markers1, x0, scale);
768 template <
typename T>
770 Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
const Form<T>& a,
771 const Eigen::Ref<
const Eigen::Matrix<T, Eigen::Dynamic, 1>>& bc_values1,
772 const std::vector<bool>& bc_markers1,
773 const Eigen::Ref<
const Eigen::Matrix<T, Eigen::Dynamic, 1>>& x0,
776 if (a.integrals().num_integrals(fem::IntegralType::cell) > 0)
777 _lift_bc_cells(b, a, bc_values1, bc_markers1, x0, scale);
778 if (a.integrals().num_integrals(fem::IntegralType::exterior_facet) > 0)
779 _lift_bc_exterior_facets(b, a, bc_values1, bc_markers1, x0, scale);