40#ifdef DEAL_II_WITH_MPI
60 check_primary_dof_list(
62 const std::vector<types::global_dof_index> &primary_dof_list)
64 const unsigned int N = primary_dof_list.size();
67 for (
unsigned int i = 0; i <
N; ++i)
68 for (
unsigned int j = 0; j <
N; ++j)
69 tmp(i, j) = face_interpolation_matrix(primary_dof_list[i], j);
80 double diagonal_sum = 0;
81 for (
unsigned int i = 0; i <
N; ++i)
82 diagonal_sum += std::fabs(tmp(i, i));
83 const double typical_diagonal_element = diagonal_sum /
N;
87 std::vector<unsigned int> p(N);
88 for (
unsigned int i = 0; i <
N; ++i)
91 for (
unsigned int j = 0; j <
N; ++j)
95 double max = std::fabs(tmp(j, j));
97 for (
unsigned int i = j + 1; i <
N; ++i)
99 if (std::fabs(tmp(i, j)) > max)
101 max = std::fabs(tmp(i, j));
108 if (max < 1.e-12 * typical_diagonal_element)
114 for (
unsigned int k = 0; k <
N; ++k)
115 std::swap(tmp(j, k), tmp(r, k));
117 std::swap(p[j], p[r]);
121 const double hr = 1. / tmp(j, j);
123 for (
unsigned int k = 0; k <
N; ++k)
127 for (
unsigned int i = 0; i <
N; ++i)
131 tmp(i, k) -= tmp(i, j) * tmp(j, k) * hr;
134 for (
unsigned int i = 0; i <
N; ++i)
170 template <
int dim,
int spacedim>
172 select_primary_dofs_for_face_restriction(
176 std::vector<bool> &primary_dof_mask)
182 const unsigned int face_no = 0;
217 std::vector<types::global_dof_index> primary_dof_list;
218 unsigned int index = 0;
223 unsigned int dofs_added = 0;
232 primary_dof_list.push_back(index + i);
235 if (check_primary_dof_list(face_interpolation_matrix,
236 primary_dof_list) ==
true)
241 primary_dof_list.pop_back();
254 unsigned int dofs_added = 0;
260 primary_dof_list.push_back(index + i);
261 if (check_primary_dof_list(face_interpolation_matrix,
262 primary_dof_list) ==
true)
265 primary_dof_list.pop_back();
277 unsigned int dofs_added = 0;
283 primary_dof_list.push_back(index + i);
284 if (check_primary_dof_list(face_interpolation_matrix,
285 primary_dof_list) ==
true)
288 primary_dof_list.pop_back();
299 std::fill(primary_dof_mask.begin(), primary_dof_mask.end(),
false);
300 for (
const auto dof : primary_dof_list)
301 primary_dof_mask[dof] = true;
310 template <
int dim,
int spacedim>
312 ensure_existence_of_primary_dof_mask(
316 std::unique_ptr<std::vector<bool>> &primary_dof_mask)
322 const unsigned int face_no = 0;
324 if (primary_dof_mask ==
nullptr)
328 select_primary_dofs_for_face_restriction(fe1,
330 face_interpolation_matrix,
342 template <
int dim,
int spacedim>
344 ensure_existence_of_face_matrix(
353 const unsigned int face_no = 0;
355 if (matrix ==
nullptr)
357 matrix = std::make_unique<FullMatrix<double>>(
368 template <
int dim,
int spacedim>
370 ensure_existence_of_subface_matrix(
373 const unsigned int subface,
380 const unsigned int face_no = 0;
382 if (matrix ==
nullptr)
384 matrix = std::make_unique<FullMatrix<double>>(
401 ensure_existence_of_split_face_matrix(
403 const std::vector<bool> &primary_dof_mask,
408 Assert(std::count(primary_dof_mask.begin(),
409 primary_dof_mask.end(),
411 static_cast<signed int>(face_interpolation_matrix.
n()),
414 if (split_matrix ==
nullptr)
416 split_matrix = std::make_unique<
419 const unsigned int n_primary_dofs = face_interpolation_matrix.
n();
420 const unsigned int n_dofs = face_interpolation_matrix.
m();
426 split_matrix->first.reinit(n_primary_dofs, n_primary_dofs);
427 split_matrix->second.reinit(n_dofs - n_primary_dofs,
430 unsigned int nth_primary_dof = 0, nth_dependent_dof = 0;
432 for (
unsigned int i = 0; i < n_dofs; ++i)
433 if (primary_dof_mask[i] ==
true)
435 for (
unsigned int j = 0; j < n_primary_dofs; ++j)
436 split_matrix->first(nth_primary_dof, j) =
437 face_interpolation_matrix(i, j);
442 for (
unsigned int j = 0; j < n_primary_dofs; ++j)
443 split_matrix->second(nth_dependent_dof, j) =
444 face_interpolation_matrix(i, j);
453 split_matrix->first.gauss_jordan();
463 template <
int dim,
int spacedim>
485 template <
typename number1,
typename number2>
488 const std::vector<types::global_dof_index> &primary_dofs,
489 const std::vector<types::global_dof_index> &dependent_dofs,
493 Assert(face_constraints.
n() == primary_dofs.size(),
495 Assert(face_constraints.
m() == dependent_dofs.size(),
497 face_constraints.
m()));
499 const unsigned int n_primary_dofs = primary_dofs.size();
500 const unsigned int n_dependent_dofs = dependent_dofs.size();
504 for (
unsigned int row = 0; row != n_dependent_dofs; ++row)
507 for (
unsigned int col = 0; col != n_primary_dofs; ++col)
518 boost::container::small_vector<std::pair<size_type, size_type>, 25>
520 sorted_primary_dofs.reserve(n_primary_dofs);
521 for (
unsigned int i = 0; i < n_primary_dofs; ++i)
522 sorted_primary_dofs.emplace_back(primary_dofs[i], i);
523 std::sort(sorted_primary_dofs.begin(), sorted_primary_dofs.end());
525 boost::container::small_vector<std::pair<size_type, number2>, 25>
527 entries.reserve(n_primary_dofs);
528 for (
unsigned int row = 0; row != n_dependent_dofs; ++row)
548 bool is_trivial_constraint =
false;
550 for (
unsigned int i = 0; i < n_primary_dofs; ++i)
551 if (face_constraints(row, i) == 1.0)
552 if (dependent_dofs[row] == primary_dofs[i])
554 is_trivial_constraint =
true;
556 for (
unsigned int ii = 0; ii < n_primary_dofs; ++ii)
558 Assert(face_constraints(row, ii) == 0.0,
564 if (is_trivial_constraint ==
true)
579 for (
const auto &[dof_index, unsorted_index] :
581 if (
std::
fabs(face_constraints(row, unsorted_index)) >= 1
e-14)
582 entries.emplace_back(dof_index,
583 face_constraints(row, unsorted_index));
594 template <
typename number,
int spacedim>
605 template <
typename number,
int spacedim>
608 const ::DoFHandler<1, spacedim> & ,
610 std::integral_constant<int, 1>)
617 template <
typename number,
int spacedim>
622 std::integral_constant<int, 1>)
629 template <
int dim_,
int spacedim,
typename number>
634 std::integral_constant<int, 2>)
636 const unsigned int dim = 2;
638 std::vector<types::global_dof_index> dofs_on_mother;
639 std::vector<types::global_dof_index> dofs_on_children;
646 boost::container::small_vector<
647 std::pair<typename AffineConstraints<number>::size_type, number>,
663 if (cell->is_artificial())
666 for (
const unsigned int face : cell->face_indices())
667 if (cell->face(face)->has_children())
673 Assert(cell->face(face)->n_active_fe_indices() == 1,
675 Assert(cell->face(face)->fe_index_is_active(
676 cell->active_fe_index()) ==
true,
678 for (
unsigned int c = 0; c < cell->face(face)->n_children();
680 if (!cell->neighbor_child_on_subface(face, c)
682 Assert(cell->face(face)->child(c)->n_active_fe_indices() ==
688 for (
unsigned int c = 0; c < cell->face(face)->n_children();
690 if (!cell->neighbor_child_on_subface(face, c)
692 Assert(cell->face(face)->child(c)->fe_index_is_active(
693 cell->active_fe_index()) ==
true,
700 const unsigned int n_dofs_on_mother =
707 dofs_on_mother.resize(n_dofs_on_mother);
710 dofs_on_children.clear();
711 dofs_on_children.reserve(n_dofs_on_children);
721 this_face = cell->face(face);
725 unsigned int next_index = 0;
726 for (
unsigned int vertex = 0; vertex < 2; ++vertex)
729 dofs_on_mother[next_index++] =
730 this_face->vertex_dof_index(vertex, dof, fe_index);
732 dofs_on_mother[next_index++] =
733 this_face->dof_index(dof, fe_index);
737 dofs_on_children.push_back(
738 this_face->child(0)->vertex_dof_index(1, dof, fe_index));
739 for (
unsigned int child = 0; child < 2; ++child)
742 if (cell->neighbor_child_on_subface(face, child)
747 dofs_on_children.push_back(
748 this_face->child(child)->dof_index(dof, fe_index));
751 Assert(dofs_on_children.size() <= n_dofs_on_children,
755 for (
unsigned int row = 0; row != dofs_on_children.size();
758 constraint_entries.clear();
759 constraint_entries.reserve(dofs_on_mother.size());
760 for (
unsigned int i = 0; i != dofs_on_mother.size(); ++i)
761 constraint_entries.emplace_back(dofs_on_mother[i],
775 if (!cell->at_boundary(face) &&
776 !cell->neighbor(face)->is_artificial())
778 Assert(cell->face(face)->n_active_fe_indices() == 1,
780 Assert(cell->face(face)->fe_index_is_active(
781 cell->active_fe_index()) ==
true,
790 template <
int dim_,
int spacedim,
typename number>
795 std::integral_constant<int, 3>)
797 const unsigned int dim = 3;
799 std::vector<types::global_dof_index> dofs_on_mother;
800 std::vector<types::global_dof_index> dofs_on_children;
807 boost::container::small_vector<
808 std::pair<typename AffineConstraints<number>::size_type, number>,
824 if (cell->is_artificial())
827 for (
const unsigned int face : cell->face_indices())
828 if (cell->face(face)->has_children())
833 if (cell->get_fe().n_dofs_per_face(face) == 0)
836 Assert(cell->face(face)->refinement_case() ==
845 Assert(cell->face(face)->fe_index_is_active(
846 cell->active_fe_index()) ==
true,
848 for (
unsigned int c = 0; c < cell->face(face)->n_children();
850 if (!cell->neighbor_child_on_subface(face, c)
853 cell->face(face)->child(c)->n_active_fe_indices(), 1);
858 for (
unsigned int c = 0; c < cell->face(face)->n_children();
860 if (!cell->neighbor_child_on_subface(face, c)
863 Assert(cell->face(face)->child(c)->fe_index_is_active(
864 cell->active_fe_index()) ==
true,
866 for (
unsigned int e = 0; e < 4; ++e)
871 ->n_active_fe_indices() == 1,
876 ->fe_index_is_active(
877 cell->active_fe_index()) ==
true,
881 for (
unsigned int e = 0; e < 4; ++e)
883 Assert(cell->face(face)->line(e)->n_active_fe_indices() ==
886 Assert(cell->face(face)->line(e)->fe_index_is_active(
887 cell->active_fe_index()) ==
true,
896 const unsigned int n_dofs_on_children =
903 dofs_on_mother.resize(n_dofs_on_mother);
906 dofs_on_children.clear();
907 dofs_on_children.reserve(n_dofs_on_children);
917 this_face = cell->face(face);
921 unsigned int next_index = 0;
922 for (
unsigned int vertex = 0; vertex < 4; ++vertex)
925 dofs_on_mother[next_index++] =
926 this_face->vertex_dof_index(vertex, dof, fe_index);
927 for (
unsigned int line = 0; line < 4; ++line)
929 dofs_on_mother[next_index++] =
930 this_face->line(line)->dof_index(dof, fe_index);
933 dofs_on_mother[next_index++] =
934 this_face->dof_index(dof, fe_index);
943 ((this_face->child(0)->vertex_index(3) ==
944 this_face->child(1)->vertex_index(2)) &&
945 (this_face->child(0)->vertex_index(3) ==
946 this_face->child(2)->vertex_index(1)) &&
947 (this_face->child(0)->vertex_index(3) ==
948 this_face->child(3)->vertex_index(0))),
952 dofs_on_children.push_back(
953 this_face->child(0)->vertex_dof_index(3, dof));
956 for (
unsigned int line = 0; line < 4; ++line)
959 dofs_on_children.push_back(
960 this_face->line(line)->child(0)->vertex_dof_index(
967 dofs_on_children.push_back(
968 this_face->child(0)->line(1)->dof_index(dof, fe_index));
970 dofs_on_children.push_back(
971 this_face->child(2)->line(1)->dof_index(dof, fe_index));
973 dofs_on_children.push_back(
974 this_face->child(0)->line(3)->dof_index(dof, fe_index));
976 dofs_on_children.push_back(
977 this_face->child(1)->line(3)->dof_index(dof, fe_index));
980 for (
unsigned int line = 0; line < 4; ++line)
981 for (
unsigned int child = 0; child < 2; ++child)
985 dofs_on_children.push_back(
986 this_face->line(line)->child(child)->dof_index(
991 for (
unsigned int child = 0; child < 4; ++child)
994 if (cell->neighbor_child_on_subface(face, child)
999 dofs_on_children.push_back(
1000 this_face->child(child)->dof_index(dof, fe_index));
1004 Assert(dofs_on_children.size() <= n_dofs_on_children,
1012 for (
unsigned int row = 0; row != dofs_on_children.size();
1017 constraint_entries.clear();
1018 constraint_entries.reserve(dofs_on_mother.size());
1019 for (
unsigned int i = 0; i != dofs_on_mother.size(); ++i)
1020 constraint_entries.emplace_back(dofs_on_mother[i],
1035 if (!cell->at_boundary(face) &&
1036 !cell->neighbor(face)->is_artificial())
1038 Assert(cell->face(face)->n_active_fe_indices() == 1,
1040 Assert(cell->face(face)->fe_index_is_active(
1041 cell->active_fe_index()) ==
true,
1050 template <
int dim_,
int spacedim,
typename number>
1055 std::integral_constant<int, 2>)
1062 const unsigned int dim = 2;
1064 std::vector<types::global_dof_index> face_dof_indices;
1065 std::map<types::global_dof_index, std::set<types::global_dof_index>>
1072 if (cell->is_artificial())
1076 for (
const unsigned int f : cell->face_indices())
1080 if (!cell->face(f)->has_children())
1083 Assert(cell->face(f)->n_active_fe_indices() == 1,
1085 Assert(cell->face(f)->fe_index_is_active(
1086 cell->active_fe_index()) ==
true,
1091 for (
unsigned int c = 0; c < cell->face(f)->n_children(); ++c)
1093 if (cell->neighbor_child_on_subface(f, c)
1097 Assert(cell->face(f)->child(c)->n_active_fe_indices() ==
1101 Assert(cell->face(f)->child(c)->fe_index_is_active(
1102 cell->active_fe_index()) ==
true,
1111 face_dof_indices.resize(n_dofs);
1113 cell->face(f)->get_dof_indices(face_dof_indices);
1114 const std::vector<types::global_dof_index> dof_on_mother_face =
1117 cell->face(f)->child(0)->get_dof_indices(face_dof_indices);
1118 const std::vector<types::global_dof_index> dof_on_child_face_0 =
1121 cell->face(f)->child(1)->get_dof_indices(face_dof_indices);
1122 const std::vector<types::global_dof_index> dof_on_child_face_1 =
1131 const bool direction_mother = (cell->face(f)->vertex_index(0) >
1132 cell->face(f)->vertex_index(1)) ?
1135 const bool direction_child_0 =
1136 (cell->face(f)->child(0)->vertex_index(0) >
1137 cell->face(f)->child(0)->vertex_index(1)) ?
1140 const bool direction_child_1 =
1141 (cell->face(f)->child(1)->vertex_index(0) >
1142 cell->face(f)->child(1)->vertex_index(1)) ?
1146 for (
unsigned int row = 0; row < n_dofs; ++row)
1148 constraints.
add_line(dof_on_child_face_0[row]);
1149 constraints.
add_line(dof_on_child_face_1[row]);
1152 for (
unsigned int row = 0; row < n_dofs; ++row)
1154 for (
unsigned int dof_i_on_mother = 0;
1155 dof_i_on_mother < n_dofs;
1162 unsigned int shift_0 =
1163 (direction_mother == direction_child_0) ? 0 : n_dofs;
1164 constraints.
add_entry(dof_on_child_face_0[row],
1165 dof_on_mother_face[dof_i_on_mother],
1169 unsigned int shift_1 =
1170 (direction_mother == direction_child_1) ? 0 : n_dofs;
1171 constraints.
add_entry(dof_on_child_face_1[row],
1172 dof_on_mother_face[dof_i_on_mother],
1182 template <
int dim_,
int spacedim,
typename number>
1187 std::integral_constant<int, 3>)
1194 const unsigned int dim = 3;
1196 std::vector<types::global_dof_index> dofs_on_mother;
1197 std::vector<types::global_dof_index> dofs_on_children;
1203 if (cell->is_artificial())
1207 for (
const unsigned int face : cell->face_indices())
1210 if (cell->face(face)->has_children() ==
false)
1213 if (cell->get_fe().n_dofs_per_face(face) == 0)
1216 Assert(cell->face(face)->refinement_case() ==
1222 Assert(cell->face(face)->fe_index_is_active(
1223 cell->active_fe_index()) ==
true,
1228 for (
unsigned int c = 0; c < cell->face(face)->n_children();
1231 if (cell->neighbor_child_on_subface(face, c)
1236 cell->face(face)->child(c)->n_active_fe_indices(), 1);
1238 Assert(cell->face(face)->child(c)->fe_index_is_active(
1239 cell->active_fe_index()) ==
true,
1242 for (
unsigned int e = 0;
1243 e < GeometryInfo<dim>::vertices_per_face;
1249 ->n_active_fe_indices() == 1,
1255 ->fe_index_is_active(
1256 cell->active_fe_index()) ==
true,
1261 for (
unsigned int e = 0;
1262 e < GeometryInfo<dim>::vertices_per_face;
1265 Assert(cell->face(face)->line(e)->n_active_fe_indices() ==
1269 Assert(cell->face(face)->line(e)->fe_index_is_active(
1270 cell->active_fe_index()) ==
true,
1277 const unsigned int fe_index = cell->active_fe_index();
1280 unsigned int degree(fe.
degree);
1285 dofs_on_mother.resize(n_dofs_on_mother);
1287 const unsigned int n_lines_on_mother =
1302 const unsigned int n_internal_lines_on_children = 4;
1314 const unsigned int n_external_lines_on_children = 8;
1316 const unsigned int n_lines_on_children =
1317 n_internal_lines_on_children + n_external_lines_on_children;
1320 const unsigned int n_children_per_face =
1322 const unsigned int n_children_per_line =
1328 const unsigned int n_dofs_on_children =
1332 dofs_on_children.clear();
1333 dofs_on_children.reserve(n_dofs_on_children);
1343 unsigned int next_index = 0;
1349 for (
unsigned int line = 0;
1350 line < GeometryInfo<dim>::lines_per_face;
1353 dofs_on_mother[next_index++] =
1354 this_face->line(line)->dof_index(dof, fe_index);
1358 dofs_on_mother[next_index++] =
1359 this_face->dof_index(dof, fe_index);
1378 dofs_on_children.push_back(
1379 this_face->child(0)->line(1)->dof_index(dof, fe_index));
1382 dofs_on_children.push_back(
1383 this_face->child(2)->line(1)->dof_index(dof, fe_index));
1386 dofs_on_children.push_back(
1387 this_face->child(0)->line(3)->dof_index(dof, fe_index));
1390 dofs_on_children.push_back(
1391 this_face->child(1)->line(3)->dof_index(dof, fe_index));
1396 for (
unsigned int line = 0;
1397 line < GeometryInfo<dim>::lines_per_face;
1399 for (
unsigned int child = 0; child < n_children_per_line;
1402 dofs_on_children.push_back(
1403 this_face->line(line)->child(child)->dof_index(dof,
1407 for (
unsigned int child = 0; child < n_children_per_face; ++child)
1410 if (cell->neighbor_child_on_subface(face, child)
1416 dofs_on_children.push_back(
1417 this_face->child(child)->dof_index(dof, fe_index));
1422 Assert(dofs_on_children.size() <= n_dofs_on_children,
1432 std::vector<bool> direction_mother(
1434 for (
unsigned int line = 0;
1435 line < GeometryInfo<dim>::lines_per_face;
1437 if (this_face->line(line)->vertex_index(0) >
1438 this_face->line(line)->vertex_index(1))
1439 direction_mother[line] =
true;
1442 std::vector<bool> direction_child_intern(
1443 n_internal_lines_on_children,
false);
1448 unsigned int center = this_face->child(0)->vertex_index(3);
1451 for (
unsigned int line = 0; line < n_internal_lines_on_children;
1455 direction_child_intern[line] =
1456 this_face->line(line)->child(0)->vertex_index(1) <
1463 direction_child_intern[line] =
1464 this_face->line(line)->child(0)->vertex_index(1) >
1471 std::vector<bool> direction_child(n_external_lines_on_children,
1473 for (
unsigned int line = 0;
1474 line < GeometryInfo<dim>::lines_per_face;
1477 if (this_face->line(line)->child(0)->vertex_index(0) >
1478 this_face->line(line)->child(0)->vertex_index(1))
1479 direction_child[2 * line] =
true;
1480 if (this_face->line(line)->child(1)->vertex_index(0) >
1481 this_face->line(line)->child(1)->vertex_index(1))
1482 direction_child[2 * line + 1] =
true;
1487 bool mother_flip_x =
false;
1488 bool mother_flip_y =
false;
1489 bool mother_flip_xy =
false;
1490 std::vector<bool> child_flip_x(n_children_per_face,
false);
1491 std::vector<bool> child_flip_y(n_children_per_face,
false);
1492 std::vector<bool> child_flip_xy(n_children_per_face,
false);
1495 [2] = {{1, 2}, {0, 3}, {3, 0}, {2, 1}};
1500 unsigned int current_glob = cell->face(face)->vertex_index(0);
1501 unsigned int current_max = 0;
1502 for (
unsigned int v = 1;
1503 v < GeometryInfo<dim>::vertices_per_face;
1505 if (current_glob < this_face->vertex_index(v))
1508 current_glob = this_face->vertex_index(v);
1513 if (current_max < 2)
1514 mother_flip_y =
true;
1518 if (current_max % 2 == 0)
1519 mother_flip_x =
true;
1522 if (this_face->vertex_index(
1523 vertices_adjacent_on_face[current_max][0]) <
1524 this_face->vertex_index(
1525 vertices_adjacent_on_face[current_max][1]))
1526 mother_flip_xy =
true;
1531 for (
unsigned int child = 0; child < n_children_per_face; ++child)
1533 unsigned int current_max = 0;
1534 unsigned int current_glob =
1535 this_face->child(child)->vertex_index(0);
1537 for (
unsigned int v = 1;
1538 v < GeometryInfo<dim>::vertices_per_face;
1540 if (current_glob < this_face->child(child)->vertex_index(v))
1543 current_glob = this_face->child(child)->vertex_index(v);
1546 if (current_max < 2)
1547 child_flip_y[child] =
true;
1549 if (current_max % 2 == 0)
1550 child_flip_x[child] =
true;
1552 if (this_face->child(child)->vertex_index(
1553 vertices_adjacent_on_face[current_max][0]) <
1554 this_face->child(child)->vertex_index(
1555 vertices_adjacent_on_face[current_max][1]))
1556 child_flip_xy[child] =
true;
1558 child_flip_xy[child] = mother_flip_xy;
1562 std::vector<std::vector<double>> constraints_matrix(
1565 std::vector<double>(dofs_on_mother.size(), 0));
1570 for (
unsigned int line = 0; line < n_internal_lines_on_children;
1574 unsigned int line_mother = line / 2;
1575 unsigned int row_mother =
1579 for (
unsigned int i = 0;
1582 constraints_matrix[row + row_start][i] =
1586 for (
unsigned int line = 0; line < n_internal_lines_on_children;
1590 unsigned int line_mother = line / 2;
1591 unsigned int row_mother =
1595 for (
unsigned int i =
1597 i < dofs_on_mother.size();
1599 constraints_matrix[row + row_start][i] =
1604 unsigned int row_offset =
1606 for (
unsigned int line = 0; line < n_external_lines_on_children;
1610 unsigned int line_mother = line / 2;
1611 unsigned int row_mother =
1613 for (
unsigned int row = row_offset;
1616 for (
unsigned int i = 0; i < dofs_on_mother.size(); ++i)
1617 constraints_matrix[row + row_start][i] =
1623 for (
unsigned int face = 0; face < n_children_per_face; ++face)
1626 for (
unsigned int row = row_offset;
1629 for (
unsigned int i = 0; i < dofs_on_mother.size(); ++i)
1630 constraints_matrix[row + row_start][i] =
1639 for (
unsigned int i = 0;
1644 unsigned int tmp_i = i % degree;
1647 for (
unsigned int j = 0;
1652 unsigned int tmp_j = j % degree;
1654 if ((line_i < 2 && line_j < 2) ||
1655 (line_i >= 2 && line_j >= 2))
1657 if (direction_child_intern[line_i] !=
1658 direction_mother[line_j])
1660 if ((tmp_i + tmp_j) % 2 == 1)
1662 constraints_matrix[i][j] *= -1.0;
1668 if (direction_mother[line_i])
1670 if ((tmp_i + tmp_j) % 2 == 1)
1672 constraints_matrix[i][j] *= -1.0;
1680 for (
unsigned int i =
1686 unsigned int tmp_i = i % degree;
1689 for (
unsigned int j = 0;
1694 unsigned int tmp_j = j % degree;
1696 if (direction_child[line_i] != direction_mother[line_j])
1698 if ((tmp_i + tmp_j) % 2 == 1)
1700 constraints_matrix[i][j] *= -1.0;
1718 unsigned int tmp_i = i % degree;
1720 unsigned int start_j =
1723 for (
unsigned int block = 0; block < n_blocks; ++block)
1726 for (
unsigned int jy = 0; jy < degree - 1; ++jy)
1727 for (
unsigned int jx = 0; jx < degree - 1; ++jx)
1729 unsigned int j = start_j + jx + (jy * (degree - 1));
1730 if (direction_child_intern[line_i] != mother_flip_y)
1732 if ((jy + tmp_i) % 2 == 0)
1734 constraints_matrix[i][j] *= -1.0;
1739 start_j += (degree - 1) * (degree - 1);
1742 for (
unsigned int jy = 0; jy < degree - 1; ++jy)
1743 for (
unsigned int jx = 0; jx < degree - 1; ++jx)
1745 unsigned int j = start_j + jx + (jy * (degree - 1));
1747 if (direction_child_intern[line_i] != mother_flip_y)
1749 if ((jy + tmp_i) % 2 == 0)
1751 constraints_matrix[i][j] *= -1.0;
1755 start_j += (degree - 1) * (degree - 1);
1759 start_j += degree - 1;
1763 start_j += degree - 1;
1773 unsigned int tmp_i = i % degree;
1775 unsigned int start_j =
1778 for (
unsigned int block = 0; block < n_blocks; block++)
1781 for (
unsigned int jy = 0; jy < degree - 1; ++jy)
1782 for (
unsigned int jx = 0; jx < degree - 1; ++jx)
1784 unsigned int j = start_j + jx + (jy * (degree - 1));
1785 if (direction_child_intern[line_i] != mother_flip_x)
1787 if ((jx + tmp_i) % 2 == 0)
1789 constraints_matrix[i][j] *= -1.0;
1794 start_j += (degree - 1) * (degree - 1);
1797 for (
unsigned int jy = 0; jy < degree - 1; ++jy)
1798 for (
unsigned int jx = 0; jx < degree - 1; ++jx)
1800 unsigned int j = start_j + jx + (jy * (degree - 1));
1801 if (direction_child_intern[line_i] != mother_flip_x)
1803 if ((jx + tmp_i) % 2 == 0)
1805 constraints_matrix[i][j] *= -1.0;
1809 start_j += (degree - 1) * (degree - 1);
1813 start_j += degree - 1;
1817 start_j += degree - 1;
1822 unsigned int degree_square = (degree - 1) * (degree - 1);
1826 for (
unsigned int child_face = 0;
1827 child_face < n_children_per_face;
1829 for (
unsigned int block = 0; block < n_blocks; ++block)
1841 for (
unsigned int iy = 0; iy < degree - 1; ++iy)
1842 for (
unsigned int ix = 0; ix < degree - 1; ++ix)
1848 for (
unsigned int jy = 0; jy < degree - 1; ++jy)
1849 for (
unsigned int jx = 0; jx < degree - 1; ++jx)
1851 if (child_flip_x[child_face] !=
1854 if ((ix + jx) % 2 == 1)
1856 constraints_matrix[i][j] *= -1.0;
1860 if (child_flip_y[child_face] !=
1863 if ((iy + jy) % 2 == 1)
1865 constraints_matrix[i][j] *= -1.0;
1875 for (
unsigned int iy = 0; iy < degree - 1; ++iy)
1876 for (
unsigned int ix = 0; ix < degree - 1; ++ix)
1881 degree_square + block * block_size;
1882 for (
unsigned int jy = 0; jy < degree - 1; ++jy)
1883 for (
unsigned int jx = 0; jx < degree - 1; ++jx)
1885 if (child_flip_x[child_face] !=
1888 if ((ix + jx) % 2 == 1)
1890 constraints_matrix[i][j] *= -1.0;
1894 if (child_flip_y[child_face] !=
1897 if ((iy + jy) % 2 == 1)
1899 constraints_matrix[i][j] *= -1.0;
1911 for (
unsigned int iy = 0; iy < degree - 1; ++iy)
1916 degree_square + block * block_size;
1917 for (
unsigned int jy = 0; jy < degree - 1; ++jy)
1918 for (
unsigned int jx = 0; jx < degree - 1; ++jx)
1920 if (child_flip_x[child_face] !=
1925 constraints_matrix[i][j] *= -1.0;
1929 if (child_flip_y[child_face] !=
1932 if ((iy + jy) % 2 == 1)
1934 constraints_matrix[i][j] *= -1.0;
1943 2 * degree_square + block * block_size;
1944 for (
unsigned int jy = 0; jy < degree - 1; ++jy)
1946 if (child_flip_y[child_face] !=
1949 if ((iy + jy) % 2 == 1)
1951 constraints_matrix[i][j] *= -1.0;
1961 for (
unsigned int ix = 0; ix < degree - 1; ++ix)
1966 degree_square + block * block_size;
1967 for (
unsigned int jy = 0; jy < degree - 1; ++jy)
1968 for (
unsigned int jx = 0; jx < degree - 1; ++jx)
1970 if (child_flip_x[child_face] !=
1973 if ((ix + jx) % 2 == 1)
1975 constraints_matrix[i][j] *= -1.0;
1979 if (child_flip_y[child_face] !=
1984 constraints_matrix[i][j] *= -1.0;
1993 2 * degree_square + (degree - 1) +
1995 for (
unsigned int jx = 0; jx < degree - 1; ++jx)
1997 if (child_flip_x[child_face] !=
2000 if ((ix + jx) % 2 == 1)
2002 constraints_matrix[i][j] *= -1.0;
2019 for (
unsigned int i = 0;
2027 std::vector<double> constraints_matrix_old(
2028 dofs_on_mother.size(), 0);
2029 for (
unsigned int j = 0; j < dofs_on_mother.size(); ++j)
2031 constraints_matrix_old[j] = constraints_matrix[i][j];
2034 unsigned int j_start =
2036 for (
unsigned block = 0; block < n_blocks; block++)
2039 for (
unsigned int jy = 0; jy < degree - 1; ++jy)
2040 for (
unsigned int jx = 0; jx < degree - 1; ++jx)
2042 unsigned int j_old =
2043 j_start + jx + (jy * (degree - 1));
2044 unsigned int j_new =
2045 j_start + jy + (jx * (degree - 1));
2046 constraints_matrix[i][j_new] =
2047 constraints_matrix_old[j_old];
2049 j_start += degree_square;
2052 for (
unsigned int jy = 0; jy < degree - 1; ++jy)
2053 for (
unsigned int jx = 0; jx < degree - 1; ++jx)
2055 unsigned int j_old =
2056 j_start + jx + (jy * (degree - 1));
2057 unsigned int j_new =
2058 j_start + jy + (jx * (degree - 1));
2059 constraints_matrix[i][j_new] =
2060 -constraints_matrix_old[j_old];
2062 j_start += degree_square;
2065 for (
unsigned int j = j_start;
2066 j < j_start + (degree - 1);
2069 constraints_matrix[i][j] =
2070 constraints_matrix_old[j + (degree - 1)];
2071 constraints_matrix[i][j + (degree - 1)] =
2072 constraints_matrix_old[j];
2074 j_start += 2 * (degree - 1);
2081 const unsigned int deg = degree - 1;
2084 std::vector<std::vector<double>> constraints_matrix_old(
2087 for (
unsigned int i = 0;
2091 constraints_matrix_old[i][j] = constraints_matrix
2096 for (
unsigned int child = 0; child < n_children_per_face;
2099 if (!child_flip_xy[child])
2102 unsigned int i_start_new =
2107 unsigned int j_start =
2110 for (
unsigned int block = 0; block < n_blocks; block++)
2113 for (
unsigned int ix = 0; ix < deg; ++ix)
2115 for (
unsigned int iy = 0; iy < deg; ++iy)
2117 for (
unsigned int j = 0;
2120 constraints_matrix[i_start_new + iy +
2121 (ix * deg)][j + j_start] =
2122 constraints_matrix_old[i_start_old + ix +
2126 i_start_new += deg * deg;
2127 i_start_old += deg * deg;
2130 for (
unsigned int ix = 0; ix < deg; ++ix)
2132 for (
unsigned int iy = 0; iy < deg; ++iy)
2134 for (
unsigned int j = 0;
2137 constraints_matrix[i_start_new + iy +
2138 (ix * deg)][j + j_start] =
2139 -constraints_matrix_old[i_start_old + ix +
2143 i_start_new += deg * deg;
2144 i_start_old += deg * deg;
2147 for (
unsigned int ix = 0; ix < deg; ++ix)
2151 constraints_matrix[i_start_new + ix][j +
2153 constraints_matrix_old[i_start_old + ix + deg]
2157 constraints_matrix[i_start_new + ix +
2159 constraints_matrix_old[i_start_old + ix][j];
2162 i_start_new += 2 * deg;
2163 i_start_old += 2 * deg;
2168 for (
unsigned int i = 0;
2172 constraints_matrix_old[i][j] = constraints_matrix
2179 unsigned int i_start =
2182 unsigned int j_start_new =
2184 unsigned int j_start_old = 0;
2186 for (
unsigned int block = 0; block < n_blocks; ++block)
2189 for (
unsigned int jx = 0; jx < deg; ++jx)
2191 for (
unsigned int jy = 0; jy < deg; ++jy)
2193 for (
unsigned int i = 0;
2197 constraints_matrix[i + i_start][j_start_new +
2200 constraints_matrix_old[i][j_start_old + jx +
2204 j_start_new += deg * deg;
2205 j_start_old += deg * deg;
2208 for (
unsigned int jx = 0; jx < deg; ++jx)
2210 for (
unsigned int jy = 0; jy < deg; ++jy)
2212 for (
unsigned int i = 0;
2216 constraints_matrix[i + i_start][j_start_new +
2219 -constraints_matrix_old[i][j_start_old +
2223 j_start_new += deg * deg;
2224 j_start_old += deg * deg;
2227 for (
unsigned int jx = 0; jx < deg; ++jx)
2229 for (
unsigned int i = 0;
2233 constraints_matrix[i + i_start][j_start_new +
2235 constraints_matrix_old[i][j_start_old + jx +
2237 constraints_matrix[i + i_start][j_start_new +
2239 constraints_matrix_old[i][j_start_old + jx];
2242 j_start_new += 2 * deg;
2243 j_start_old += 2 * deg;
2253 for (
unsigned int line = 0; line < n_internal_lines_on_children;
2260 constraints.
add_line(dofs_on_children[row_start + row]);
2261 for (
unsigned int i = 0; i < dofs_on_mother.size(); ++i)
2264 dofs_on_children[row_start + row],
2266 constraints_matrix[row_start + row][i]);
2269 dofs_on_children[row_start + row], 0.);
2274 for (
unsigned int line = 0; line < n_external_lines_on_children;
2277 unsigned int row_start =
2282 constraints.
add_line(dofs_on_children[row_start + row]);
2283 for (
unsigned int i = 0; i < dofs_on_mother.size(); ++i)
2286 dofs_on_children[row_start + row],
2288 constraints_matrix[row_start + row][i]);
2291 dofs_on_children[row_start + row], 0.);
2296 for (
unsigned int f = 0; f < n_children_per_face; ++f)
2298 unsigned int row_start =
2304 constraints.
add_line(dofs_on_children[row_start + row]);
2306 for (
unsigned int i = 0; i < dofs_on_mother.size(); ++i)
2309 dofs_on_children[row_start + row],
2311 constraints_matrix[row_start + row][i]);
2315 dofs_on_children[row_start + row], 0.);
2323 template <
int dim,
int spacedim,
typename number>
2341 std::vector<types::global_dof_index> primary_dofs;
2342 std::vector<types::global_dof_index> dependent_dofs;
2343 std::vector<types::global_dof_index> scratch_dofs;
2349 n_finite_elements(dof_handler), n_finite_elements(dof_handler));
2351 subface_interpolation_matrices(
2352 n_finite_elements(dof_handler),
2353 n_finite_elements(dof_handler),
2363 split_face_interpolation_matrices(n_finite_elements(dof_handler),
2364 n_finite_elements(dof_handler));
2370 n_finite_elements(dof_handler), n_finite_elements(dof_handler));
2381 if (cell->is_artificial())
2384 for (
const unsigned int face : cell->face_indices())
2385 if (cell->face(face)->has_children())
2390 if (cell->get_fe().n_dofs_per_face(face) == 0)
2393 Assert(cell->face(face)->refinement_case() ==
2404 Assert(cell->face(face)->n_active_fe_indices() == 1,
2406 Assert(cell->face(face)->fe_index_is_active(
2407 cell->active_fe_index()) ==
true,
2409 for (
unsigned int c = 0; c < cell->face(face)->n_children();
2411 if (!cell->neighbor_child_on_subface(face, c)
2413 Assert(cell->face(face)->child(c)->n_active_fe_indices() ==
2429 std::set<types::fe_index> fe_ind_face_subface;
2430 fe_ind_face_subface.insert(cell->active_fe_index());
2433 for (
unsigned int c = 0;
2434 c < cell->face(face)->n_active_descendants();
2437 const auto subcell =
2438 cell->neighbor_child_on_subface(face, c);
2439 if (!subcell->is_artificial())
2441 mother_face_dominates =
2442 mother_face_dominates &
2443 (cell->get_fe().compare_for_domination(
2444 subcell->get_fe(), 1));
2445 fe_ind_face_subface.insert(
2446 subcell->active_fe_index());
2450 switch (mother_face_dominates)
2462 primary_dofs.resize(
2463 cell->get_fe().n_dofs_per_face(face));
2465 cell->face(face)->get_dof_indices(
2466 primary_dofs, cell->active_fe_index());
2472 for (
unsigned int c = 0;
2473 c < cell->face(face)->n_children();
2476 if (cell->neighbor_child_on_subface(face, c)
2481 active_face_iterator subface =
2482 cell->face(face)->child(c);
2484 Assert(subface->n_active_fe_indices() == 1,
2488 subface->nth_active_fe_index(0);
2499 if (cell->get_fe().compare_for_domination(
2500 subface->
get_fe(subface_fe_index),
2507 dependent_dofs.resize(
2508 subface->
get_fe(subface_fe_index)
2510 subface->get_dof_indices(dependent_dofs,
2516 (void)dependent_dof;
2544 ensure_existence_of_subface_matrix(
2546 subface->
get_fe(subface_fe_index),
2548 subface_interpolation_matrices
2549 [cell->active_fe_index()][subface_fe_index][c]);
2553 filter_constraints(primary_dofs,
2555 *(subface_interpolation_matrices
2556 [cell->active_fe_index()]
2557 [subface_fe_index][c]),
2577 const ::hp::FECollection<dim, spacedim>
2607 {fe_ind_face_subface.begin(),
2608 fe_ind_face_subface.end()},
2614 "Could not find a least face dominating FE."));
2617 dof_handler.
get_fe(dominating_fe_index);
2622 cell->get_fe().n_dofs_per_face(face),
2625 ensure_existence_of_face_matrix(
2628 face_interpolation_matrices[dominating_fe_index]
2629 [cell->active_fe_index()]);
2633 ensure_existence_of_primary_dof_mask(
2636 (*face_interpolation_matrices
2637 [dominating_fe_index][cell->active_fe_index()]),
2638 primary_dof_masks[dominating_fe_index]
2639 [cell->active_fe_index()]);
2641 ensure_existence_of_split_face_matrix(
2642 *face_interpolation_matrices[dominating_fe_index]
2643 [cell->active_fe_index()],
2644 (*primary_dof_masks[dominating_fe_index]
2645 [cell->active_fe_index()]),
2646 split_face_interpolation_matrices
2647 [dominating_fe_index][cell->active_fe_index()]);
2650 &restrict_mother_to_virtual_primary_inv =
2651 (split_face_interpolation_matrices
2652 [dominating_fe_index][cell->active_fe_index()]
2656 &restrict_mother_to_virtual_dependent =
2657 (split_face_interpolation_matrices
2658 [dominating_fe_index][cell->active_fe_index()]
2663 constraint_matrix.reinit(
2664 cell->get_fe().n_dofs_per_face(face) -
2667 restrict_mother_to_virtual_dependent.
mmult(
2669 restrict_mother_to_virtual_primary_inv);
2673 scratch_dofs.resize(
2674 cell->get_fe().n_dofs_per_face(face));
2675 cell->face(face)->get_dof_indices(
2676 scratch_dofs, cell->active_fe_index());
2679 primary_dofs.clear();
2680 dependent_dofs.clear();
2681 for (
unsigned int i = 0;
2682 i < cell->get_fe().n_dofs_per_face(face);
2684 if ((*primary_dof_masks[dominating_fe_index]
2686 ->active_fe_index()])[i] ==
2688 primary_dofs.push_back(scratch_dofs[i]);
2690 dependent_dofs.push_back(scratch_dofs[i]);
2695 cell->get_fe().n_dofs_per_face(face) -
2698 filter_constraints(primary_dofs,
2707 for (
unsigned int sf = 0;
2708 sf < cell->face(face)->n_children();
2713 if (cell->neighbor_child_on_subface(face, sf)
2714 ->is_artificial() ||
2715 (dim == 2 && cell->is_ghost() &&
2716 cell->neighbor_child_on_subface(face, sf)
2722 ->n_active_fe_indices() == 1,
2726 cell->face(face)->child(sf)->nth_active_fe_index(
2729 dof_handler.
get_fe(subface_fe_index);
2736 ensure_existence_of_subface_matrix(
2740 subface_interpolation_matrices
2741 [dominating_fe_index][subface_fe_index][sf]);
2744 &restrict_subface_to_virtual = *(
2745 subface_interpolation_matrices
2746 [dominating_fe_index][subface_fe_index][sf]);
2748 constraint_matrix.reinit(
2752 restrict_subface_to_virtual.
mmult(
2754 restrict_mother_to_virtual_primary_inv);
2756 dependent_dofs.resize(
2758 cell->face(face)->child(sf)->get_dof_indices(
2759 dependent_dofs, subface_fe_index);
2761 filter_constraints(primary_dofs,
2784 Assert(cell->face(face)->fe_index_is_active(
2785 cell->active_fe_index()) ==
true,
2792 if (!cell->at_boundary(face) &&
2793 cell->neighbor(face)->is_artificial())
2799 !cell->face(face)->at_boundary() &&
2800 (cell->neighbor(face)->active_fe_index() !=
2801 cell->active_fe_index()) &&
2802 (!cell->face(face)->has_children() &&
2803 !cell->neighbor_is_coarser(face)))
2806 spacedim>::level_cell_iterator
2807 neighbor = cell->neighbor(face);
2811 cell->get_fe().compare_for_domination(neighbor->get_fe(),
2818 primary_dofs.resize(
2819 cell->get_fe().n_dofs_per_face(face));
2820 cell->face(face)->get_dof_indices(
2821 primary_dofs, cell->active_fe_index());
2826 if (primary_dofs.empty())
2829 dependent_dofs.resize(
2830 neighbor->get_fe().n_dofs_per_face(face));
2831 cell->face(face)->get_dof_indices(
2832 dependent_dofs, neighbor->active_fe_index());
2836 ensure_existence_of_face_matrix(
2839 face_interpolation_matrices
2840 [cell->active_fe_index()]
2841 [neighbor->active_fe_index()]);
2847 *(face_interpolation_matrices
2848 [cell->active_fe_index()]
2849 [neighbor->active_fe_index()]),
2892 if (cell < neighbor)
2907 cell->active_fe_index();
2909 neighbor->active_fe_index();
2910 std::set<types::fe_index> fes;
2911 fes.insert(this_fe_index);
2912 fes.insert(neighbor_fe_index);
2913 const ::hp::FECollection<dim, spacedim>
2919 {fes.begin(), fes.end()}, 1);
2924 "Could not find the dominating FE for " +
2925 cell->get_fe().get_name() +
" and " +
2926 neighbor->get_fe().get_name() +
2927 " inside FECollection."));
2930 fe_collection[dominating_fe_index];
2938 cell->get_fe().n_dofs_per_face(face),
2941 ensure_existence_of_face_matrix(
2944 face_interpolation_matrices
2945 [dominating_fe_index][cell->active_fe_index()]);
2949 ensure_existence_of_primary_dof_mask(
2952 (*face_interpolation_matrices
2953 [dominating_fe_index]
2954 [cell->active_fe_index()]),
2955 primary_dof_masks[dominating_fe_index]
2956 [cell->active_fe_index()]);
2958 ensure_existence_of_split_face_matrix(
2959 *face_interpolation_matrices
2960 [dominating_fe_index][cell->active_fe_index()],
2961 (*primary_dof_masks[dominating_fe_index]
2962 [cell->active_fe_index()]),
2963 split_face_interpolation_matrices
2964 [dominating_fe_index][cell->active_fe_index()]);
2967 double> &restrict_mother_to_virtual_primary_inv =
2968 (split_face_interpolation_matrices
2969 [dominating_fe_index][cell->active_fe_index()]
2973 double> &restrict_mother_to_virtual_dependent =
2974 (split_face_interpolation_matrices
2975 [dominating_fe_index][cell->active_fe_index()]
2980 constraint_matrix.reinit(
2981 cell->get_fe().n_dofs_per_face(face) -
2984 restrict_mother_to_virtual_dependent.
mmult(
2986 restrict_mother_to_virtual_primary_inv);
2990 scratch_dofs.resize(
2991 cell->get_fe().n_dofs_per_face(face));
2992 cell->face(face)->get_dof_indices(
2993 scratch_dofs, cell->active_fe_index());
2996 primary_dofs.clear();
2997 dependent_dofs.clear();
2998 for (
unsigned int i = 0;
2999 i < cell->get_fe().n_dofs_per_face(face);
3001 if ((*primary_dof_masks[dominating_fe_index]
3002 [cell->active_fe_index()])
3004 primary_dofs.push_back(scratch_dofs[i]);
3006 dependent_dofs.push_back(scratch_dofs[i]);
3012 dependent_dofs.size(),
3013 cell->get_fe().n_dofs_per_face(face) -
3016 filter_constraints(primary_dofs,
3025 neighbor->get_fe().n_dofs_per_face(face),
3028 ensure_existence_of_face_matrix(
3031 face_interpolation_matrices
3032 [dominating_fe_index]
3033 [neighbor->active_fe_index()]);
3036 &restrict_secondface_to_virtual =
3037 *(face_interpolation_matrices
3038 [dominating_fe_index]
3039 [neighbor->active_fe_index()]);
3041 constraint_matrix.reinit(
3042 neighbor->get_fe().n_dofs_per_face(face),
3045 restrict_secondface_to_virtual.
mmult(
3047 restrict_mother_to_virtual_primary_inv);
3049 dependent_dofs.resize(
3050 neighbor->get_fe().n_dofs_per_face(face));
3051 cell->face(face)->get_dof_indices(
3052 dependent_dofs, neighbor->active_fe_index());
3054 filter_constraints(primary_dofs,
3080 template <
int dim,
int spacedim,
typename number>
3087 "The given DoFHandler does not have any DoFs. Did you forget to "
3088 "call dof_handler.distribute_dofs()?"));
3098 dof_handler, constraints, std::integral_constant<int, dim>());
3103 dof_handler, constraints, std::integral_constant<int, dim>());
3110 template <
typename FaceIterator,
typename number>
3113 const FaceIterator &face_1,
3119 const number periodicity_factor,
3120 const unsigned int level)
3122 static const int dim = FaceIterator::AccessorType::dimension;
3123 static const int spacedim = FaceIterator::AccessorType::space_dimension;
3138 face_1->get_fe(face_1->nth_active_fe_index(0)).n_unique_faces(), 1);
3140 face_2->get_fe(face_2->nth_active_fe_index(0)).n_unique_faces(), 1);
3141 const unsigned int face_no = 0;
3149 if ((!use_mg) && face_2->has_children())
3151 Assert(face_2->n_children() ==
3155 const unsigned int dofs_per_face =
3156 face_1->get_fe(face_1->nth_active_fe_index(0))
3157 .n_dofs_per_face(face_no);
3161 std::vector<types::global_dof_index> dofs_1(dofs_per_face);
3162 face_1->get_dof_indices(dofs_1, face_1->nth_active_fe_index(0));
3163 for (
unsigned int i = 0; i < dofs_per_face; ++i)
3173 for (
unsigned int c = 0; c < face_2->n_children(); ++c)
3178 const auto &fe = face_1->get_fe(face_1->nth_active_fe_index(0));
3179 fe.get_subface_interpolation_matrix(fe,
3181 subface_interpolation,
3183 subface_interpolation.
mmult(child_transformation, transformation);
3187 child_transformation,
3190 combined_orientation,
3191 periodicity_factor);
3203 Assert(face_1->get_fe(face_1_index) == face_2->get_fe(face_2_index),
3205 "Matching periodic cells need to use the same finite element"));
3211 "The number of components in the mask has to be either "
3212 "zero or equal to the number of components in the finite "
3217 std::vector<types::global_dof_index> dofs_1(dofs_per_face);
3218 std::vector<types::global_dof_index> dofs_2(dofs_per_face);
3221 face_1->get_mg_dof_indices(
level, dofs_1, face_1_index);
3223 face_1->get_dof_indices(dofs_1, face_1_index);
3226 face_2->get_mg_dof_indices(
level, dofs_2, face_2_index);
3228 face_2->get_dof_indices(dofs_2, face_2_index);
3246 for (
unsigned int i = 0; i < dofs_per_face; ++i)
3255 if (
const auto tria =
dynamic_cast<
3257 &face_1->get_triangulation()))
3258 if (tria->with_artificial_cells() &&
3260 for (
unsigned int i = 0; i < dofs_per_face; ++i)
3284 std::map<unsigned int, unsigned int> cell_to_rotated_face_index;
3287 for (
unsigned int i = 0; i < dofs_per_face; ++i)
3301 boost::container::small_vector<
3302 std::pair<typename AffineConstraints<number>::size_type, number>,
3310 for (
unsigned int i = 0; i < dofs_per_face; ++i)
3331 bool is_identity_constrained =
false;
3333 number constraint_factor = periodicity_factor;
3335 constexpr double eps = 1.e-13;
3336 for (
unsigned int jj = 0; jj < dofs_per_face; ++jj)
3338 const auto entry = transformation(i, jj);
3341 if (is_identity_constrained)
3346 is_identity_constrained =
false;
3349 is_identity_constrained =
true;
3351 constraint_factor = entry * periodicity_factor;
3360 if (!is_identity_constrained)
3367 constraint_entries.clear();
3368 constraint_entries.reserve(dofs_per_face);
3370 for (
unsigned int jj = 0; jj < dofs_per_face; ++jj)
3374 const unsigned int j =
3376 jj, 0, combined_orientation)];
3378 if (
std::abs(transformation(i, jj)) > eps)
3379 constraint_entries.emplace_back(dofs_1[j],
3380 transformation(i, jj));
3397 const unsigned int j =
3399 target, 0, combined_orientation)];
3401 auto dof_left = dofs_1[j];
3402 auto dof_right = dofs_2[i];
3408 (dof_left < dof_right &&
3411 std::swap(dof_left, dof_right);
3412 constraint_factor = 1. / constraint_factor;
3428 bool constraints_are_cyclic =
true;
3429 number cycle_constraint_factor = constraint_factor;
3431 for (
auto test_dof = dof_right; test_dof != dof_left;)
3435 constraints_are_cyclic =
false;
3439 const auto &constraint_entries =
3441 if (constraint_entries.size() == 1)
3443 test_dof = constraint_entries[0].first;
3444 cycle_constraint_factor *= constraint_entries[0].second;
3448 constraints_are_cyclic =
false;
3463 if (constraints_are_cyclic)
3465 if (
std::abs(cycle_constraint_factor - number(1.)) > eps)
3471 dof_left, {{dof_right, constraint_factor}}, 0.);
3493 ExcMessage(
"The periodicity constraint is too large. "
3494 "The parameter periodicity_factor might "
3495 "be too large or too small."));
3509 template <
int dim,
int spacedim>
3511 compute_transformation(
3514 const std::vector<unsigned int> &first_vector_components)
3519 const unsigned int face_no = 0;
3525 if (matrix.m() == n_dofs_per_face)
3532 if (first_vector_components.empty() &&
matrix.m() == 0)
3549 using DoFTuple = std::array<unsigned int, spacedim>;
3554 for (
unsigned int i = 0; i < n_dofs_per_face; ++i)
3556 std::vector<unsigned int>::const_iterator comp_it =
3557 std::find(first_vector_components.begin(),
3558 first_vector_components.end(),
3560 if (comp_it != first_vector_components.end())
3562 const unsigned int first_vector_component = *comp_it;
3565 DoFTuple vector_dofs;
3567 unsigned int n_found = 1;
3572 "Error: the finite element does not have enough components "
3573 "to define rotated periodic boundaries."));
3575 for (
unsigned int k = 0; k < n_dofs_per_face; ++k)
3576 if ((k != i) && (quadrature.point(k) == quadrature.point(i)) &&
3578 first_vector_component) &&
3580 first_vector_component + spacedim))
3584 first_vector_component] = k;
3592 for (
unsigned int i = 0; i < spacedim; ++i)
3594 transformation[vector_dofs[i]][vector_dofs[i]] = 0.;
3595 for (
unsigned int j = 0; j < spacedim; ++j)
3596 transformation[vector_dofs[i]][vector_dofs[j]] =
3601 return transformation;
3609 template <
typename FaceIterator,
typename number>
3612 const FaceIterator &face_1,
3618 const std::vector<unsigned int> &first_vector_components,
3619 const number periodicity_factor)
3621 static const int dim = FaceIterator::AccessorType::dimension;
3622 static const int spacedim = FaceIterator::AccessorType::space_dimension;
3626 const auto [orientation, rotation, flip] =
3630 (orientation ==
true && flip ==
false && rotation ==
false),
3632 "The supplied orientation (orientation, rotation, flip) "
3633 "is invalid for 1d"));
3635 Assert((dim != 2) || (flip ==
false && rotation ==
false),
3637 "The supplied orientation (orientation, rotation, flip) "
3638 "is invalid for 2d"));
3641 ExcMessage(
"face_1 and face_2 are equal! Cannot constrain DoFs "
3642 "on the very same face"));
3644 Assert(face_1->at_boundary() && face_2->at_boundary(),
3645 ExcMessage(
"Faces for periodicity constraints must be on the "
3648 Assert(matrix.m() == matrix.n(),
3650 "The supplied (rotation or interpolation) matrix must "
3651 "be a square matrix"));
3653 Assert(first_vector_components.empty() || matrix.m() == spacedim,
3654 ExcMessage(
"first_vector_components is nonempty, so matrix must "
3655 "be a rotation matrix exactly of size spacedim"));
3657 if (!face_1->has_children())
3662 face_1->get_fe(face_1->nth_active_fe_index(0)).n_unique_faces(),
3664 const unsigned int face_no = 0;
3667 const unsigned int n_dofs_per_face =
3668 face_1->get_fe(face_1->nth_active_fe_index(0))
3669 .n_dofs_per_face(face_no);
3671 Assert(matrix.m() == 0 ||
3672 (first_vector_components.empty() &&
3673 matrix.m() == n_dofs_per_face) ||
3674 (!first_vector_components.empty() &&
3675 matrix.m() == spacedim),
3677 "The matrix must have either size 0 or spacedim "
3678 "(if first_vector_components is nonempty) "
3679 "or the size must be equal to the # of DoFs on the face "
3680 "(if first_vector_components is empty)."));
3683 if (!face_2->has_children())
3688 face_2->get_fe(face_2->nth_active_fe_index(0)).n_unique_faces(),
3690 const unsigned int face_no = 0;
3693 const unsigned int n_dofs_per_face =
3694 face_2->get_fe(face_2->nth_active_fe_index(0))
3695 .n_dofs_per_face(face_no);
3697 Assert(matrix.m() == 0 ||
3698 (first_vector_components.empty() &&
3699 matrix.m() == n_dofs_per_face) ||
3700 (!first_vector_components.empty() &&
3701 matrix.m() == spacedim),
3703 "The matrix must have either size 0 or spacedim "
3704 "(if first_vector_components is nonempty) "
3705 "or the size must be equal to the # of DoFs on the face "
3706 "(if first_vector_components is empty)."));
3710 if (face_1->has_children() && face_2->has_children())
3715 Assert(face_1->n_children() ==
3721 for (
unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face;
3727 const unsigned int face_no = dim == 2 ? 2 : 4;
3731 const auto reference_cell = ReferenceCells::get_hypercube<dim>();
3732 const unsigned int j =
3733 reference_cell.child_cell_on_face(face_no,
3735 combined_orientation);
3741 combined_orientation,
3743 first_vector_components,
3744 periodicity_factor);
3754 face_1->has_children() ?
3755 face_2->get_fe(face_2->nth_active_fe_index(0)) :
3756 face_1->get_fe(face_1->nth_active_fe_index(0));
3761 const unsigned int face_no = 0;
3767 if (n_dofs_per_face == 0)
3771 compute_transformation(fe, matrix, first_vector_components);
3773 if (!face_2->has_children())
3777 if (first_vector_components.empty() && matrix.m() == 0)
3784 combined_orientation,
3785 periodicity_factor);
3790 inverse.
invert(transformation);
3797 combined_orientation,
3798 periodicity_factor);
3812 const auto face_reference_cell = face_1->reference_cell();
3819 face_reference_cell.get_inverse_combined_orientation(
3820 combined_orientation),
3821 periodicity_factor);
3828 template <
int dim,
int spacedim,
typename number>
3835 const std::vector<unsigned int> &first_vector_components,
3836 const number periodicity_factor)
3839 for (
auto &pair : periodic_faces)
3842 const FaceIterator face_1 = pair.cell[0]->face(pair.face_idx[0]);
3843 const FaceIterator face_2 = pair.cell[1]->face(pair.face_idx[1]);
3845 Assert(face_1->at_boundary() && face_2->at_boundary(),
3858 first_vector_components,
3859 periodicity_factor);
3867 template <
int dim,
int spacedim,
typename number>
3872 const unsigned int direction,
3875 const number periodicity_factor)
3880 ExcMessage(
"The boundary indicators b_id1 and b_id2 must be "
3881 "different to denote different boundaries."));
3889 dof_handler, b_id1, b_id2, direction, matched_faces);
3891 make_periodicity_constraints<dim, spacedim>(matched_faces,
3894 std::vector<unsigned int>(),
3895 periodicity_factor);
3900 template <
int dim,
int spacedim,
typename number>
3904 const unsigned int direction,
3907 const number periodicity_factor)
3923 make_periodicity_constraints<dim, spacedim>(matched_faces,
3926 std::vector<unsigned int>(),
3927 periodicity_factor);
3941 template <
int dim,
int spacedim>
3946#ifdef DEAL_II_WITH_MPI
3947 std::vector<::LinearAlgebra::distributed::Vector<double>>
3962 template <
int dim,
int spacedim>
3964 compute_intergrid_weights_3(
3967 const unsigned int coarse_component,
4026 for (
unsigned int local_dof = 0; local_dof < copy_data.
dofs_per_cell;
4032 const unsigned int local_parameter_dof =
4041 coarse_to_fine_grid_map[cell]->set_dof_values_by_interpolation(
4042 parameter_dofs[local_parameter_dof],
4054 template <
int dim,
int spacedim>
4056 copy_intergrid_weights_3(
4057 const Assembler::CopyData<dim, spacedim> ©_data,
4058 const unsigned int coarse_component,
4060 const std::vector<types::global_dof_index> &weight_mapping,
4061 const bool is_called_in_parallel,
4062 std::vector<std::map<types::global_dof_index, float>> &weights)
4064 unsigned int pos = 0;
4065 for (
unsigned int local_dof = 0; local_dof < copy_data.dofs_per_cell;
4092 i < copy_data.global_parameter_representation[pos].size();
4098 if (copy_data.global_parameter_representation[pos](i) != 0)
4101 wi = copy_data.parameter_dof_indices[local_dof],
4102 wj = weight_mapping[i];
4104 copy_data.global_parameter_representation[pos](i);
4107 else if (!is_called_in_parallel)
4112 Assert(copy_data.global_parameter_representation[pos](i) ==
4128 template <
int dim,
int spacedim>
4130 compute_intergrid_weights_2(
4132 const unsigned int coarse_component,
4135 const std::vector<types::global_dof_index> &weight_mapping,
4136 std::vector<std::map<types::global_dof_index, float>> &weights)
4138 Assembler::CopyData<dim, spacedim> copy_data;
4140 unsigned int n_interesting_dofs = 0;
4141 for (
unsigned int local_dof = 0;
4142 local_dof < coarse_grid.
get_fe().n_dofs_per_cell();
4146 ++n_interesting_dofs;
4148 copy_data.global_parameter_representation.resize(n_interesting_dofs);
4150 bool is_called_in_parallel =
false;
4151 for (std::size_t i = 0;
4152 i < copy_data.global_parameter_representation.size();
4155#ifdef DEAL_II_WITH_MPI
4156 MPI_Comm communicator = MPI_COMM_SELF;
4159 const typename ::parallel::TriangulationBase<dim,
4161 &tria =
dynamic_cast<const typename ::parallel::
4162 TriangulationBase<dim, spacedim> &
>(
4163 coarse_to_fine_grid_map.get_destination_grid()
4164 .get_triangulation());
4165 communicator = tria.get_mpi_communicator();
4166 is_called_in_parallel =
true;
4168 catch (std::bad_cast &)
4174 const IndexSet locally_relevant_dofs =
4176 coarse_to_fine_grid_map.get_destination_grid());
4178 copy_data.global_parameter_representation[i].reinit(
4179 coarse_to_fine_grid_map.get_destination_grid()
4180 .locally_owned_dofs(),
4181 locally_relevant_dofs,
4185 copy_data.global_parameter_representation[i].reinit(n_fine_dofs);
4192 &coarse_to_fine_grid_map,
4196 const Assembler::Scratch &,
4197 Assembler::CopyData<dim, spacedim> ©_data) {
4198 compute_intergrid_weights_3<dim, spacedim>(cell,
4201 coarse_grid.get_fe(),
4202 coarse_to_fine_grid_map,
4210 is_called_in_parallel,
4211 &weights](
const Assembler::CopyData<dim, spacedim> ©_data) {
4212 copy_intergrid_weights_3<dim, spacedim>(copy_data,
4214 coarse_grid.get_fe(),
4216 is_called_in_parallel,
4224 Assembler::Scratch(),
4227#ifdef DEAL_II_WITH_MPI
4228 for (std::size_t i = 0;
4229 i < copy_data.global_parameter_representation.size();
4231 copy_data.global_parameter_representation[i].update_ghost_values();
4242 template <
int dim,
int spacedim>
4244 compute_intergrid_weights_1(
4246 const unsigned int coarse_component,
4248 const unsigned int fine_component,
4250 std::vector<std::map<types::global_dof_index, float>> &weights,
4251 std::vector<types::global_dof_index> &weight_mapping)
4255 &fine_fe = fine_grid.
get_fe();
4259 n_fine_dofs = fine_grid.
n_dofs();
4262 const unsigned int fine_dofs_per_cell = fine_fe.n_dofs_per_cell();
4266 const unsigned int coarse_dofs_per_cell_component =
4280 Assert(&coarse_to_fine_grid_map.get_source_grid() == &coarse_grid,
4282 Assert(&coarse_to_fine_grid_map.get_destination_grid() == &fine_grid,
4294 fine_fe.component_to_base_index(fine_component).first),
4301 for (
const auto &cell : coarse_grid.active_cell_iterators())
4327 std::vector<::Vector<double>> parameter_dofs(
4328 coarse_dofs_per_cell_component,
4333 for (
unsigned int local_coarse_dof = 0;
4334 local_coarse_dof < coarse_dofs_per_cell_component;
4336 for (
unsigned int fine_dof = 0; fine_dof < fine_fe.n_dofs_per_cell();
4338 if (fine_fe.system_to_component_index(fine_dof) ==
4339 std::make_pair(fine_component, local_coarse_dof))
4341 parameter_dofs[local_coarse_dof](fine_dof) = 1.;
4348 unsigned int n_parameters_on_fine_grid = 0;
4352 std::vector<bool> dof_is_interesting(fine_grid.
n_dofs(),
false);
4353 std::vector<types::global_dof_index> local_dof_indices(
4354 fine_fe.n_dofs_per_cell());
4356 for (
const auto &cell : fine_grid.active_cell_iterators() |
4359 cell->get_dof_indices(local_dof_indices);
4360 for (
unsigned int i = 0; i < fine_fe.n_dofs_per_cell(); ++i)
4361 if (fine_fe.system_to_component_index(i).first ==
4363 dof_is_interesting[local_dof_indices[i]] =
true;
4366 n_parameters_on_fine_grid = std::count(dof_is_interesting.begin(),
4367 dof_is_interesting.end(),
4374 weights.resize(n_coarse_dofs);
4376 weight_mapping.clear();
4380 std::vector<types::global_dof_index> local_dof_indices(
4381 fine_fe.n_dofs_per_cell());
4382 unsigned int next_free_index = 0;
4383 for (
const auto &cell : fine_grid.active_cell_iterators() |
4386 cell->get_dof_indices(local_dof_indices);
4387 for (
unsigned int i = 0; i < fine_fe.n_dofs_per_cell(); ++i)
4390 if ((fine_fe.system_to_component_index(i).first ==
4392 (weight_mapping[local_dof_indices[i]] ==
4395 weight_mapping[local_dof_indices[i]] = next_free_index;
4400 Assert(next_free_index == n_parameters_on_fine_grid,
4412 compute_intergrid_weights_2(coarse_grid,
4414 coarse_to_fine_grid_map,
4435 for (
unsigned int col = 0; col < n_parameters_on_fine_grid; ++col)
4440 if (weights[row].find(col) != weights[row].end())
4441 sum += weights[row][col];
4442 Assert((std::fabs(sum - 1) < 1.e-12) ||
4449 return n_parameters_on_fine_grid;
4458 template <
int dim,
int spacedim>
4462 const unsigned int coarse_component,
4464 const unsigned int fine_component,
4470 ExcMessage(
"This function is not yet implemented for DoFHandlers "
4471 "using hp-capabilities."));
4492 std::vector<std::map<types::global_dof_index, float>> weights;
4498 std::vector<types::global_dof_index> weight_mapping;
4500 const unsigned int n_parameters_on_fine_grid =
4501 internal::compute_intergrid_weights_1(coarse_grid,
4505 coarse_to_fine_grid_map,
4508 (void)n_parameters_on_fine_grid;
4512 n_fine_dofs = fine_grid.
n_dofs();
4520 mask[coarse_component] =
true;
4522 coarse_dof_is_parameter =
4523 extract_dofs<dim, spacedim>(coarse_grid,
ComponentMask(mask));
4535 std::vector<types::global_dof_index> representants(
4538 parameter_dof < n_coarse_dofs;
4540 if (coarse_dof_is_parameter.
is_element(parameter_dof))
4547 std::map<types::global_dof_index, float>::const_iterator i =
4548 weights[parameter_dof].begin();
4549 for (; i != weights[parameter_dof].end(); ++i)
4559 for (; global_dof < weight_mapping.size(); ++global_dof)
4560 if (weight_mapping[global_dof] ==
4566 representants[parameter_dof] = global_dof;
4582 std::vector<std::pair<types::global_dof_index, double>> constraint_line;
4601 std::map<types::global_dof_index, float>::const_iterator col_entry =
4603 for (; first_used_row < n_coarse_dofs; ++first_used_row)
4605 col_entry = weights[first_used_row].find(col);
4606 if (col_entry != weights[first_used_row].end())
4610 Assert(col_entry != weights[first_used_row].end(),
4613 if ((col_entry->second == 1) &&
4614 (representants[first_used_row] == global_dof))
4622 constraint_line.clear();
4624 row < n_coarse_dofs;
4627 const std::map<types::global_dof_index, float>::const_iterator j =
4628 weights[row].find(col);
4629 if ((j != weights[row].end()) && (j->second != 0))
4630 constraint_line.emplace_back(representants[row], j->second);
4639 template <
int dim,
int spacedim>
4643 const unsigned int coarse_component,
4645 const unsigned int fine_component,
4647 std::vector<std::map<types::global_dof_index, float>>
4648 &transfer_representation)
4652 ExcMessage(
"This function is not yet implemented for DoFHandlers "
4653 "using hp-capabilities."));
4674 std::vector<std::map<types::global_dof_index, float>> weights;
4680 std::vector<types::global_dof_index> weight_mapping;
4682 internal::compute_intergrid_weights_1(coarse_grid,
4686 coarse_to_fine_grid_map,
4692 std::count_if(weight_mapping.begin(),
4693 weight_mapping.end(),
4695 return dof != numbers::invalid_dof_index;
4699 std::vector<types::global_dof_index> inverse_weight_mapping(
4708 Assert((inverse_weight_mapping[parameter_dof] ==
4712 inverse_weight_mapping[parameter_dof] = i;
4719 transfer_representation.clear();
4720 transfer_representation.resize(n_rows);
4725 std::map<types::global_dof_index, float>::const_iterator j =
4727 for (; j != weights[i].end(); ++j)
4732 transfer_representation[p][i] = j->second;
4739 template <
int dim,
int spacedim,
typename number>
4748 ExcMessage(
"The number of components in the mask has to be either "
4749 "zero or equal to the number of components in the finite "
4758 std::vector<types::global_dof_index> face_dofs;
4761 std::vector<types::global_dof_index> cell_dofs;
4769 std::set<types::global_dof_index> dofs_already_treated;
4772 if (!cell->is_artificial() && cell->at_boundary())
4778 cell->get_dof_indices(cell_dofs);
4780 for (
const auto face_no : cell->face_indices())
4783 cell->face(face_no);
4787 if (face->at_boundary() &&
4789 (face->boundary_id() == boundary_id)))
4793 face->get_dof_indices(face_dofs, cell->active_fe_index());
4798 if (dofs_already_treated.find(face_dof) ==
4799 dofs_already_treated.end())
4803 const std::vector<types::global_dof_index>::iterator
4804 it_index_on_cell = std::find(cell_dofs.begin(),
4807 Assert(it_index_on_cell != cell_dofs.end(),
4809 const unsigned int index_on_cell =
4810 std::distance(cell_dofs.begin(), it_index_on_cell);
4812 cell->get_fe().get_nonzero_components(index_on_cell);
4815 for (
unsigned int c = 0; c < n_components; ++c)
4816 if (nonzero_component_array[c] && component_mask[c])
4832 Assert(zero_boundary_constraints
4833 .is_inhomogeneously_constrained(
4840 dofs_already_treated.insert(face_dof);
4849 template <
int dim,
int spacedim,
typename number>
4858 zero_boundary_constraints,
4869#include "dofs/dof_tools_constraints.inst"
void add_line(const size_type line_n)
void add_constraint(const size_type constrained_dof, const ArrayView< const std::pair< size_type, number > > &dependencies, const number inhomogeneity=0)
void add_entry(const size_type constrained_dof_index, const size_type column, const number weight)
const IndexSet & get_local_lines() const
void set_inhomogeneity(const size_type constrained_dof_index, const number value)
bool is_constrained(const size_type line_n) const
const std::vector< std::pair< size_type, number > > * get_constraint_entries(const size_type line_n) const
void constrain_dof_to_zero(const size_type constrained_dof)
bool represents_n_components(const unsigned int n) const
unsigned int n_selected_components(const unsigned int overall_number_of_components=numbers::invalid_unsigned_int) const
const hp::FECollection< dim, spacedim > & get_fe_collection() const
bool has_active_dofs() const
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) const
const Triangulation< dim, spacedim > & get_triangulation() const
bool has_hp_capabilities() const
types::global_dof_index n_dofs() const
unsigned int n_dofs_per_vertex() const
const unsigned int degree
unsigned int n_dofs_per_cell() const
unsigned int n_dofs_per_line() const
unsigned int n_dofs_per_face(unsigned int face_no=0, unsigned int child=0) const
unsigned int n_components() const
unsigned int n_unique_faces() const
unsigned int n_dofs_per_quad(unsigned int face_no=0) const
virtual std::string get_name() const =0
virtual const FiniteElement< dim, spacedim > & base_element(const unsigned int index) const
std::pair< unsigned int, unsigned int > component_to_base_index(const unsigned int component) const
const std::vector< Point< dim - 1 > > & get_unit_face_support_points(const unsigned int face_no=0) const
virtual void get_subface_interpolation_matrix(const FiniteElement< dim, spacedim > &source, const unsigned int subface, FullMatrix< double > &matrix, const unsigned int face_no=0) const
std::pair< unsigned int, unsigned int > system_to_component_index(const unsigned int index) const
const FullMatrix< double > & constraints(const ::internal::SubfaceCase< dim > &subface_case=::internal::SubfaceCase< dim >::case_isotropic) const
virtual unsigned int face_to_cell_index(const unsigned int face_dof_index, const unsigned int face, const types::geometric_orientation combined_orientation=numbers::default_geometric_orientation) const
virtual void get_face_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix, const unsigned int face_no=0) const
std::pair< unsigned int, unsigned int > face_system_to_component_index(const unsigned int index, const unsigned int face_no=0) const
void mmult(FullMatrix< number2 > &C, const FullMatrix< number2 > &B, const bool adding=false) const
void invert(const FullMatrix< number2 > &M)
bool is_element(const size_type index) const
bool get_anisotropic_refinement_flag() const
unsigned int n_cells() const
unsigned int size() const
unsigned int find_dominating_fe_extended(const std::set< unsigned int > &fes, const unsigned int codim=0) const
bool hp_constraints_are_implemented() const
unsigned int max_dofs_per_face() const
unsigned int n_components() const
unsigned int max_dofs_per_cell() const
#define DEAL_II_NAMESPACE_OPEN
constexpr bool running_in_debug_mode()
#define DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_ASSERT_UNREACHABLE()
IteratorRange< active_cell_iterator > active_cell_iterators() const
static ::ExceptionBase & ExcInvalidIterator()
static ::ExceptionBase & ExcGridNotCoarser()
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcFiniteElementsDontMatch()
static ::ExceptionBase & ExcNoComponentSelected()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcGridsDontMatch()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
typename ActiveSelector::cell_iterator cell_iterator
typename ActiveSelector::line_iterator line_iterator
typename ActiveSelector::face_iterator face_iterator
typename ActiveSelector::active_cell_iterator active_cell_iterator
void compute_intergrid_transfer_representation(const DoFHandler< dim, spacedim > &coarse_grid, const unsigned int coarse_component, const DoFHandler< dim, spacedim > &fine_grid, const unsigned int fine_component, const InterGridMap< DoFHandler< dim, spacedim > > &coarse_to_fine_grid_map, std::vector< std::map< types::global_dof_index, float > > &transfer_representation)
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
void compute_intergrid_constraints(const DoFHandler< dim, spacedim > &coarse_grid, const unsigned int coarse_component, const DoFHandler< dim, spacedim > &fine_grid, const unsigned int fine_component, const InterGridMap< DoFHandler< dim, spacedim > > &coarse_to_fine_grid_map, AffineConstraints< double > &constraints)
void make_zero_boundary_constraints(const DoFHandler< dim, spacedim > &dof, const types::boundary_id boundary_id, AffineConstraints< number > &zero_boundary_constraints, const ComponentMask &component_mask={})
Expression fabs(const Expression &x)
types::global_dof_index size_type
@ either_element_can_dominate
@ other_element_dominates
@ neither_element_dominates
@ matrix
Contents is actually a matrix.
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
T sum(const T &t, const MPI_Comm mpi_communicator)
void run(const std::vector< std::vector< Iterator > > &colored_iterators, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
std::tuple< bool, bool, bool > split_face_orientation(const types::geometric_orientation combined_face_orientation)
constexpr types::global_dof_index invalid_dof_index
constexpr unsigned int invalid_unsigned_int
constexpr types::boundary_id invalid_boundary_id
constexpr types::fe_index invalid_fe_index
typename type_identity< T >::type type_identity_t
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
unsigned short int fe_index
unsigned char geometric_orientation
static unsigned int n_children(const RefinementCase< dim > &refinement_case)