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,
1090 for (
unsigned int c = 0; c < cell->face(f)->n_children(); ++c)
1092 if (cell->neighbor_child_on_subface(f, c)->is_artificial())
1095 Assert(cell->face(f)->child(c)->n_active_fe_indices() == 1,
1098 Assert(cell->face(f)->child(c)->fe_index_is_active(
1099 cell->active_fe_index()) ==
true,
1108 face_dof_indices.resize(n_dofs);
1110 cell->face(f)->get_dof_indices(face_dof_indices);
1111 const std::vector<types::global_dof_index> dof_on_mother_face =
1114 cell->face(f)->child(0)->get_dof_indices(face_dof_indices);
1115 const std::vector<types::global_dof_index> dof_on_child_face_0 =
1118 cell->face(f)->child(1)->get_dof_indices(face_dof_indices);
1119 const std::vector<types::global_dof_index> dof_on_child_face_1 =
1128 const bool direction_mother = (cell->face(f)->vertex_index(0) >
1129 cell->face(f)->vertex_index(1)) ?
1132 const bool direction_child_0 =
1133 (cell->face(f)->child(0)->vertex_index(0) >
1134 cell->face(f)->child(0)->vertex_index(1)) ?
1137 const bool direction_child_1 =
1138 (cell->face(f)->child(1)->vertex_index(0) >
1139 cell->face(f)->child(1)->vertex_index(1)) ?
1143 for (
unsigned int row = 0; row < n_dofs; ++row)
1145 constraints.
add_line(dof_on_child_face_0[row]);
1146 constraints.
add_line(dof_on_child_face_1[row]);
1149 for (
unsigned int row = 0; row < n_dofs; ++row)
1151 for (
unsigned int dof_i_on_mother = 0;
1152 dof_i_on_mother < n_dofs;
1159 unsigned int shift_0 =
1160 (direction_mother == direction_child_0) ? 0 : n_dofs;
1161 constraints.
add_entry(dof_on_child_face_0[row],
1162 dof_on_mother_face[dof_i_on_mother],
1166 unsigned int shift_1 =
1167 (direction_mother == direction_child_1) ? 0 : n_dofs;
1168 constraints.
add_entry(dof_on_child_face_1[row],
1169 dof_on_mother_face[dof_i_on_mother],
1179 template <
int dim_,
int spacedim,
typename number>
1184 std::integral_constant<int, 3>)
1191 const unsigned int dim = 3;
1193 std::vector<types::global_dof_index> dofs_on_mother;
1194 std::vector<types::global_dof_index> dofs_on_children;
1200 if (cell->is_artificial())
1204 for (
const unsigned int face : cell->face_indices())
1207 if (cell->face(face)->has_children() ==
false)
1210 if (cell->get_fe().n_dofs_per_face(face) == 0)
1213 Assert(cell->face(face)->refinement_case() ==
1219 Assert(cell->face(face)->fe_index_is_active(
1220 cell->active_fe_index()) ==
true,
1225 for (
unsigned int c = 0; c < cell->face(face)->n_children(); ++c)
1227 if (cell->neighbor_child_on_subface(face, c)->is_artificial())
1231 cell->face(face)->child(c)->n_active_fe_indices(), 1);
1233 Assert(cell->face(face)->child(c)->fe_index_is_active(
1234 cell->active_fe_index()) ==
true,
1237 for (
unsigned int e = 0;
1238 e < GeometryInfo<dim>::vertices_per_face;
1244 ->n_active_fe_indices() == 1,
1248 cell->face(face)->child(c)->line(e)->fe_index_is_active(
1249 cell->active_fe_index()) ==
true,
1254 for (
unsigned int e = 0; e < GeometryInfo<dim>::vertices_per_face;
1257 Assert(cell->face(face)->line(e)->n_active_fe_indices() == 1,
1260 Assert(cell->face(face)->line(e)->fe_index_is_active(
1261 cell->active_fe_index()) ==
true,
1268 const unsigned int fe_index = cell->active_fe_index();
1271 unsigned int degree(fe.
degree);
1276 dofs_on_mother.resize(n_dofs_on_mother);
1278 const unsigned int n_lines_on_mother =
1293 const unsigned int n_internal_lines_on_children = 4;
1305 const unsigned int n_external_lines_on_children = 8;
1307 const unsigned int n_lines_on_children =
1308 n_internal_lines_on_children + n_external_lines_on_children;
1311 const unsigned int n_children_per_face =
1313 const unsigned int n_children_per_line =
1319 const unsigned int n_dofs_on_children =
1323 dofs_on_children.clear();
1324 dofs_on_children.reserve(n_dofs_on_children);
1334 unsigned int next_index = 0;
1340 for (
unsigned int line = 0;
1341 line < GeometryInfo<dim>::lines_per_face;
1344 dofs_on_mother[next_index++] =
1345 this_face->line(line)->dof_index(dof, fe_index);
1349 dofs_on_mother[next_index++] =
1350 this_face->dof_index(dof, fe_index);
1369 dofs_on_children.push_back(
1370 this_face->child(0)->line(1)->dof_index(dof, fe_index));
1373 dofs_on_children.push_back(
1374 this_face->child(2)->line(1)->dof_index(dof, fe_index));
1377 dofs_on_children.push_back(
1378 this_face->child(0)->line(3)->dof_index(dof, fe_index));
1381 dofs_on_children.push_back(
1382 this_face->child(1)->line(3)->dof_index(dof, fe_index));
1387 for (
unsigned int line = 0;
1388 line < GeometryInfo<dim>::lines_per_face;
1390 for (
unsigned int child = 0; child < n_children_per_line;
1393 dofs_on_children.push_back(
1394 this_face->line(line)->child(child)->dof_index(dof,
1398 for (
unsigned int child = 0; child < n_children_per_face; ++child)
1401 if (cell->neighbor_child_on_subface(face, child)
1407 dofs_on_children.push_back(
1408 this_face->child(child)->dof_index(dof, fe_index));
1413 Assert(dofs_on_children.size() <= n_dofs_on_children,
1423 std::vector<bool> direction_mother(
1425 for (
unsigned int line = 0;
1426 line < GeometryInfo<dim>::lines_per_face;
1428 if (this_face->line(line)->vertex_index(0) >
1429 this_face->line(line)->vertex_index(1))
1430 direction_mother[line] =
true;
1433 std::vector<bool> direction_child_intern(
1434 n_internal_lines_on_children,
false);
1439 unsigned int center = this_face->child(0)->vertex_index(3);
1442 for (
unsigned int line = 0; line < n_internal_lines_on_children;
1446 direction_child_intern[line] =
1447 this_face->line(line)->child(0)->vertex_index(1) <
1454 direction_child_intern[line] =
1455 this_face->line(line)->child(0)->vertex_index(1) >
1462 std::vector<bool> direction_child(n_external_lines_on_children,
1464 for (
unsigned int line = 0;
1465 line < GeometryInfo<dim>::lines_per_face;
1468 if (this_face->line(line)->child(0)->vertex_index(0) >
1469 this_face->line(line)->child(0)->vertex_index(1))
1470 direction_child[2 * line] =
true;
1471 if (this_face->line(line)->child(1)->vertex_index(0) >
1472 this_face->line(line)->child(1)->vertex_index(1))
1473 direction_child[2 * line + 1] =
true;
1478 bool mother_flip_x =
false;
1479 bool mother_flip_y =
false;
1480 bool mother_flip_xy =
false;
1481 std::vector<bool> child_flip_x(n_children_per_face,
false);
1482 std::vector<bool> child_flip_y(n_children_per_face,
false);
1483 std::vector<bool> child_flip_xy(n_children_per_face,
false);
1486 [2] = {{1, 2}, {0, 3}, {3, 0}, {2, 1}};
1491 unsigned int current_glob = cell->face(face)->vertex_index(0);
1492 unsigned int current_max = 0;
1493 for (
unsigned int v = 1;
1494 v < GeometryInfo<dim>::vertices_per_face;
1496 if (current_glob < this_face->vertex_index(v))
1499 current_glob = this_face->vertex_index(v);
1504 if (current_max < 2)
1505 mother_flip_y =
true;
1509 if (current_max % 2 == 0)
1510 mother_flip_x =
true;
1513 if (this_face->vertex_index(
1514 vertices_adjacent_on_face[current_max][0]) <
1515 this_face->vertex_index(
1516 vertices_adjacent_on_face[current_max][1]))
1517 mother_flip_xy =
true;
1522 for (
unsigned int child = 0; child < n_children_per_face; ++child)
1524 unsigned int current_max = 0;
1525 unsigned int current_glob =
1526 this_face->child(child)->vertex_index(0);
1528 for (
unsigned int v = 1;
1529 v < GeometryInfo<dim>::vertices_per_face;
1531 if (current_glob < this_face->child(child)->vertex_index(v))
1534 current_glob = this_face->child(child)->vertex_index(v);
1537 if (current_max < 2)
1538 child_flip_y[child] =
true;
1540 if (current_max % 2 == 0)
1541 child_flip_x[child] =
true;
1543 if (this_face->child(child)->vertex_index(
1544 vertices_adjacent_on_face[current_max][0]) <
1545 this_face->child(child)->vertex_index(
1546 vertices_adjacent_on_face[current_max][1]))
1547 child_flip_xy[child] =
true;
1549 child_flip_xy[child] = mother_flip_xy;
1553 std::vector<std::vector<double>> constraints_matrix(
1556 std::vector<double>(dofs_on_mother.size(), 0));
1561 for (
unsigned int line = 0; line < n_internal_lines_on_children;
1565 unsigned int line_mother = line / 2;
1566 unsigned int row_mother =
1570 for (
unsigned int i = 0;
1573 constraints_matrix[row + row_start][i] =
1577 for (
unsigned int line = 0; line < n_internal_lines_on_children;
1581 unsigned int line_mother = line / 2;
1582 unsigned int row_mother =
1586 for (
unsigned int i =
1588 i < dofs_on_mother.size();
1590 constraints_matrix[row + row_start][i] =
1595 unsigned int row_offset =
1597 for (
unsigned int line = 0; line < n_external_lines_on_children;
1601 unsigned int line_mother = line / 2;
1602 unsigned int row_mother =
1604 for (
unsigned int row = row_offset;
1607 for (
unsigned int i = 0; i < dofs_on_mother.size(); ++i)
1608 constraints_matrix[row + row_start][i] =
1614 for (
unsigned int face = 0; face < n_children_per_face; ++face)
1617 for (
unsigned int row = row_offset;
1620 for (
unsigned int i = 0; i < dofs_on_mother.size(); ++i)
1621 constraints_matrix[row + row_start][i] =
1630 for (
unsigned int i = 0;
1635 unsigned int tmp_i = i % degree;
1638 for (
unsigned int j = 0;
1643 unsigned int tmp_j = j % degree;
1645 if ((line_i < 2 && line_j < 2) ||
1646 (line_i >= 2 && line_j >= 2))
1648 if (direction_child_intern[line_i] !=
1649 direction_mother[line_j])
1651 if ((tmp_i + tmp_j) % 2 == 1)
1653 constraints_matrix[i][j] *= -1.0;
1659 if (direction_mother[line_i])
1661 if ((tmp_i + tmp_j) % 2 == 1)
1663 constraints_matrix[i][j] *= -1.0;
1671 for (
unsigned int i =
1677 unsigned int tmp_i = i % degree;
1680 for (
unsigned int j = 0;
1685 unsigned int tmp_j = j % degree;
1687 if (direction_child[line_i] != direction_mother[line_j])
1689 if ((tmp_i + tmp_j) % 2 == 1)
1691 constraints_matrix[i][j] *= -1.0;
1709 unsigned int tmp_i = i % degree;
1711 unsigned int start_j =
1714 for (
unsigned int block = 0; block < n_blocks; ++block)
1717 for (
unsigned int jy = 0; jy < degree - 1; ++jy)
1718 for (
unsigned int jx = 0; jx < degree - 1; ++jx)
1720 unsigned int j = start_j + jx + (jy * (degree - 1));
1721 if (direction_child_intern[line_i] != mother_flip_y)
1723 if ((jy + tmp_i) % 2 == 0)
1725 constraints_matrix[i][j] *= -1.0;
1730 start_j += (degree - 1) * (degree - 1);
1733 for (
unsigned int jy = 0; jy < degree - 1; ++jy)
1734 for (
unsigned int jx = 0; jx < degree - 1; ++jx)
1736 unsigned int j = start_j + jx + (jy * (degree - 1));
1738 if (direction_child_intern[line_i] != mother_flip_y)
1740 if ((jy + tmp_i) % 2 == 0)
1742 constraints_matrix[i][j] *= -1.0;
1746 start_j += (degree - 1) * (degree - 1);
1750 start_j += degree - 1;
1754 start_j += degree - 1;
1764 unsigned int tmp_i = i % degree;
1766 unsigned int start_j =
1769 for (
unsigned int block = 0; block < n_blocks; block++)
1772 for (
unsigned int jy = 0; jy < degree - 1; ++jy)
1773 for (
unsigned int jx = 0; jx < degree - 1; ++jx)
1775 unsigned int j = start_j + jx + (jy * (degree - 1));
1776 if (direction_child_intern[line_i] != mother_flip_x)
1778 if ((jx + tmp_i) % 2 == 0)
1780 constraints_matrix[i][j] *= -1.0;
1785 start_j += (degree - 1) * (degree - 1);
1788 for (
unsigned int jy = 0; jy < degree - 1; ++jy)
1789 for (
unsigned int jx = 0; jx < degree - 1; ++jx)
1791 unsigned int j = start_j + jx + (jy * (degree - 1));
1792 if (direction_child_intern[line_i] != mother_flip_x)
1794 if ((jx + tmp_i) % 2 == 0)
1796 constraints_matrix[i][j] *= -1.0;
1800 start_j += (degree - 1) * (degree - 1);
1804 start_j += degree - 1;
1808 start_j += degree - 1;
1813 unsigned int degree_square = (degree - 1) * (degree - 1);
1817 for (
unsigned int child_face = 0;
1818 child_face < n_children_per_face;
1820 for (
unsigned int block = 0; block < n_blocks; ++block)
1832 for (
unsigned int iy = 0; iy < degree - 1; ++iy)
1833 for (
unsigned int ix = 0; ix < degree - 1; ++ix)
1839 for (
unsigned int jy = 0; jy < degree - 1; ++jy)
1840 for (
unsigned int jx = 0; jx < degree - 1; ++jx)
1842 if (child_flip_x[child_face] !=
1845 if ((ix + jx) % 2 == 1)
1847 constraints_matrix[i][j] *= -1.0;
1851 if (child_flip_y[child_face] !=
1854 if ((iy + jy) % 2 == 1)
1856 constraints_matrix[i][j] *= -1.0;
1866 for (
unsigned int iy = 0; iy < degree - 1; ++iy)
1867 for (
unsigned int ix = 0; ix < degree - 1; ++ix)
1872 degree_square + block * block_size;
1873 for (
unsigned int jy = 0; jy < degree - 1; ++jy)
1874 for (
unsigned int jx = 0; jx < degree - 1; ++jx)
1876 if (child_flip_x[child_face] !=
1879 if ((ix + jx) % 2 == 1)
1881 constraints_matrix[i][j] *= -1.0;
1885 if (child_flip_y[child_face] !=
1888 if ((iy + jy) % 2 == 1)
1890 constraints_matrix[i][j] *= -1.0;
1902 for (
unsigned int iy = 0; iy < degree - 1; ++iy)
1907 degree_square + block * block_size;
1908 for (
unsigned int jy = 0; jy < degree - 1; ++jy)
1909 for (
unsigned int jx = 0; jx < degree - 1; ++jx)
1911 if (child_flip_x[child_face] !=
1916 constraints_matrix[i][j] *= -1.0;
1920 if (child_flip_y[child_face] !=
1923 if ((iy + jy) % 2 == 1)
1925 constraints_matrix[i][j] *= -1.0;
1934 2 * degree_square + block * block_size;
1935 for (
unsigned int jy = 0; jy < degree - 1; ++jy)
1937 if (child_flip_y[child_face] !=
1940 if ((iy + jy) % 2 == 1)
1942 constraints_matrix[i][j] *= -1.0;
1952 for (
unsigned int ix = 0; ix < degree - 1; ++ix)
1957 degree_square + block * block_size;
1958 for (
unsigned int jy = 0; jy < degree - 1; ++jy)
1959 for (
unsigned int jx = 0; jx < degree - 1; ++jx)
1961 if (child_flip_x[child_face] !=
1964 if ((ix + jx) % 2 == 1)
1966 constraints_matrix[i][j] *= -1.0;
1970 if (child_flip_y[child_face] !=
1975 constraints_matrix[i][j] *= -1.0;
1984 2 * degree_square + (degree - 1) +
1986 for (
unsigned int jx = 0; jx < degree - 1; ++jx)
1988 if (child_flip_x[child_face] !=
1991 if ((ix + jx) % 2 == 1)
1993 constraints_matrix[i][j] *= -1.0;
2010 for (
unsigned int i = 0;
2018 std::vector<double> constraints_matrix_old(
2019 dofs_on_mother.size(), 0);
2020 for (
unsigned int j = 0; j < dofs_on_mother.size(); ++j)
2022 constraints_matrix_old[j] = constraints_matrix[i][j];
2025 unsigned int j_start =
2027 for (
unsigned block = 0; block < n_blocks; block++)
2030 for (
unsigned int jy = 0; jy < degree - 1; ++jy)
2031 for (
unsigned int jx = 0; jx < degree - 1; ++jx)
2033 unsigned int j_old =
2034 j_start + jx + (jy * (degree - 1));
2035 unsigned int j_new =
2036 j_start + jy + (jx * (degree - 1));
2037 constraints_matrix[i][j_new] =
2038 constraints_matrix_old[j_old];
2040 j_start += degree_square;
2043 for (
unsigned int jy = 0; jy < degree - 1; ++jy)
2044 for (
unsigned int jx = 0; jx < degree - 1; ++jx)
2046 unsigned int j_old =
2047 j_start + jx + (jy * (degree - 1));
2048 unsigned int j_new =
2049 j_start + jy + (jx * (degree - 1));
2050 constraints_matrix[i][j_new] =
2051 -constraints_matrix_old[j_old];
2053 j_start += degree_square;
2056 for (
unsigned int j = j_start;
2057 j < j_start + (degree - 1);
2060 constraints_matrix[i][j] =
2061 constraints_matrix_old[j + (degree - 1)];
2062 constraints_matrix[i][j + (degree - 1)] =
2063 constraints_matrix_old[j];
2065 j_start += 2 * (degree - 1);
2072 const unsigned int deg = degree - 1;
2075 std::vector<std::vector<double>> constraints_matrix_old(
2078 for (
unsigned int i = 0;
2082 constraints_matrix_old[i][j] = constraints_matrix
2087 for (
unsigned int child = 0; child < n_children_per_face;
2090 if (!child_flip_xy[child])
2093 unsigned int i_start_new =
2098 unsigned int j_start =
2101 for (
unsigned int block = 0; block < n_blocks; block++)
2104 for (
unsigned int ix = 0; ix < deg; ++ix)
2106 for (
unsigned int iy = 0; iy < deg; ++iy)
2108 for (
unsigned int j = 0;
2111 constraints_matrix[i_start_new + iy +
2112 (ix * deg)][j + j_start] =
2113 constraints_matrix_old[i_start_old + ix +
2117 i_start_new += deg * deg;
2118 i_start_old += deg * deg;
2121 for (
unsigned int ix = 0; ix < deg; ++ix)
2123 for (
unsigned int iy = 0; iy < deg; ++iy)
2125 for (
unsigned int j = 0;
2128 constraints_matrix[i_start_new + iy +
2129 (ix * deg)][j + j_start] =
2130 -constraints_matrix_old[i_start_old + ix +
2134 i_start_new += deg * deg;
2135 i_start_old += deg * deg;
2138 for (
unsigned int ix = 0; ix < deg; ++ix)
2142 constraints_matrix[i_start_new + ix][j +
2144 constraints_matrix_old[i_start_old + ix + deg]
2148 constraints_matrix[i_start_new + ix +
2150 constraints_matrix_old[i_start_old + ix][j];
2153 i_start_new += 2 * deg;
2154 i_start_old += 2 * deg;
2159 for (
unsigned int i = 0;
2163 constraints_matrix_old[i][j] = constraints_matrix
2170 unsigned int i_start =
2173 unsigned int j_start_new =
2175 unsigned int j_start_old = 0;
2177 for (
unsigned int block = 0; block < n_blocks; ++block)
2180 for (
unsigned int jx = 0; jx < deg; ++jx)
2182 for (
unsigned int jy = 0; jy < deg; ++jy)
2184 for (
unsigned int i = 0;
2188 constraints_matrix[i + i_start][j_start_new +
2191 constraints_matrix_old[i][j_start_old + jx +
2195 j_start_new += deg * deg;
2196 j_start_old += deg * deg;
2199 for (
unsigned int jx = 0; jx < deg; ++jx)
2201 for (
unsigned int jy = 0; jy < deg; ++jy)
2203 for (
unsigned int i = 0;
2207 constraints_matrix[i + i_start][j_start_new +
2210 -constraints_matrix_old[i][j_start_old +
2214 j_start_new += deg * deg;
2215 j_start_old += deg * deg;
2218 for (
unsigned int jx = 0; jx < deg; ++jx)
2220 for (
unsigned int i = 0;
2224 constraints_matrix[i + i_start][j_start_new +
2226 constraints_matrix_old[i][j_start_old + jx +
2228 constraints_matrix[i + i_start][j_start_new +
2230 constraints_matrix_old[i][j_start_old + jx];
2233 j_start_new += 2 * deg;
2234 j_start_old += 2 * deg;
2244 for (
unsigned int line = 0; line < n_internal_lines_on_children;
2251 constraints.
add_line(dofs_on_children[row_start + row]);
2252 for (
unsigned int i = 0; i < dofs_on_mother.size(); ++i)
2255 dofs_on_children[row_start + row],
2257 constraints_matrix[row_start + row][i]);
2260 dofs_on_children[row_start + row], 0.);
2265 for (
unsigned int line = 0; line < n_external_lines_on_children;
2268 unsigned int row_start =
2273 constraints.
add_line(dofs_on_children[row_start + row]);
2274 for (
unsigned int i = 0; i < dofs_on_mother.size(); ++i)
2277 dofs_on_children[row_start + row],
2279 constraints_matrix[row_start + row][i]);
2282 dofs_on_children[row_start + row], 0.);
2287 for (
unsigned int f = 0; f < n_children_per_face; ++f)
2289 unsigned int row_start =
2295 constraints.
add_line(dofs_on_children[row_start + row]);
2297 for (
unsigned int i = 0; i < dofs_on_mother.size(); ++i)
2300 dofs_on_children[row_start + row],
2302 constraints_matrix[row_start + row][i]);
2306 dofs_on_children[row_start + row], 0.);
2314 template <
int dim,
int spacedim,
typename number>
2332 std::vector<types::global_dof_index> primary_dofs;
2333 std::vector<types::global_dof_index> dependent_dofs;
2334 std::vector<types::global_dof_index> scratch_dofs;
2340 n_finite_elements(dof_handler), n_finite_elements(dof_handler));
2342 subface_interpolation_matrices(
2343 n_finite_elements(dof_handler),
2344 n_finite_elements(dof_handler),
2354 split_face_interpolation_matrices(n_finite_elements(dof_handler),
2355 n_finite_elements(dof_handler));
2361 n_finite_elements(dof_handler), n_finite_elements(dof_handler));
2372 if (cell->is_artificial())
2375 for (
const unsigned int face : cell->face_indices())
2376 if (cell->face(face)->has_children())
2381 if (cell->get_fe().n_dofs_per_face(face) == 0)
2384 Assert(cell->face(face)->refinement_case() ==
2395 Assert(cell->face(face)->n_active_fe_indices() == 1,
2397 Assert(cell->face(face)->fe_index_is_active(
2398 cell->active_fe_index()) ==
true,
2400 for (
unsigned int c = 0; c < cell->face(face)->n_children();
2402 if (!cell->neighbor_child_on_subface(face, c)
2404 Assert(cell->face(face)->child(c)->n_active_fe_indices() ==
2420 std::set<types::fe_index> fe_ind_face_subface;
2421 fe_ind_face_subface.insert(cell->active_fe_index());
2424 for (
unsigned int c = 0;
2425 c < cell->face(face)->n_active_descendants();
2428 const auto subcell =
2429 cell->neighbor_child_on_subface(face, c);
2430 if (!subcell->is_artificial())
2432 mother_face_dominates =
2433 mother_face_dominates &
2434 (cell->get_fe().compare_for_domination(
2435 subcell->get_fe(), 1));
2436 fe_ind_face_subface.insert(
2437 subcell->active_fe_index());
2441 switch (mother_face_dominates)
2453 primary_dofs.resize(
2454 cell->get_fe().n_dofs_per_face(face));
2456 cell->face(face)->get_dof_indices(
2457 primary_dofs, cell->active_fe_index());
2463 for (
unsigned int c = 0;
2464 c < cell->face(face)->n_children();
2467 if (cell->neighbor_child_on_subface(face, c)
2472 active_face_iterator subface =
2473 cell->face(face)->child(c);
2475 Assert(subface->n_active_fe_indices() == 1,
2479 subface->nth_active_fe_index(0);
2490 if (cell->get_fe().compare_for_domination(
2491 subface->
get_fe(subface_fe_index),
2498 dependent_dofs.resize(
2499 subface->
get_fe(subface_fe_index)
2501 subface->get_dof_indices(dependent_dofs,
2507 (void)dependent_dof;
2535 ensure_existence_of_subface_matrix(
2537 subface->
get_fe(subface_fe_index),
2539 subface_interpolation_matrices
2540 [cell->active_fe_index()][subface_fe_index][c]);
2544 filter_constraints(primary_dofs,
2546 *(subface_interpolation_matrices
2547 [cell->active_fe_index()]
2548 [subface_fe_index][c]),
2568 const ::hp::FECollection<dim, spacedim>
2598 {fe_ind_face_subface.begin(),
2599 fe_ind_face_subface.end()},
2605 "Could not find a least face dominating FE."));
2608 dof_handler.
get_fe(dominating_fe_index);
2613 cell->get_fe().n_dofs_per_face(face),
2616 ensure_existence_of_face_matrix(
2619 face_interpolation_matrices[dominating_fe_index]
2620 [cell->active_fe_index()]);
2624 ensure_existence_of_primary_dof_mask(
2627 (*face_interpolation_matrices
2628 [dominating_fe_index][cell->active_fe_index()]),
2629 primary_dof_masks[dominating_fe_index]
2630 [cell->active_fe_index()]);
2632 ensure_existence_of_split_face_matrix(
2633 *face_interpolation_matrices[dominating_fe_index]
2634 [cell->active_fe_index()],
2635 (*primary_dof_masks[dominating_fe_index]
2636 [cell->active_fe_index()]),
2637 split_face_interpolation_matrices
2638 [dominating_fe_index][cell->active_fe_index()]);
2641 &restrict_mother_to_virtual_primary_inv =
2642 (split_face_interpolation_matrices
2643 [dominating_fe_index][cell->active_fe_index()]
2647 &restrict_mother_to_virtual_dependent =
2648 (split_face_interpolation_matrices
2649 [dominating_fe_index][cell->active_fe_index()]
2654 constraint_matrix.reinit(
2655 cell->get_fe().n_dofs_per_face(face) -
2658 restrict_mother_to_virtual_dependent.
mmult(
2660 restrict_mother_to_virtual_primary_inv);
2664 scratch_dofs.resize(
2665 cell->get_fe().n_dofs_per_face(face));
2666 cell->face(face)->get_dof_indices(
2667 scratch_dofs, cell->active_fe_index());
2670 primary_dofs.clear();
2671 dependent_dofs.clear();
2672 for (
unsigned int i = 0;
2673 i < cell->get_fe().n_dofs_per_face(face);
2675 if ((*primary_dof_masks[dominating_fe_index]
2677 ->active_fe_index()])[i] ==
2679 primary_dofs.push_back(scratch_dofs[i]);
2681 dependent_dofs.push_back(scratch_dofs[i]);
2686 cell->get_fe().n_dofs_per_face(face) -
2689 filter_constraints(primary_dofs,
2698 for (
unsigned int sf = 0;
2699 sf < cell->face(face)->n_children();
2704 if (cell->neighbor_child_on_subface(face, sf)
2705 ->is_artificial() ||
2706 (dim == 2 && cell->is_ghost() &&
2707 cell->neighbor_child_on_subface(face, sf)
2713 ->n_active_fe_indices() == 1,
2717 cell->face(face)->child(sf)->nth_active_fe_index(
2720 dof_handler.
get_fe(subface_fe_index);
2727 ensure_existence_of_subface_matrix(
2731 subface_interpolation_matrices
2732 [dominating_fe_index][subface_fe_index][sf]);
2735 &restrict_subface_to_virtual = *(
2736 subface_interpolation_matrices
2737 [dominating_fe_index][subface_fe_index][sf]);
2739 constraint_matrix.reinit(
2743 restrict_subface_to_virtual.
mmult(
2745 restrict_mother_to_virtual_primary_inv);
2747 dependent_dofs.resize(
2749 cell->face(face)->child(sf)->get_dof_indices(
2750 dependent_dofs, subface_fe_index);
2752 filter_constraints(primary_dofs,
2775 Assert(cell->face(face)->fe_index_is_active(
2776 cell->active_fe_index()) ==
true,
2783 if (!cell->at_boundary(face) &&
2784 cell->neighbor(face)->is_artificial())
2790 !cell->face(face)->at_boundary() &&
2791 (cell->neighbor(face)->active_fe_index() !=
2792 cell->active_fe_index()) &&
2793 (!cell->face(face)->has_children() &&
2794 !cell->neighbor_is_coarser(face)))
2797 spacedim>::level_cell_iterator
2798 neighbor = cell->neighbor(face);
2802 cell->get_fe().compare_for_domination(neighbor->get_fe(),
2809 primary_dofs.resize(
2810 cell->get_fe().n_dofs_per_face(face));
2811 cell->face(face)->get_dof_indices(
2812 primary_dofs, cell->active_fe_index());
2817 if (primary_dofs.empty())
2820 dependent_dofs.resize(
2821 neighbor->get_fe().n_dofs_per_face(face));
2822 cell->face(face)->get_dof_indices(
2823 dependent_dofs, neighbor->active_fe_index());
2827 ensure_existence_of_face_matrix(
2830 face_interpolation_matrices
2831 [cell->active_fe_index()]
2832 [neighbor->active_fe_index()]);
2838 *(face_interpolation_matrices
2839 [cell->active_fe_index()]
2840 [neighbor->active_fe_index()]),
2883 if (cell < neighbor)
2898 cell->active_fe_index();
2900 neighbor->active_fe_index();
2901 std::set<types::fe_index> fes;
2902 fes.insert(this_fe_index);
2903 fes.insert(neighbor_fe_index);
2904 const ::hp::FECollection<dim, spacedim>
2910 {fes.begin(), fes.end()}, 1);
2915 "Could not find the dominating FE for " +
2916 cell->get_fe().get_name() +
" and " +
2917 neighbor->get_fe().get_name() +
2918 " inside FECollection."));
2921 fe_collection[dominating_fe_index];
2929 cell->get_fe().n_dofs_per_face(face),
2932 ensure_existence_of_face_matrix(
2935 face_interpolation_matrices
2936 [dominating_fe_index][cell->active_fe_index()]);
2940 ensure_existence_of_primary_dof_mask(
2943 (*face_interpolation_matrices
2944 [dominating_fe_index]
2945 [cell->active_fe_index()]),
2946 primary_dof_masks[dominating_fe_index]
2947 [cell->active_fe_index()]);
2949 ensure_existence_of_split_face_matrix(
2950 *face_interpolation_matrices
2951 [dominating_fe_index][cell->active_fe_index()],
2952 (*primary_dof_masks[dominating_fe_index]
2953 [cell->active_fe_index()]),
2954 split_face_interpolation_matrices
2955 [dominating_fe_index][cell->active_fe_index()]);
2958 double> &restrict_mother_to_virtual_primary_inv =
2959 (split_face_interpolation_matrices
2960 [dominating_fe_index][cell->active_fe_index()]
2964 double> &restrict_mother_to_virtual_dependent =
2965 (split_face_interpolation_matrices
2966 [dominating_fe_index][cell->active_fe_index()]
2971 constraint_matrix.reinit(
2972 cell->get_fe().n_dofs_per_face(face) -
2975 restrict_mother_to_virtual_dependent.
mmult(
2977 restrict_mother_to_virtual_primary_inv);
2981 scratch_dofs.resize(
2982 cell->get_fe().n_dofs_per_face(face));
2983 cell->face(face)->get_dof_indices(
2984 scratch_dofs, cell->active_fe_index());
2987 primary_dofs.clear();
2988 dependent_dofs.clear();
2989 for (
unsigned int i = 0;
2990 i < cell->get_fe().n_dofs_per_face(face);
2992 if ((*primary_dof_masks[dominating_fe_index]
2993 [cell->active_fe_index()])
2995 primary_dofs.push_back(scratch_dofs[i]);
2997 dependent_dofs.push_back(scratch_dofs[i]);
3003 dependent_dofs.size(),
3004 cell->get_fe().n_dofs_per_face(face) -
3007 filter_constraints(primary_dofs,
3016 neighbor->get_fe().n_dofs_per_face(face),
3019 ensure_existence_of_face_matrix(
3022 face_interpolation_matrices
3023 [dominating_fe_index]
3024 [neighbor->active_fe_index()]);
3027 &restrict_secondface_to_virtual =
3028 *(face_interpolation_matrices
3029 [dominating_fe_index]
3030 [neighbor->active_fe_index()]);
3032 constraint_matrix.reinit(
3033 neighbor->get_fe().n_dofs_per_face(face),
3036 restrict_secondface_to_virtual.
mmult(
3038 restrict_mother_to_virtual_primary_inv);
3040 dependent_dofs.resize(
3041 neighbor->get_fe().n_dofs_per_face(face));
3042 cell->face(face)->get_dof_indices(
3043 dependent_dofs, neighbor->active_fe_index());
3045 filter_constraints(primary_dofs,
3071 template <
int dim,
int spacedim,
typename number>
3078 "The given DoFHandler does not have any DoFs. Did you forget to "
3079 "call dof_handler.distribute_dofs()?"));
3089 dof_handler, constraints, std::integral_constant<int, dim>());
3094 dof_handler, constraints, std::integral_constant<int, dim>());
3101 template <
typename FaceIterator,
typename number>
3104 const FaceIterator &face_1,
3109 const unsigned char combined_orientation,
3110 const number periodicity_factor,
3111 const unsigned int level)
3113 static const int dim = FaceIterator::AccessorType::dimension;
3114 static const int spacedim = FaceIterator::AccessorType::space_dimension;
3129 face_1->get_fe(face_1->nth_active_fe_index(0)).n_unique_faces(), 1);
3131 face_2->get_fe(face_2->nth_active_fe_index(0)).n_unique_faces(), 1);
3132 const unsigned int face_no = 0;
3140 if ((!use_mg) && face_2->has_children())
3142 Assert(face_2->n_children() ==
3146 const unsigned int dofs_per_face =
3147 face_1->get_fe(face_1->nth_active_fe_index(0))
3148 .n_dofs_per_face(face_no);
3152 std::vector<types::global_dof_index> dofs_1(dofs_per_face);
3153 face_1->get_dof_indices(dofs_1, face_1->nth_active_fe_index(0));
3154 for (
unsigned int i = 0; i < dofs_per_face; ++i)
3164 for (
unsigned int c = 0; c < face_2->n_children(); ++c)
3169 const auto &fe = face_1->get_fe(face_1->nth_active_fe_index(0));
3170 fe.get_subface_interpolation_matrix(fe,
3172 subface_interpolation,
3174 subface_interpolation.
mmult(child_transformation, transformation);
3178 child_transformation,
3181 combined_orientation,
3182 periodicity_factor);
3194 Assert(face_1->get_fe(face_1_index) == face_2->get_fe(face_2_index),
3196 "Matching periodic cells need to use the same finite element"));
3202 "The number of components in the mask has to be either "
3203 "zero or equal to the number of components in the finite "
3208 std::vector<types::global_dof_index> dofs_1(dofs_per_face);
3209 std::vector<types::global_dof_index> dofs_2(dofs_per_face);
3212 face_1->get_mg_dof_indices(
level, dofs_1, face_1_index);
3214 face_1->get_dof_indices(dofs_1, face_1_index);
3217 face_2->get_mg_dof_indices(
level, dofs_2, face_2_index);
3219 face_2->get_dof_indices(dofs_2, face_2_index);
3237 for (
unsigned int i = 0; i < dofs_per_face; ++i)
3246 if (
const auto tria =
dynamic_cast<
3248 &face_1->get_triangulation()))
3249 if (tria->with_artificial_cells() &&
3251 for (
unsigned int i = 0; i < dofs_per_face; ++i)
3275 std::map<unsigned int, unsigned int> cell_to_rotated_face_index;
3278 for (
unsigned int i = 0; i < dofs_per_face; ++i)
3292 boost::container::small_vector<
3293 std::pair<typename AffineConstraints<number>::size_type, number>,
3301 for (
unsigned int i = 0; i < dofs_per_face; ++i)
3322 bool is_identity_constrained =
false;
3324 number constraint_factor = periodicity_factor;
3326 constexpr double eps = 1.e-13;
3327 for (
unsigned int jj = 0; jj < dofs_per_face; ++jj)
3329 const auto entry = transformation(i, jj);
3332 if (is_identity_constrained)
3337 is_identity_constrained =
false;
3340 is_identity_constrained =
true;
3342 constraint_factor = entry * periodicity_factor;
3351 if (!is_identity_constrained)
3358 constraint_entries.clear();
3359 constraint_entries.reserve(dofs_per_face);
3361 for (
unsigned int jj = 0; jj < dofs_per_face; ++jj)
3365 const unsigned int j =
3367 jj, 0, combined_orientation)];
3369 if (
std::abs(transformation(i, jj)) > eps)
3370 constraint_entries.emplace_back(dofs_1[j],
3371 transformation(i, jj));
3388 const unsigned int j =
3390 target, 0, combined_orientation)];
3392 auto dof_left = dofs_1[j];
3393 auto dof_right = dofs_2[i];
3399 (dof_left < dof_right &&
3402 std::swap(dof_left, dof_right);
3403 constraint_factor = 1. / constraint_factor;
3419 bool constraints_are_cyclic =
true;
3420 number cycle_constraint_factor = constraint_factor;
3422 for (
auto test_dof = dof_right; test_dof != dof_left;)
3426 constraints_are_cyclic =
false;
3430 const auto &constraint_entries =
3432 if (constraint_entries.size() == 1)
3434 test_dof = constraint_entries[0].first;
3435 cycle_constraint_factor *= constraint_entries[0].second;
3439 constraints_are_cyclic =
false;
3454 if (constraints_are_cyclic)
3456 if (
std::abs(cycle_constraint_factor - number(1.)) > eps)
3462 dof_left, {{dof_right, constraint_factor}}, 0.);
3484 ExcMessage(
"The periodicity constraint is too large. "
3485 "The parameter periodicity_factor might "
3486 "be too large or too small."));
3500 template <
int dim,
int spacedim>
3502 compute_transformation(
3505 const std::vector<unsigned int> &first_vector_components)
3510 const unsigned int face_no = 0;
3516 if (matrix.m() == n_dofs_per_face)
3523 if (first_vector_components.empty() &&
matrix.m() == 0)
3540 using DoFTuple = std::array<unsigned int, spacedim>;
3545 for (
unsigned int i = 0; i < n_dofs_per_face; ++i)
3547 std::vector<unsigned int>::const_iterator comp_it =
3548 std::find(first_vector_components.begin(),
3549 first_vector_components.end(),
3551 if (comp_it != first_vector_components.end())
3553 const unsigned int first_vector_component = *comp_it;
3556 DoFTuple vector_dofs;
3558 unsigned int n_found = 1;
3563 "Error: the finite element does not have enough components "
3564 "to define rotated periodic boundaries."));
3566 for (
unsigned int k = 0; k < n_dofs_per_face; ++k)
3567 if ((k != i) && (quadrature.point(k) == quadrature.point(i)) &&
3569 first_vector_component) &&
3571 first_vector_component + spacedim))
3575 first_vector_component] = k;
3583 for (
unsigned int i = 0; i < spacedim; ++i)
3585 transformation[vector_dofs[i]][vector_dofs[i]] = 0.;
3586 for (
unsigned int j = 0; j < spacedim; ++j)
3587 transformation[vector_dofs[i]][vector_dofs[j]] =
3592 return transformation;
3600 template <
typename FaceIterator,
typename number>
3603 const FaceIterator &face_1,
3607 const unsigned char combined_orientation,
3609 const std::vector<unsigned int> &first_vector_components,
3610 const number periodicity_factor)
3612 static const int dim = FaceIterator::AccessorType::dimension;
3613 static const int spacedim = FaceIterator::AccessorType::space_dimension;
3616 const auto [orientation, rotation, flip] =
3620 (orientation ==
true && flip ==
false && rotation ==
false),
3621 ExcMessage(
"The supplied orientation (orientation, rotation, flip) "
3622 "is invalid for 1d"));
3624 Assert((dim != 2) || (flip ==
false && rotation ==
false),
3625 ExcMessage(
"The supplied orientation (orientation, rotation, flip) "
3626 "is invalid for 2d"));
3629 ExcMessage(
"face_1 and face_2 are equal! Cannot constrain DoFs "
3630 "on the very same face"));
3632 Assert(face_1->at_boundary() && face_2->at_boundary(),
3633 ExcMessage(
"Faces for periodicity constraints must be on the "
3636 Assert(matrix.m() == matrix.n(),
3637 ExcMessage(
"The supplied (rotation or interpolation) matrix must "
3638 "be a square matrix"));
3640 Assert(first_vector_components.empty() || matrix.m() == spacedim,
3641 ExcMessage(
"first_vector_components is nonempty, so matrix must "
3642 "be a rotation matrix exactly of size spacedim"));
3644 if (!face_1->has_children())
3649 face_1->get_fe(face_1->nth_active_fe_index(0)).n_unique_faces(), 1);
3650 const unsigned int face_no = 0;
3653 const unsigned int n_dofs_per_face =
3654 face_1->get_fe(face_1->nth_active_fe_index(0))
3655 .n_dofs_per_face(face_no);
3657 Assert(matrix.m() == 0 ||
3658 (first_vector_components.empty() &&
3659 matrix.m() == n_dofs_per_face) ||
3660 (!first_vector_components.empty() && matrix.m() == spacedim),
3662 "The matrix must have either size 0 or spacedim "
3663 "(if first_vector_components is nonempty) "
3664 "or the size must be equal to the # of DoFs on the face "
3665 "(if first_vector_components is empty)."));
3668 if (!face_2->has_children())
3673 face_2->get_fe(face_2->nth_active_fe_index(0)).n_unique_faces(), 1);
3674 const unsigned int face_no = 0;
3677 const unsigned int n_dofs_per_face =
3678 face_2->get_fe(face_2->nth_active_fe_index(0))
3679 .n_dofs_per_face(face_no);
3681 Assert(matrix.m() == 0 ||
3682 (first_vector_components.empty() &&
3683 matrix.m() == n_dofs_per_face) ||
3684 (!first_vector_components.empty() && matrix.m() == spacedim),
3686 "The matrix must have either size 0 or spacedim "
3687 "(if first_vector_components is nonempty) "
3688 "or the size must be equal to the # of DoFs on the face "
3689 "(if first_vector_components is empty)."));
3693 if (face_1->has_children() && face_2->has_children())
3698 Assert(face_1->n_children() ==
3704 for (
unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face;
3710 const unsigned int face_no = dim == 2 ? 2 : 4;
3714 const auto reference_cell = ReferenceCells::get_hypercube<dim>();
3715 const unsigned int j =
3716 reference_cell.child_cell_on_face(face_no,
3718 combined_orientation);
3724 combined_orientation,
3726 first_vector_components,
3727 periodicity_factor);
3737 face_1->has_children() ?
3738 face_2->get_fe(face_2->nth_active_fe_index(0)) :
3739 face_1->get_fe(face_1->nth_active_fe_index(0));
3744 const unsigned int face_no = 0;
3750 if (n_dofs_per_face == 0)
3754 compute_transformation(fe, matrix, first_vector_components);
3756 if (!face_2->has_children())
3760 if (first_vector_components.empty() && matrix.m() == 0)
3767 combined_orientation,
3768 periodicity_factor);
3773 inverse.
invert(transformation);
3780 combined_orientation,
3781 periodicity_factor);
3795 const auto face_reference_cell = face_1->reference_cell();
3802 face_reference_cell.get_inverse_combined_orientation(
3803 combined_orientation),
3804 periodicity_factor);
3811 template <
int dim,
int spacedim,
typename number>
3818 const std::vector<unsigned int> &first_vector_components,
3819 const number periodicity_factor)
3822 for (
auto &pair : periodic_faces)
3825 const FaceIterator face_1 = pair.cell[0]->face(pair.face_idx[0]);
3826 const FaceIterator face_2 = pair.cell[1]->face(pair.face_idx[1]);
3828 Assert(face_1->at_boundary() && face_2->at_boundary(),
3841 first_vector_components,
3842 periodicity_factor);
3850 template <
int dim,
int spacedim,
typename number>
3855 const unsigned int direction,
3858 const number periodicity_factor)
3863 ExcMessage(
"The boundary indicators b_id1 and b_id2 must be "
3864 "different to denote different boundaries."));
3872 dof_handler, b_id1, b_id2, direction, matched_faces);
3874 make_periodicity_constraints<dim, spacedim>(matched_faces,
3877 std::vector<unsigned int>(),
3878 periodicity_factor);
3883 template <
int dim,
int spacedim,
typename number>
3887 const unsigned int direction,
3890 const number periodicity_factor)
3906 make_periodicity_constraints<dim, spacedim>(matched_faces,
3909 std::vector<unsigned int>(),
3910 periodicity_factor);
3924 template <
int dim,
int spacedim>
3929#ifdef DEAL_II_WITH_MPI
3930 std::vector<::LinearAlgebra::distributed::Vector<double>>
3945 template <
int dim,
int spacedim>
3947 compute_intergrid_weights_3(
3950 const unsigned int coarse_component,
4009 for (
unsigned int local_dof = 0; local_dof < copy_data.
dofs_per_cell;
4015 const unsigned int local_parameter_dof =
4024 coarse_to_fine_grid_map[cell]->set_dof_values_by_interpolation(
4025 parameter_dofs[local_parameter_dof],
4037 template <
int dim,
int spacedim>
4039 copy_intergrid_weights_3(
4040 const Assembler::CopyData<dim, spacedim> ©_data,
4041 const unsigned int coarse_component,
4043 const std::vector<types::global_dof_index> &weight_mapping,
4044 const bool is_called_in_parallel,
4045 std::vector<std::map<types::global_dof_index, float>> &weights)
4047 unsigned int pos = 0;
4048 for (
unsigned int local_dof = 0; local_dof < copy_data.dofs_per_cell;
4075 i < copy_data.global_parameter_representation[pos].size();
4081 if (copy_data.global_parameter_representation[pos](i) != 0)
4084 wi = copy_data.parameter_dof_indices[local_dof],
4085 wj = weight_mapping[i];
4087 copy_data.global_parameter_representation[pos](i);
4090 else if (!is_called_in_parallel)
4095 Assert(copy_data.global_parameter_representation[pos](i) ==
4111 template <
int dim,
int spacedim>
4113 compute_intergrid_weights_2(
4115 const unsigned int coarse_component,
4118 const std::vector<types::global_dof_index> &weight_mapping,
4119 std::vector<std::map<types::global_dof_index, float>> &weights)
4121 Assembler::CopyData<dim, spacedim> copy_data;
4123 unsigned int n_interesting_dofs = 0;
4124 for (
unsigned int local_dof = 0;
4125 local_dof < coarse_grid.
get_fe().n_dofs_per_cell();
4129 ++n_interesting_dofs;
4131 copy_data.global_parameter_representation.resize(n_interesting_dofs);
4133 bool is_called_in_parallel =
false;
4134 for (std::size_t i = 0;
4135 i < copy_data.global_parameter_representation.size();
4138#ifdef DEAL_II_WITH_MPI
4139 MPI_Comm communicator = MPI_COMM_SELF;
4142 const typename ::parallel::TriangulationBase<dim,
4144 &tria =
dynamic_cast<const typename ::parallel::
4145 TriangulationBase<dim, spacedim> &
>(
4146 coarse_to_fine_grid_map.get_destination_grid()
4147 .get_triangulation());
4148 communicator = tria.get_mpi_communicator();
4149 is_called_in_parallel =
true;
4151 catch (std::bad_cast &)
4157 const IndexSet locally_relevant_dofs =
4159 coarse_to_fine_grid_map.get_destination_grid());
4161 copy_data.global_parameter_representation[i].reinit(
4162 coarse_to_fine_grid_map.get_destination_grid()
4163 .locally_owned_dofs(),
4164 locally_relevant_dofs,
4168 copy_data.global_parameter_representation[i].reinit(n_fine_dofs);
4175 &coarse_to_fine_grid_map,
4179 const Assembler::Scratch &,
4180 Assembler::CopyData<dim, spacedim> ©_data) {
4181 compute_intergrid_weights_3<dim, spacedim>(cell,
4184 coarse_grid.get_fe(),
4185 coarse_to_fine_grid_map,
4193 is_called_in_parallel,
4194 &weights](
const Assembler::CopyData<dim, spacedim> ©_data) {
4195 copy_intergrid_weights_3<dim, spacedim>(copy_data,
4197 coarse_grid.get_fe(),
4199 is_called_in_parallel,
4207 Assembler::Scratch(),
4210#ifdef DEAL_II_WITH_MPI
4211 for (std::size_t i = 0;
4212 i < copy_data.global_parameter_representation.size();
4214 copy_data.global_parameter_representation[i].update_ghost_values();
4225 template <
int dim,
int spacedim>
4227 compute_intergrid_weights_1(
4229 const unsigned int coarse_component,
4231 const unsigned int fine_component,
4233 std::vector<std::map<types::global_dof_index, float>> &weights,
4234 std::vector<types::global_dof_index> &weight_mapping)
4238 &fine_fe = fine_grid.
get_fe();
4242 n_fine_dofs = fine_grid.
n_dofs();
4245 const unsigned int fine_dofs_per_cell = fine_fe.n_dofs_per_cell();
4249 const unsigned int coarse_dofs_per_cell_component =
4263 Assert(&coarse_to_fine_grid_map.get_source_grid() == &coarse_grid,
4265 Assert(&coarse_to_fine_grid_map.get_destination_grid() == &fine_grid,
4277 fine_fe.component_to_base_index(fine_component).first),
4283 for (
const auto &cell : coarse_grid.active_cell_iterators())
4309 std::vector<::Vector<double>> parameter_dofs(
4310 coarse_dofs_per_cell_component,
4315 for (
unsigned int local_coarse_dof = 0;
4316 local_coarse_dof < coarse_dofs_per_cell_component;
4318 for (
unsigned int fine_dof = 0; fine_dof < fine_fe.n_dofs_per_cell();
4320 if (fine_fe.system_to_component_index(fine_dof) ==
4321 std::make_pair(fine_component, local_coarse_dof))
4323 parameter_dofs[local_coarse_dof](fine_dof) = 1.;
4330 unsigned int n_parameters_on_fine_grid = 0;
4334 std::vector<bool> dof_is_interesting(fine_grid.
n_dofs(),
false);
4335 std::vector<types::global_dof_index> local_dof_indices(
4336 fine_fe.n_dofs_per_cell());
4338 for (
const auto &cell : fine_grid.active_cell_iterators() |
4341 cell->get_dof_indices(local_dof_indices);
4342 for (
unsigned int i = 0; i < fine_fe.n_dofs_per_cell(); ++i)
4343 if (fine_fe.system_to_component_index(i).first ==
4345 dof_is_interesting[local_dof_indices[i]] =
true;
4348 n_parameters_on_fine_grid = std::count(dof_is_interesting.begin(),
4349 dof_is_interesting.end(),
4356 weights.resize(n_coarse_dofs);
4358 weight_mapping.clear();
4362 std::vector<types::global_dof_index> local_dof_indices(
4363 fine_fe.n_dofs_per_cell());
4364 unsigned int next_free_index = 0;
4365 for (
const auto &cell : fine_grid.active_cell_iterators() |
4368 cell->get_dof_indices(local_dof_indices);
4369 for (
unsigned int i = 0; i < fine_fe.n_dofs_per_cell(); ++i)
4372 if ((fine_fe.system_to_component_index(i).first ==
4374 (weight_mapping[local_dof_indices[i]] ==
4377 weight_mapping[local_dof_indices[i]] = next_free_index;
4382 Assert(next_free_index == n_parameters_on_fine_grid,
4394 compute_intergrid_weights_2(coarse_grid,
4396 coarse_to_fine_grid_map,
4416 for (
unsigned int col = 0; col < n_parameters_on_fine_grid; ++col)
4420 if (weights[row].find(col) != weights[row].
end())
4421 sum += weights[row][col];
4422 Assert((std::fabs(sum - 1) < 1.e-12) ||
4429 return n_parameters_on_fine_grid;
4438 template <
int dim,
int spacedim>
4442 const unsigned int coarse_component,
4444 const unsigned int fine_component,
4450 ExcMessage(
"This function is not yet implemented for DoFHandlers "
4451 "using hp-capabilities."));
4472 std::vector<std::map<types::global_dof_index, float>> weights;
4478 std::vector<types::global_dof_index> weight_mapping;
4480 const unsigned int n_parameters_on_fine_grid =
4481 internal::compute_intergrid_weights_1(coarse_grid,
4485 coarse_to_fine_grid_map,
4488 (void)n_parameters_on_fine_grid;
4492 n_fine_dofs = fine_grid.
n_dofs();
4500 mask[coarse_component] =
true;
4502 coarse_dof_is_parameter =
4503 extract_dofs<dim, spacedim>(coarse_grid,
ComponentMask(mask));
4515 std::vector<types::global_dof_index> representants(
4518 parameter_dof < n_coarse_dofs;
4520 if (coarse_dof_is_parameter.
is_element(parameter_dof))
4527 std::map<types::global_dof_index, float>::const_iterator i =
4528 weights[parameter_dof].begin();
4529 for (; i != weights[parameter_dof].end(); ++i)
4539 for (; global_dof < weight_mapping.size(); ++global_dof)
4540 if (weight_mapping[global_dof] ==
4546 representants[parameter_dof] = global_dof;
4562 std::vector<std::pair<types::global_dof_index, double>> constraint_line;
4581 std::map<types::global_dof_index, float>::const_iterator col_entry =
4583 for (; first_used_row < n_coarse_dofs; ++first_used_row)
4585 col_entry = weights[first_used_row].find(col);
4586 if (col_entry != weights[first_used_row].end())
4590 Assert(col_entry != weights[first_used_row].end(),
4593 if ((col_entry->second == 1) &&
4594 (representants[first_used_row] == global_dof))
4602 constraint_line.clear();
4604 row < n_coarse_dofs;
4607 const std::map<types::global_dof_index, float>::const_iterator j =
4608 weights[row].find(col);
4609 if ((j != weights[row].end()) && (j->second != 0))
4610 constraint_line.emplace_back(representants[row], j->second);
4619 template <
int dim,
int spacedim>
4623 const unsigned int coarse_component,
4625 const unsigned int fine_component,
4627 std::vector<std::map<types::global_dof_index, float>>
4628 &transfer_representation)
4632 ExcMessage(
"This function is not yet implemented for DoFHandlers "
4633 "using hp-capabilities."));
4654 std::vector<std::map<types::global_dof_index, float>> weights;
4660 std::vector<types::global_dof_index> weight_mapping;
4662 internal::compute_intergrid_weights_1(coarse_grid,
4666 coarse_to_fine_grid_map,
4672 std::count_if(weight_mapping.begin(),
4673 weight_mapping.end(),
4675 return dof != numbers::invalid_dof_index;
4679 std::vector<types::global_dof_index> inverse_weight_mapping(
4688 Assert((inverse_weight_mapping[parameter_dof] ==
4692 inverse_weight_mapping[parameter_dof] = i;
4699 transfer_representation.clear();
4700 transfer_representation.resize(n_rows);
4705 std::map<types::global_dof_index, float>::const_iterator j =
4707 for (; j != weights[i].end(); ++j)
4712 transfer_representation[p][i] = j->second;
4719 template <
int dim,
int spacedim,
typename number>
4728 ExcMessage(
"The number of components in the mask has to be either "
4729 "zero or equal to the number of components in the finite "
4738 std::vector<types::global_dof_index> face_dofs;
4741 std::vector<types::global_dof_index> cell_dofs;
4749 std::set<types::global_dof_index> dofs_already_treated;
4752 if (!cell->is_artificial() && cell->at_boundary())
4758 cell->get_dof_indices(cell_dofs);
4760 for (
const auto face_no : cell->face_indices())
4763 cell->face(face_no);
4767 if (face->at_boundary() &&
4769 (face->boundary_id() == boundary_id)))
4773 face->get_dof_indices(face_dofs, cell->active_fe_index());
4778 if (dofs_already_treated.find(face_dof) ==
4779 dofs_already_treated.end())
4783 const std::vector<types::global_dof_index>::iterator
4784 it_index_on_cell = std::find(cell_dofs.begin(),
4787 Assert(it_index_on_cell != cell_dofs.end(),
4789 const unsigned int index_on_cell =
4790 std::distance(cell_dofs.begin(), it_index_on_cell);
4792 cell->get_fe().get_nonzero_components(index_on_cell);
4795 for (
unsigned int c = 0; c < n_components; ++c)
4796 if (nonzero_component_array[c] && component_mask[c])
4812 Assert(zero_boundary_constraints
4813 .is_inhomogeneously_constrained(
4820 dofs_already_treated.insert(face_dof);
4829 template <
int dim,
int spacedim,
typename number>
4838 zero_boundary_constraints,
4849#include "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
virtual unsigned int face_to_cell_index(const unsigned int face_dof_index, const unsigned int face, const unsigned char combined_orientation=ReferenceCell::default_combined_face_orientation()) 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 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
#define DEAL_II_NAMESPACE_CLOSE
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={})
#define DEAL_II_ASSERT_UNREACHABLE()
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)
VectorType::value_type * end(VectorType &V)
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 unsigned char combined_face_orientation)
const types::boundary_id invalid_boundary_id
const types::fe_index invalid_fe_index
static const unsigned int invalid_unsigned_int
const types::global_dof_index invalid_dof_index
typename type_identity< T >::type type_identity_t
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
unsigned short int fe_index
static unsigned int n_children(const RefinementCase< dim > &refinement_case)