50 zero_indices(
unsigned int (&indices)[dim])
52 for (
unsigned int d = 0;
d < dim; ++
d)
63 increment_indices(
unsigned int (&indices)[dim],
const unsigned int dofs1d)
66 for (
int d = 0;
d < dim - 1; ++
d)
67 if (indices[
d] == dofs1d)
83 template <
class PolynomialType,
int xdim,
int xspacedim>
90 template <
int spacedim>
99 template <
int spacedim>
104 const unsigned int dim = 2;
106 unsigned int q_deg = fe.
degree;
155 std::vector<
Point<dim - 1>> constraint_points;
157 constraint_points.emplace_back(0.5);
161 const unsigned int n = q_deg - 1;
162 const double step = 1. / q_deg;
164 for (
unsigned int i = 1; i <= n; ++i)
165 constraint_points.push_back(
169 for (
unsigned int i = 1; i <= n; ++i)
170 constraint_points.push_back(
183 const std::vector<unsigned int> &index_map_inverse =
185 const std::vector<unsigned int> face_index_map =
192 for (
unsigned int i = 0; i < constraint_points.size(); ++i)
193 for (
unsigned int j = 0; j < q_deg + 1; ++j)
196 p[0] = constraint_points[i](0);
198 fe.
poly_space.compute_value(index_map_inverse[j], p);
210 template <
int spacedim>
215 const unsigned int dim = 3;
217 unsigned int q_deg = fe.
degree;
237 std::vector<
Point<dim - 1>> constraint_points;
240 constraint_points.emplace_back(0.5, 0.5);
243 constraint_points.emplace_back(0, 0.5);
244 constraint_points.emplace_back(1, 0.5);
245 constraint_points.emplace_back(0.5, 0);
246 constraint_points.emplace_back(0.5, 1);
250 const unsigned int n = q_deg - 1;
251 const double step = 1. / q_deg;
252 std::vector<
Point<dim - 2>> line_support_points(n);
253 for (
unsigned int i = 0; i < n; ++i)
254 line_support_points[i](0) = (i + 1) * step;
255 Quadrature<dim - 2> qline(line_support_points);
258 std::vector<
Point<dim - 1>> p_line(n);
264 for (
unsigned int i = 0; i < n; ++i)
268 for (
unsigned int i = 0; i < n; ++i)
272 for (
unsigned int i = 0; i < n; ++i)
276 for (
unsigned int i = 0; i < n; ++i)
280 for (
unsigned int face = 0;
283 for (
unsigned int subface = 0;
291 constraint_points.insert(constraint_points.end(),
297 std::vector<
Point<dim - 1>> inner_points(n * n);
298 for (
unsigned int i = 0, iy = 1; iy <= n; ++iy)
299 for (
unsigned int ix = 1; ix <= n; ++ix)
303 for (
unsigned int child = 0;
306 for (
const auto &inner_point : inner_points)
307 constraint_points.push_back(
314 const unsigned int pnts = (q_deg + 1) * (q_deg + 1);
318 const std::vector<unsigned int> &index_map_inverse =
320 const std::vector<unsigned int> face_index_map =
330 for (
unsigned int i = 0; i < constraint_points.size(); ++i)
332 const double interval =
static_cast<double>(q_deg * 2);
333 bool mirror[dim - 1];
345 for (
unsigned int k = 0; k < dim - 1; ++k)
347 const int coord_int =
348 static_cast<int>(constraint_points[i](k) * interval + 0.25);
349 constraint_point(k) = 1. * coord_int / interval;
371 mirror[k] = (constraint_point(k) > 0.5);
373 constraint_point(k) = 1.0 - constraint_point(k);
376 for (
unsigned int j = 0; j < pnts; ++j)
378 unsigned int indices[2] = {j % (q_deg + 1), j / (q_deg + 1)};
380 for (
unsigned int k = 0; k < 2; ++k)
382 indices[k] = q_deg - indices[k];
384 const unsigned int new_index =
385 indices[1] * (q_deg + 1) + indices[0];
388 fe.
poly_space.compute_value(index_map_inverse[new_index],
405 template <
class PolynomialType,
int dim,
int spacedim>
409 const std::vector<bool> & restriction_is_additive_flags)
413 restriction_is_additive_flags,
423 template <
class PolynomialType,
int dim,
int spacedim>
426 const std::vector<
Point<1>> &points)
429 ExcMessage(
"The first support point has to be zero."));
430 Assert(points.back()[0] == 1,
431 ExcMessage(
"The last support point has to be one."));
436 const unsigned int q_dofs_per_cell =
437 Utilities::fixed_power<dim>(q_degree + 1);
438 Assert(q_dofs_per_cell == this->dofs_per_cell ||
439 q_dofs_per_cell + 1 == this->dofs_per_cell ||
440 q_dofs_per_cell + dim == this->dofs_per_cell,
444 std::vector<unsigned int> renumber =
445 FETools::hierarchic_to_lexicographic_numbering<dim>(q_degree);
446 for (
unsigned int i = q_dofs_per_cell; i < this->dofs_per_cell; ++i)
447 renumber.push_back(i);
448 this->poly_space.set_numbering(renumber);
469 template <
class PolynomialType,
int dim,
int spacedim>
481 Assert(interpolation_matrix.
m() == this->dofs_per_cell,
483 this->dofs_per_cell));
489 const unsigned int q_dofs_per_cell =
490 Utilities::fixed_power<dim>(q_degree + 1);
491 const unsigned int source_q_dofs_per_cell =
492 Utilities::fixed_power<dim>(source_fe->degree + 1);
498 for (
unsigned int j = 0; j < q_dofs_per_cell; ++j)
501 const Point<dim> p = this->unit_support_points[j];
505 Assert(std::abs(this->poly_space.compute_value(j, p) - 1.) < 1
e-13,
508 for (
unsigned int i = 0; i < source_q_dofs_per_cell; ++i)
509 interpolation_matrix(j, i) =
510 source_fe->poly_space.compute_value(i, p);
514 if (q_dofs_per_cell < this->dofs_per_cell)
517 for (
unsigned int i = 0; i < source_q_dofs_per_cell; ++i)
518 interpolation_matrix(q_dofs_per_cell, i) = 0.;
519 for (
unsigned int j = 0; j < q_dofs_per_cell; ++j)
520 interpolation_matrix(j, source_q_dofs_per_cell) = 0.;
521 interpolation_matrix(q_dofs_per_cell, source_q_dofs_per_cell) = 1.;
525 const double eps = 2
e-13 * q_degree * dim;
526 for (
unsigned int i = 0; i < this->dofs_per_cell; ++i)
527 for (
unsigned int j = 0; j < source_fe->dofs_per_cell; ++j)
529 interpolation_matrix(i, j) = 0.;
533 for (
unsigned int i = 0; i < this->dofs_per_cell; ++i)
536 for (
unsigned int j = 0; j < source_fe->dofs_per_cell; ++j)
537 sum += interpolation_matrix(i, j);
562 spacedim>::ExcInterpolationNotImplemented()));
567 template <
class PolynomialType,
int dim,
int spacedim>
574 get_subface_interpolation_matrix(source_fe,
576 interpolation_matrix);
581 template <
class PolynomialType,
int dim,
int spacedim>
585 const unsigned int subface,
599 Assert(interpolation_matrix.
n() == this->dofs_per_face,
601 this->dofs_per_face));
609 this->dofs_per_face <= source_fe->dofs_per_face,
611 spacedim>::ExcInterpolationNotImplemented()));
615 source_fe->get_unit_face_support_points());
620 double eps = 2
e-13 * q_degree * (dim - 1);
631 for (
unsigned int i = 0; i < source_fe->dofs_per_face; ++i)
633 const Point<dim> &p = subface_quadrature.point(i);
635 for (
unsigned int j = 0; j < this->dofs_per_face; ++j)
637 double matrix_entry =
638 this->shape_value(this->face_to_cell_index(j, 0), p);
648 interpolation_matrix(i, j) = matrix_entry;
654 for (
unsigned int j = 0; j < source_fe->dofs_per_face; ++j)
658 for (
unsigned int i = 0; i < this->dofs_per_face; ++i)
659 sum += interpolation_matrix(j, i);
664 else if (
dynamic_cast<const FE_Nothing<dim> *
>(&x_source_fe) !=
nullptr)
672 spacedim>::ExcInterpolationNotImplemented()));
677 template <
class PolynomialType,
int dim,
int spacedim>
686 template <
class PolynomialType,
int dim,
int spacedim>
687 std::vector<std::pair<unsigned int, unsigned int>>
696 &fe_other) !=
nullptr)
698 return std::vector<std::pair<unsigned int, unsigned int>>(
699 1, std::make_pair(0
U, 0
U));
705 return std::vector<std::pair<unsigned int, unsigned int>>();
716 return std::vector<std::pair<unsigned int, unsigned int>>();
721 return std::vector<std::pair<unsigned int, unsigned int>>();
727 template <
class PolynomialType,
int dim,
int spacedim>
728 std::vector<std::pair<unsigned int, unsigned int>>
746 const unsigned int p = this->degree;
747 const unsigned int q = fe_q_other->degree;
749 std::vector<std::pair<unsigned int, unsigned int>> identities;
751 const std::vector<unsigned int> &index_map_inverse =
752 this->poly_space.get_numbering_inverse();
753 const std::vector<unsigned int> &index_map_inverse_other =
754 fe_q_other->poly_space.get_numbering_inverse();
756 for (
unsigned int i = 0; i < p - 1; ++i)
757 for (
unsigned int j = 0; j < q - 1; ++j)
759 this->unit_support_points[index_map_inverse[i + 1]][0] -
760 fe_q_other->unit_support_points[index_map_inverse_other[j + 1]]
762 identities.emplace_back(i, j);
770 return std::vector<std::pair<unsigned int, unsigned int>>();
781 return std::vector<std::pair<unsigned int, unsigned int>>();
786 return std::vector<std::pair<unsigned int, unsigned int>>();
792 template <
class PolynomialType,
int dim,
int spacedim>
793 std::vector<std::pair<unsigned int, unsigned int>>
808 const unsigned int p = this->degree;
809 const unsigned int q = fe_q_other->degree;
811 std::vector<std::pair<unsigned int, unsigned int>> identities;
813 const std::vector<unsigned int> &index_map_inverse =
814 this->poly_space.get_numbering_inverse();
815 const std::vector<unsigned int> &index_map_inverse_other =
816 fe_q_other->poly_space.get_numbering_inverse();
818 for (
unsigned int i1 = 0; i1 < p - 1; ++i1)
819 for (
unsigned int i2 = 0; i2 < p - 1; ++i2)
820 for (
unsigned int j1 = 0; j1 < q - 1; ++j1)
821 for (
unsigned int j2 = 0; j2 < q - 1; ++j2)
823 this->unit_support_points[index_map_inverse[i1 + 1]][0] -
825 ->unit_support_points[index_map_inverse_other[j1 + 1]]
828 this->unit_support_points[index_map_inverse[i2 + 1]][0] -
830 ->unit_support_points[index_map_inverse_other[j2 + 1]]
832 identities.emplace_back(i1 * (p - 1) + i2, j1 * (q - 1) + j2);
840 return std::vector<std::pair<unsigned int, unsigned int>>();
851 return std::vector<std::pair<unsigned int, unsigned int>>();
856 return std::vector<std::pair<unsigned int, unsigned int>>();
868 template <
class PolynomialType,
int dim,
int spacedim>
871 const std::vector<
Point<1>> &points)
873 const std::vector<unsigned int> &index_map_inverse =
874 this->poly_space.get_numbering_inverse();
885 this->unit_support_points.resize(support_quadrature.
size());
886 for (
unsigned int k = 0; k < support_quadrature.
size(); ++k)
887 this->unit_support_points[index_map_inverse[k]] =
888 support_quadrature.
point(k);
893 template <
class PolynomialType,
int dim,
int spacedim>
896 const std::vector<
Point<1>> &points)
902 this->unit_face_support_points.resize(
903 Utilities::fixed_power<dim - 1>(q_degree + 1));
906 const std::vector<unsigned int> face_index_map =
914 const Quadrature<dim - 1> support_quadrature(support_1d);
918 this->unit_face_support_points.resize(support_quadrature.size());
919 for (
unsigned int k = 0; k < support_quadrature.size(); ++k)
920 this->unit_face_support_points[face_index_map[k]] =
921 support_quadrature.point(k);
926 template <
class PolynomialType,
int dim,
int spacedim>
935 Assert(this->adjust_quad_dof_index_for_face_orientation_table.n_elements() ==
936 8 * this->dofs_per_quad,
939 const unsigned int n = q_degree - 1;
960 for (
unsigned int local = 0; local < this->dofs_per_quad; ++local)
964 unsigned int i = local % n, j = local / n;
967 this->adjust_quad_dof_index_for_face_orientation_table(local, 0) =
970 this->adjust_quad_dof_index_for_face_orientation_table(local, 1) =
971 i + (n - 1 - j) * n - local;
973 this->adjust_quad_dof_index_for_face_orientation_table(local, 2) =
974 (n - 1 - j) + (n - 1 - i) * n - local;
976 this->adjust_quad_dof_index_for_face_orientation_table(local, 3) =
977 (n - 1 - i) + j * n - local;
979 this->adjust_quad_dof_index_for_face_orientation_table(local, 4) = 0;
981 this->adjust_quad_dof_index_for_face_orientation_table(local, 5) =
982 j + (n - 1 - i) * n - local;
984 this->adjust_quad_dof_index_for_face_orientation_table(local, 6) =
985 (n - 1 - i) + (n - 1 - j) * n - local;
987 this->adjust_quad_dof_index_for_face_orientation_table(local, 7) =
988 (n - 1 - j) + i * n - local;
992 for (
unsigned int i = 0; i < this->dofs_per_line; ++i)
993 this->adjust_line_dof_index_for_line_orientation_table[i] =
994 this->dofs_per_line - 1 - i - i;
999 template <
class PolynomialType,
int dim,
int spacedim>
1002 const unsigned int face_index,
1003 const unsigned int face,
1004 const bool face_orientation,
1005 const bool face_flip,
1006 const bool face_rotation)
const
1022 if (face_index < this->first_face_line_index)
1027 const unsigned int face_vertex = face_index / this->dofs_per_vertex;
1028 const unsigned int dof_index_on_vertex =
1029 face_index % this->dofs_per_vertex;
1034 face, face_vertex, face_orientation, face_flip, face_rotation) *
1035 this->dofs_per_vertex +
1036 dof_index_on_vertex);
1038 else if (face_index < this->first_face_quad_index)
1043 const unsigned int index = face_index - this->first_face_line_index;
1045 const unsigned int face_line = index / this->dofs_per_line;
1046 const unsigned int dof_index_on_line = index % this->dofs_per_line;
1050 unsigned int adjusted_dof_index_on_line = 0;
1060 if (face_flip ==
false)
1061 adjusted_dof_index_on_line = dof_index_on_line;
1063 adjusted_dof_index_on_line =
1064 this->dofs_per_line - dof_index_on_line - 1;
1075 Assert((this->dofs_per_line <= 1) ||
1076 ((face_orientation ==
true) && (face_flip ==
false) &&
1077 (face_rotation ==
false)),
1079 adjusted_dof_index_on_line = dof_index_on_line;
1086 return (this->first_line_index +
1088 face, face_line, face_orientation, face_flip, face_rotation) *
1089 this->dofs_per_line +
1090 adjusted_dof_index_on_line);
1098 const unsigned int index = face_index - this->first_face_quad_index;
1103 Assert((this->dofs_per_quad <= 1) ||
1104 ((face_orientation ==
true) && (face_flip ==
false) &&
1105 (face_rotation ==
false)),
1107 return (this->first_quad_index + face * this->dofs_per_quad + index);
1113 template <
class PolynomialType,
int dim,
int spacedim>
1114 std::vector<unsigned int>
1116 const unsigned int degree)
1119 AssertThrow(degree > 0,
typename FEQ::ExcFEQCannotHaveDegree0());
1120 std::vector<unsigned int> dpo(dim + 1, 1
U);
1121 for (
unsigned int i = 1; i < dpo.size(); ++i)
1122 dpo[i] = dpo[i - 1] * (degree - 1);
1128 template <
class PolynomialType,
int dim,
int spacedim>
1131 const std::vector<
Point<1>> &points)
1133 Implementation::initialize_constraints(points, *
this);
1138 template <
class PolynomialType,
int dim,
int spacedim>
1141 const unsigned int child,
1148 "Prolongation matrices are only available for refined cells!"));
1152 if (this->prolongation[refinement_case - 1][child].n() == 0)
1154 std::lock_guard<std::mutex> lock(this->mutex);
1157 if (this->prolongation[refinement_case - 1][child].n() ==
1158 this->dofs_per_cell)
1159 return this->prolongation[refinement_case - 1][child];
1162 const unsigned int q_dofs_per_cell =
1163 Utilities::fixed_power<dim>(q_degree + 1);
1171 const double eps = 1
e-15 * q_degree * dim;
1176 for (
unsigned int i = 0; i < q_dofs_per_cell; ++i)
1179 i, this->unit_support_points[i])) <
eps,
1181 "to one or zero in a nodal point. "
1182 "This typically indicates that the "
1183 "polynomial interpolation is "
1184 "ill-conditioned such that round-off "
1185 "prevents the sum to be one."));
1186 for (
unsigned int j = 0; j < q_dofs_per_cell; ++j)
1189 i, this->unit_support_points[j])) <
eps,
1191 "The Lagrange polynomial does not evaluate "
1192 "to one or zero in a nodal point. "
1193 "This typically indicates that the "
1194 "polynomial interpolation is "
1195 "ill-conditioned such that round-off "
1196 "prevents the sum to be one."));
1204 const unsigned int dofs1d = q_degree + 1;
1205 std::vector<Table<2, double>> subcell_evaluations(
1207 const std::vector<unsigned int> &index_map_inverse =
1208 this->poly_space.get_numbering_inverse();
1212 unsigned int step_size_diag = 0;
1214 unsigned int factor = 1;
1215 for (
unsigned int d = 0;
d < dim; ++
d)
1217 step_size_diag += factor;
1226 for (
unsigned int j = 0; j < dofs1d; ++j)
1228 const unsigned int diag_comp = index_map_inverse[j * step_size_diag];
1229 const Point<dim> p_subcell = this->unit_support_points[diag_comp];
1234 for (
unsigned int i = 0; i < dofs1d; ++i)
1235 for (
unsigned int d = 0;
d < dim; ++
d)
1240 const double cell_value =
1241 this->poly_space.compute_value(index_map_inverse[i],
point);
1261 subcell_evaluations[
d](j, i) = 0;
1263 subcell_evaluations[
d](j, i) = cell_value;
1269 unsigned int j_indices[dim];
1270 internal::FE_Q_Base::zero_indices<dim>(j_indices);
1271 for (
unsigned int j = 0; j < q_dofs_per_cell; j += dofs1d)
1273 unsigned int i_indices[dim];
1274 internal::FE_Q_Base::zero_indices<dim>(i_indices);
1275 for (
unsigned int i = 0; i < q_dofs_per_cell; i += dofs1d)
1277 double val_extra_dim = 1.;
1278 for (
unsigned int d = 1;
d < dim; ++
d)
1280 subcell_evaluations[
d](j_indices[
d - 1], i_indices[
d - 1]);
1284 for (
unsigned int jj = 0; jj < dofs1d; ++jj)
1286 const unsigned int j_ind = index_map_inverse[j + jj];
1287 for (
unsigned int ii = 0; ii < dofs1d; ++ii)
1288 prolongate(j_ind, index_map_inverse[i + ii]) =
1289 val_extra_dim * subcell_evaluations[0](jj, ii);
1294 internal::FE_Q_Base::increment_indices<dim>(i_indices, dofs1d);
1297 internal::FE_Q_Base::increment_indices<dim>(j_indices, dofs1d);
1302 if (q_dofs_per_cell < this->dofs_per_cell)
1303 prolongate(q_dofs_per_cell, q_dofs_per_cell) = 1.;
1308 for (
unsigned int row = 0; row < this->dofs_per_cell; ++row)
1311 for (
unsigned int col = 0; col < this->dofs_per_cell; ++col)
1312 sum += prolongate(row, col);
1316 "prolongation matrix do not add to one. "
1317 "This typically indicates that the "
1318 "polynomial interpolation is "
1319 "ill-conditioned such that round-off "
1320 "prevents the sum to be one."));
1326 this->prolongation[refinement_case - 1][child]));
1330 return this->prolongation[refinement_case - 1][child];
1335 template <
class PolynomialType,
int dim,
int spacedim>
1338 const unsigned int child,
1345 "Restriction matrices are only available for refined cells!"));
1349 if (this->restriction[refinement_case - 1][child].n() == 0)
1351 std::lock_guard<std::mutex> lock(this->mutex);
1354 if (this->restriction[refinement_case - 1][child].n() ==
1355 this->dofs_per_cell)
1356 return this->restriction[refinement_case - 1][child];
1359 this->dofs_per_cell);
1361 const unsigned int q_dofs_per_cell =
1362 Utilities::fixed_power<dim>(q_degree + 1);
1382 const double eps = 1
e-15 * q_degree * dim;
1383 const std::vector<unsigned int> &index_map_inverse =
1384 this->poly_space.get_numbering_inverse();
1386 const unsigned int dofs1d = q_degree + 1;
1387 std::vector<Tensor<1, dim>> evaluations1d(dofs1d);
1389 my_restriction.
reinit(this->dofs_per_cell, this->dofs_per_cell);
1391 for (
unsigned int i = 0; i < q_dofs_per_cell; ++i)
1393 unsigned int mother_dof = index_map_inverse[i];
1394 const Point<dim> p_cell = this->unit_support_points[mother_dof];
1408 for (
unsigned int j = 0; j < dofs1d; ++j)
1409 for (
unsigned int d = 0;
d < dim; ++
d)
1413 evaluations1d[j][
d] =
1414 this->poly_space.compute_value(index_map_inverse[j],
1417 unsigned int j_indices[dim];
1418 internal::FE_Q_Base::zero_indices<dim>(j_indices);
1419 double sum_check = 0;
1420 for (
unsigned int j = 0; j < q_dofs_per_cell; j += dofs1d)
1422 double val_extra_dim = 1.;
1423 for (
unsigned int d = 1;
d < dim; ++
d)
1424 val_extra_dim *= evaluations1d[j_indices[
d - 1]][
d];
1425 for (
unsigned int jj = 0; jj < dofs1d; ++jj)
1433 const double val = val_extra_dim * evaluations1d[jj][0];
1434 const unsigned int child_dof = index_map_inverse[j + jj];
1436 my_restriction(mother_dof, child_dof) = 1.;
1438 my_restriction(mother_dof, child_dof) = val;
1441 internal::FE_Q_Base::increment_indices<dim>(j_indices,
1447 "restriction matrix do not add to one. "
1448 "This typically indicates that the "
1449 "polynomial interpolation is "
1450 "ill-conditioned such that round-off "
1451 "prevents the sum to be one."));
1455 if (q_dofs_per_cell < this->dofs_per_cell)
1456 my_restriction(this->dofs_per_cell - 1, this->dofs_per_cell - 1) =
1464 this->restriction[refinement_case - 1][child]));
1467 return this->restriction[refinement_case - 1][child];
1477 template <
class PolynomialType,
int dim,
int spacedim>
1480 const unsigned int shape_index,
1481 const unsigned int face_index)
const
1490 return (((shape_index == 0) && (face_index == 0)) ||
1491 ((shape_index == 1) && (face_index == 1)));
1495 if (((dim == 2) && (shape_index >= this->first_quad_index)) ||
1496 ((dim == 3) && (shape_index >= this->first_hex_index)))
1500 if (shape_index < this->first_line_index)
1505 const unsigned int vertex_no = shape_index;
1509 for (
unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face; ++v)
1516 else if (shape_index < this->first_quad_index)
1519 const unsigned int line_index =
1520 (shape_index - this->first_line_index) / this->dofs_per_line;
1526 return (line_index == face_index);
1530 const unsigned int lines_per_face =
1533 for (
unsigned int l = 0;
l < lines_per_face; ++
l)
1543 else if (shape_index < this->first_hex_index)
1546 const unsigned int quad_index =
1547 (shape_index - this->first_quad_index) / this->dofs_per_quad;
1548 Assert(
static_cast<signed int>(quad_index) <
1558 return (quad_index == face_index);
1577 template <
typename PolynomialType,
int dim,
int spacedim>
1578 std::pair<Table<2, bool>, std::vector<unsigned int>>
1586 for (
unsigned int i = 0; i < Utilities::fixed_power<dim>(q_degree + 1); ++i)
1587 constant_modes(0, i) =
true;
1588 return std::pair<Table<2, bool>, std::vector<unsigned int>>(
1589 constant_modes, std::vector<unsigned int>(1, 0));
1595 #include "fe_q_base.inst"