57 for (
unsigned int d = 0; d < dim; ++d)
71 for (
unsigned int d = 0; d < dim - 1; ++d)
88template <
int xdim,
int xspacedim>
95 template <
int spacedim>
104 template <
int spacedim>
109 const unsigned int dim = 2;
111 unsigned int q_deg = fe.degree;
113 &fe.get_poly_space()) !=
nullptr)
114 q_deg = fe.degree - 1;
165 for (
unsigned int i = 1; i <
q_deg; ++i)
170 for (
unsigned int i = 1; i <
q_deg; ++i)
178 fe.interface_constraints.TableBase<2, double>::reinit(
179 fe.interface_constraints_size());
183 const std::vector<unsigned int> &index_map_inverse =
184 fe.get_poly_space_numbering_inverse();
188 fe.poly_space->compute_value(index_map_inverse[0],
Point<dim>()) -
193 for (
unsigned int j = 0;
j <
q_deg + 1; ++
j)
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;
219 &fe.get_poly_space()) !=
nullptr)
220 q_deg = fe.degree - 1;
250 const unsigned int n =
q_deg - 1;
251 const double step = 1. /
q_deg;
252 std::vector<Point<1>> line_support_points(n);
253 for (
unsigned int i = 0; i < n; ++i)
254 line_support_points[i][0] = (i + 1) * step;
258 auto get_points = [&](
const unsigned int face_no,
261 ReferenceCells::get_hypercube<2>(),
284 for (
unsigned int face = 0;
287 for (
unsigned int subface = 0;
291 const auto p_line = get_points(face, subface);
299 for (
unsigned int i = 0,
iy = 1;
iy <= n; ++
iy)
300 for (
unsigned int ix = 1;
ix <= n; ++
ix)
304 for (
unsigned int child = 0;
319 const std::vector<unsigned int> &index_map_inverse =
320 fe.get_poly_space_numbering_inverse();
324 fe.poly_space->compute_value(index_map_inverse[0],
Point<dim>()) -
328 fe.interface_constraints.TableBase<2, double>::reinit(
329 fe.interface_constraints_size());
333 const double interval =
static_cast<double>(
q_deg * 2);
346 for (
unsigned int k = 0;
k < dim - 1; ++
k)
377 for (
unsigned int j = 0;
j <
pnts; ++
j)
379 unsigned int indices[2] = {
j % (
q_deg + 1),
j / (
q_deg + 1)};
381 for (
unsigned int k = 0;
k < 2; ++
k)
383 indices[
k] =
q_deg - indices[
k];
386 indices[1] * (
q_deg + 1) + indices[0];
389 fe.poly_space->compute_value(index_map_inverse[
new_index],
406template <
int dim,
int spacedim>
410 const std::vector<bool> &restriction_is_additive_flags)
414 restriction_is_additive_flags,
424template <
int dim,
int spacedim>
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."));
437 Utilities::fixed_power<dim>(
q_degree + 1);
443 [
this, q_dofs_per_cell]() {
444 std::vector<unsigned int> renumber =
445 FETools::hierarchic_to_lexicographic_numbering<dim>(
q_degree);
447 renumber.push_back(i);
500template <
int dim,
int spacedim>
513 this->n_dofs_per_cell()));
520 Utilities::fixed_power<dim>(
q_degree + 1);
522 Utilities::fixed_power<dim>(
source_fe->degree + 1);
540 source_fe->poly_space->compute_value(i, p);
558 for (
unsigned int j = 0;
j <
source_fe->n_dofs_per_cell(); ++
j)
568 for (
unsigned int j = 0;
j <
source_fe->n_dofs_per_cell(); ++
j)
600template <
int dim,
int spacedim>
605 const unsigned int face_no)
const
615template <
int dim,
int spacedim>
619 const unsigned int subface,
621 const unsigned int face_no)
const
725template <
int dim,
int spacedim>
734template <
int dim,
int spacedim>
735std::vector<std::pair<unsigned int, unsigned int>>
767 else if (
fe_other.n_unique_faces() == 1 &&
fe_other.n_dofs_per_face(0) == 0)
787template <
int dim,
int spacedim>
788std::vector<std::pair<unsigned int, unsigned int>>
805 const unsigned int p = this->
degree;
808 std::vector<std::pair<unsigned int, unsigned int>>
identities;
810 const std::vector<unsigned int> &index_map_inverse =
813 fe_q_other->get_poly_space_numbering_inverse();
815 for (
unsigned int i = 0; i < p - 1; ++i)
816 for (
unsigned int j = 0;
j <
q - 1; ++
j)
845 std::vector<std::pair<unsigned int, unsigned int>>
identities;
847 for (
unsigned int i = 0; i < this->
degree - 1; ++i)
851 fe_p_other->get_unit_support_points()[
j + 3][0]) < 1e-14)
862 else if (
fe_other.n_unique_faces() == 1 &&
fe_other.n_dofs_per_face(0) == 0)
882template <
int dim,
int spacedim>
883std::vector<std::pair<unsigned int, unsigned int>>
886 const unsigned int)
const
898 const unsigned int p = this->
degree;
901 std::vector<std::pair<unsigned int, unsigned int>>
identities;
903 const std::vector<unsigned int> &index_map_inverse =
906 fe_q_other->get_poly_space_numbering_inverse();
908 for (
unsigned int i1 = 0;
i1 < p - 1; ++
i1)
909 for (
unsigned int i2 = 0;
i2 < p - 1; ++
i2)
910 for (
unsigned int j1 = 0;
j1 <
q - 1; ++
j1)
911 for (
unsigned int j2 = 0;
j2 <
q - 1; ++
j2)
930 return std::vector<std::pair<unsigned int, unsigned int>>();
932 else if (
fe_other.n_unique_faces() == 1 &&
fe_other.n_dofs_per_face(0) == 0)
941 return std::vector<std::pair<unsigned int, unsigned int>>();
946 return std::vector<std::pair<unsigned int, unsigned int>>();
958template <
int dim,
int spacedim>
961 const std::vector<
Point<1>> &points)
963 const std::vector<unsigned int> &index_map_inverse =
976 for (
unsigned int k = 0;
k < support_quadrature.size(); ++
k)
978 support_quadrature.point(
k);
983template <
int dim,
int spacedim>
986 const std::vector<
Point<1>> &points)
991 const unsigned int face_no = 0;
994 Utilities::fixed_power<dim - 1>(
q_degree + 1));
1015 for (
unsigned int k = 0;
k < support_quadrature.size(); ++
k)
1017 support_quadrature.point(
k);
1022template <
int dim,
int spacedim>
1038 const unsigned int face_no = 0;
1046 const unsigned int n =
q_degree - 1;
1067 for (
unsigned int local = 0; local < this->
n_dofs_per_quad(face_no); ++local)
1071 unsigned int i = local % n,
j = local / n;
1080 i + (n - 1 -
j) * n - local;
1084 (n - 1 -
j) + (n - 1 - i) * n - local;
1088 (n - 1 - i) +
j * n - local;
1095 j + (n - 1 - i) * n - local;
1099 (n - 1 - i) + (n - 1 -
j) * n - local;
1103 (n - 1 -
j) + i * n - local;
1109template <
int dim,
int spacedim>
1112 const unsigned int face_index,
1113 const unsigned int face,
1135 dof_index_on_vertex;
1142 const unsigned int index =
1174 combined_orientation ==
1189 adjusted_dof_index_on_line);
1197 const unsigned int index =
1212template <
int dim,
int spacedim>
1213std::vector<unsigned int>
1218 std::vector<unsigned int>
dpo(dim + 1, 1U);
1219 for (
unsigned int i = 1; i <
dpo.size(); ++i)
1226template <
int dim,
int spacedim>
1229 const std::vector<
Point<1>> &points)
1236template <
int dim,
int spacedim>
1239 const unsigned int child,
1246 "Prolongation matrices are only available for refined cells!"));
1251 if (this->
prolongation[refinement_case - 1][child].n() == 0)
1256 if (this->
prolongation[refinement_case - 1][child].n() ==
1262 Utilities::fixed_power<dim>(
q_degree + 1);
1278 i,
this->unit_support_points[i])) < eps,
1280 "to one or zero in a nodal point. "
1281 "This typically indicates that the "
1282 "polynomial interpolation is "
1283 "ill-conditioned such that round-off "
1284 "prevents the sum to be one."));
1288 i,
this->unit_support_points[
j])) < eps,
1290 "The Lagrange polynomial does not evaluate "
1291 "to one or zero in a nodal point. "
1292 "This typically indicates that the "
1293 "polynomial interpolation is "
1294 "ill-conditioned such that round-off "
1295 "prevents the sum to be one."));
1307 const std::vector<unsigned int> &index_map_inverse =
1314 unsigned int factor = 1;
1315 for (
unsigned int d = 0;
d < dim; ++
d)
1327 for (
unsigned int j = 0;
j <
dofs1d; ++
j)
1335 for (
unsigned int i = 0; i <
dofs1d; ++i)
1336 for (
unsigned int d = 0;
d < dim; ++
d)
1342 this->
poly_space->compute_value(index_map_inverse[i], point);
1371 internal::FE_Q_Base::zero_indices<dim>(
j_indices);
1375 internal::FE_Q_Base::zero_indices<dim>(
i_indices);
1379 for (
unsigned int d = 1;
d < dim; ++
d)
1387 const unsigned int j_ind = index_map_inverse[
j +
jj];
1389 prolongate(
j_ind, index_map_inverse[i +
ii]) =
1413 sum += prolongate(row, col);
1414 Assert(std::fabs(sum - 1.) <
1417 "prolongation matrix do not add to one. "
1418 "This typically indicates that the "
1419 "polynomial interpolation is "
1420 "ill-conditioned such that round-off "
1421 "prevents the sum to be one."));
1427 this->
prolongation[refinement_case - 1][child]) = std::move(prolongate);
1436template <
int dim,
int spacedim>
1439 const unsigned int child,
1446 "Restriction matrices are only available for refined cells!"));
1451 if (this->
restriction[refinement_case - 1][child].n() == 0)
1456 if (this->
restriction[refinement_case - 1][child].n() ==
1458 return this->
restriction[refinement_case - 1][child];
1463 const unsigned int q_dofs_per_cell =
1464 Utilities::fixed_power<dim>(
q_degree + 1);
1485 const std::vector<unsigned int> &index_map_inverse =
1488 const unsigned int dofs1d =
q_degree + 1;
1495 unsigned int mother_dof = index_map_inverse[i];
1510 for (
unsigned int j = 0;
j <
dofs1d; ++
j)
1511 for (
unsigned int d = 0;
d < dim; ++
d)
1516 this->
poly_space->compute_value(index_map_inverse[
j],
1520 internal::FE_Q_Base::zero_indices<dim>(
j_indices);
1525 for (
unsigned int d = 1;
d < dim; ++
d)
1536 const unsigned int child_dof = index_map_inverse[
j +
jj];
1537 if (std::fabs(val - 1.) < eps)
1539 else if (std::fabs(val) > eps)
1543 internal::FE_Q_Base::increment_indices<dim>(
j_indices,
1551 "restriction matrix do not add to one. "
1552 "This typically indicates that the "
1553 "polynomial interpolation is "
1554 "ill-conditioned such that round-off "
1555 "prevents the sum to be one."));
1572 return this->
restriction[refinement_case - 1][child];
1582template <
int dim,
int spacedim>
1586 const unsigned int face_index)
const
1595 return (((
shape_index == 0) && (face_index == 0)) ||
1625 const unsigned int line_index =
1631 if constexpr (dim == 2)
1632 return (line_index == face_index);
1633 else if constexpr (dim == 3)
1649 const unsigned int quad_index =
1660 return (quad_index == face_index);
1679template <
int dim,
int spacedim>
1680std::pair<Table<2, bool>, std::vector<unsigned int>>
1689 constant_modes(0, i) =
true;
1690 return std::pair<Table<2, bool>, std::vector<unsigned int>>(
1691 constant_modes, std::vector<unsigned int>(1, 0));
1697#include "fe/fe_q_base.inst"
std::vector< unsigned int > get_poly_space_numbering_inverse() const
virtual double shape_value(const unsigned int i, const Point< dim > &p) const override
const std::unique_ptr< ScalarPolynomialsBase< dim > > poly_space
Threads::Mutex prolongation_matrix_mutex
void initialize_unit_face_support_points(const std::vector< Point< 1 > > &points)
virtual bool hp_constraints_are_implemented() const override
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
virtual void get_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix) const override
Threads::Mutex restriction_matrix_mutex
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_quad_dof_identities(const FiniteElement< dim, spacedim > &fe_other, const unsigned int face_no=0) const override
const unsigned int q_degree
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const override
void initialize(const std::vector< Point< 1 > > &support_points_1d)
void initialize_unit_support_points(const std::vector< Point< 1 > > &points)
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 override
virtual void get_face_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix, const unsigned int face_no=0) const override
static std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
void initialize_constraints(const std::vector< Point< 1 > > &points)
void initialize_dof_index_permutations()
virtual const FullMatrix< double > & get_restriction_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const override
virtual const FullMatrix< double > & get_prolongation_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const override
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const override
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 override
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_vertex_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
FE_Q_Base(const ScalarPolynomialsBase< dim > &poly_space, const FiniteElementData< dim > &fe_data, const std::vector< bool > &restriction_is_additive_flags)
unsigned int get_first_line_index() 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 get_first_quad_index(const unsigned int quad_no=0) const
unsigned int n_dofs_per_face(unsigned int face_no=0, unsigned int child=0) const
ReferenceCell reference_cell() const
unsigned int get_first_face_line_index(const unsigned int face_no=0) const
unsigned int get_first_face_quad_index(const unsigned int face_no=0) const
unsigned int n_unique_faces() const
unsigned int n_dofs_per_quad(unsigned int face_no=0) const
unsigned int get_first_hex_index() const
std::vector< std::vector< Point< dim - 1 > > > unit_face_support_points
std::vector< std::vector< FullMatrix< double > > > restriction
std::vector< Table< 2, int > > adjust_quad_dof_index_for_face_orientation_table
std::vector< int > adjust_line_dof_index_for_line_orientation_table
std::vector< Point< dim > > unit_support_points
std::vector< std::vector< FullMatrix< double > > > prolongation
static void project_to_subface(const ReferenceCell &reference_cell, const SubQuadrature &quadrature, const unsigned int face_no, const unsigned int subface_no, std::vector< Point< dim > > &q_points, const RefinementCase< dim - 1 > &ref_case=RefinementCase< dim - 1 >::isotropic_refinement)
static void project_to_face(const ReferenceCell &reference_cell, const SubQuadrature &quadrature, const unsigned int face_no, std::vector< Point< dim > > &q_points)
unsigned int face_to_cell_vertices(const unsigned int face, const unsigned int vertex, const types::geometric_orientation face_orientation) const
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcInterpolationNotImplemented()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Task< RT > new_task(const std::function< RT()> &function)
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
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)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
T sum(const T &t, const MPI_Comm mpi_communicator)
types::geometric_orientation combined_face_orientation(const bool face_orientation, const bool face_rotation, const bool face_flip)
constexpr unsigned int invalid_unsigned_int
constexpr types::geometric_orientation default_geometric_orientation
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
static void initialize_constraints(const std::vector< Point< 1 > > &, FE_Q_Base< 3, spacedim > &fe)
static void initialize_constraints(const std::vector< Point< 1 > > &, FE_Q_Base< 1, spacedim > &)
static void initialize_constraints(const std::vector< Point< 1 > > &, FE_Q_Base< 2, spacedim > &fe)
static Point< dim > cell_to_child_coordinates(const Point< dim > &p, const unsigned int child_index, const RefinementCase< dim > refine_case=RefinementCase< dim >::isotropic_refinement)
static Point< dim > child_to_cell_coordinates(const Point< dim > &p, const unsigned int child_index, const RefinementCase< dim > refine_case=RefinementCase< dim >::isotropic_refinement)