51 get_QGaussLobatto_points(
const unsigned int degree)
56 return std::vector<Point<1>>(1,
Point<1>(0.5));
64template <
int dim,
int spacedim>
98template <
int dim,
int spacedim>
105 polynomials.
size() - 1,
108 polynomials.
size() - 1),
110 polynomials.
size() - 1)
116 polynomials.
size() - 1)
130template <
int dim,
int spacedim>
138 std::ostringstream namebuf;
140 << this->degree <<
")";
141 return namebuf.str();
146template <
int dim,
int spacedim>
150 std::vector<double> &nodal_values)
const
153 this->get_unit_support_points().size());
157 for (
unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
161 nodal_values[i] = support_point_values[i](0);
167template <
int dim,
int spacedim>
168std::unique_ptr<FiniteElement<dim, spacedim>>
171 return std::make_unique<FE_DGQ<dim, spacedim>>(*this);
180template <
int dim,
int spacedim>
181std::vector<unsigned int>
184 std::vector<unsigned int> dpo(dim + 1, 0U);
186 for (
unsigned int i = 1; i < dim; ++i)
193template <
int dim,
int spacedim>
196 const char direction)
const
198 const unsigned int n = this->degree + 1;
200 for (
unsigned int i = 1; i < dim; ++i)
209 for (
unsigned int i = n; i > 0;)
219 for (
unsigned int iz = 0; iz < ((dim > 2) ? n : 1); ++iz)
220 for (
unsigned int j = 0; j < n; ++j)
221 for (
unsigned int i = 0; i < n; ++i)
223 unsigned int k = n * i - j + n - 1 + n * n * iz;
230 for (
unsigned int iz = 0; iz < ((dim > 2) ? n : 1); ++iz)
231 for (
unsigned int iy = 0; iy < n; ++iy)
232 for (
unsigned int ix = 0; ix < n; ++ix)
234 unsigned int k = n * ix - iy + n - 1 + n * n * iz;
242 for (
unsigned int iz = 0; iz < n; ++iz)
243 for (
unsigned int iy = 0; iy < n; ++iy)
244 for (
unsigned int ix = 0; ix < n; ++ix)
246 unsigned int k = n * (n * iy - iz + n - 1) + ix;
254 for (
unsigned int iz = 0; iz < n; ++iz)
255 for (
unsigned int iy = 0; iy < n; ++iy)
256 for (
unsigned int ix = 0; ix < n; ++ix)
258 unsigned int k = n * (n * iy - iz + n - 1) + ix;
270template <
int dim,
int spacedim>
280 Assert(interpolation_matrix.
m() == this->n_dofs_per_cell(),
282 this->n_dofs_per_cell()));
283 Assert(interpolation_matrix.
n() == source_fe->n_dofs_per_cell(),
285 source_fe->n_dofs_per_cell()));
293 this->n_dofs_per_cell());
295 source_fe->n_dofs_per_cell());
297 source_fe->n_dofs_per_cell());
298 for (
unsigned int j = 0; j < this->n_dofs_per_cell(); ++j)
304 for (
unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
305 cell_interpolation(j, i) = this->poly_space->compute_value(i, p);
307 for (
unsigned int i = 0; i < source_fe->n_dofs_per_cell(); ++i)
308 source_interpolation(j, i) =
309 source_fe->poly_space->compute_value(i, p);
316 cell_interpolation.gauss_jordan();
317 cell_interpolation.mmult(interpolation_matrix, source_interpolation);
320 for (
unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
321 for (
unsigned int j = 0; j < source_fe->n_dofs_per_cell(); ++j)
322 if (std::fabs(interpolation_matrix(i, j)) < 1e-15)
323 interpolation_matrix(i, j) = 0.;
332 for (
unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
335 for (
unsigned int j = 0; j < source_fe->n_dofs_per_cell(); ++j)
336 sum += interpolation_matrix(i, j);
338 Assert(std::fabs(sum - 1) <
339 5e-14 *
std::max(this->degree, 1U) * dim,
364 spacedim>::ExcInterpolationNotImplemented()));
369template <
int dim,
int spacedim>
374 const unsigned int)
const
382 (void)interpolation_matrix;
386 typename FE::ExcInterpolationNotImplemented());
388 Assert(interpolation_matrix.
m() == 0,
390 Assert(interpolation_matrix.
n() == 0,
396template <
int dim,
int spacedim>
402 const unsigned int)
const
410 (void)interpolation_matrix;
414 typename FE::ExcInterpolationNotImplemented());
416 Assert(interpolation_matrix.
m() == 0,
418 Assert(interpolation_matrix.
n() == 0,
424template <
int dim,
int spacedim>
427 const unsigned int child,
434 "Prolongation matrices are only available for refined cells!"));
436 child, this->reference_cell().
template n_children<dim>(refinement_case));
439 if (this->prolongation[refinement_case - 1][child].n() == 0)
441 std::lock_guard<std::mutex> lock(prolongation_matrix_mutex);
444 if (this->prolongation[refinement_case - 1][child].n() ==
445 this->n_dofs_per_cell())
446 return this->prolongation[refinement_case - 1][child];
454 std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
456 isotropic_matrices.back().resize(
457 this->reference_cell().
template n_children<dim>(
460 this->n_dofs_per_cell()));
470 std::move(isotropic_matrices.back());
497 return this->prolongation[refinement_case - 1][child];
502template <
int dim,
int spacedim>
505 const unsigned int child,
512 "Restriction matrices are only available for refined cells!"));
514 child, this->reference_cell().
template n_children<dim>(refinement_case));
517 if (this->restriction[refinement_case - 1][child].n() == 0)
519 std::lock_guard<std::mutex> lock(restriction_matrix_mutex);
522 if (this->restriction[refinement_case - 1][child].n() ==
523 this->n_dofs_per_cell())
524 return this->restriction[refinement_case - 1][child];
532 std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
534 isotropic_matrices.back().resize(
535 this->reference_cell().
template n_children<dim>(
538 this->n_dofs_per_cell()));
548 std::move(isotropic_matrices.back());
575 return this->restriction[refinement_case - 1][child];
580template <
int dim,
int spacedim>
589template <
int dim,
int spacedim>
590std::vector<std::pair<unsigned int, unsigned int>>
598 return std::vector<std::pair<unsigned int, unsigned int>>();
603template <
int dim,
int spacedim>
604std::vector<std::pair<unsigned int, unsigned int>>
612 return std::vector<std::pair<unsigned int, unsigned int>>();
617template <
int dim,
int spacedim>
618std::vector<std::pair<unsigned int, unsigned int>>
621 const unsigned int)
const
627 return std::vector<std::pair<unsigned int, unsigned int>>();
632template <
int dim,
int spacedim>
636 const unsigned int codim)
const
656 if (this->degree < fe_dgq_other->degree)
658 else if (this->degree == fe_dgq_other->degree)
666 if (this->degree < fe_q_other->degree)
668 else if (this->degree == fe_q_other->degree)
676 if (this->degree < fe_bernstein_other->degree)
678 else if (this->degree == fe_bernstein_other->degree)
686 if (this->degree < fe_bubbles_other->degree)
688 else if (this->degree == fe_bubbles_other->degree)
696 if (this->degree < fe_dg0_other->degree)
698 else if (this->degree == fe_dg0_other->degree)
706 if (this->degree < fe_q_iso_q1_other->degree)
708 else if (this->degree == fe_q_iso_q1_other->degree)
716 if (this->degree < fe_hierarchical_other->degree)
718 else if (this->degree == fe_hierarchical_other->degree)
726 if (this->degree < fe_dgp_other->degree)
728 else if (this->degree == fe_dgp_other->degree)
736 if (fe_nothing->is_dominating())
751template <
int dim,
int spacedim>
754 const unsigned int face_index)
const
759 unsigned int n = this->degree + 1;
764 if (this->unit_support_points.empty())
770 bool support_points_on_boundary =
true;
771 for (
unsigned int d = 0; d < dim; ++d)
772 if (
std::abs(this->unit_support_points[0][d]) > 1e-13)
773 support_points_on_boundary =
false;
774 for (
unsigned int d = 0; d < dim; ++d)
775 if (
std::abs(this->unit_support_points.back()[d] - 1.) > 1e-13)
776 support_points_on_boundary =
false;
777 if (support_points_on_boundary ==
false)
780 unsigned int n2 = n * n;
792 return (((shape_index == 0) && (face_index == 0)) ||
793 ((shape_index == this->degree) && (face_index == 1)));
798 if (face_index == 0 && (shape_index % n) == 0)
800 if (face_index == 1 && (shape_index % n) == this->degree)
802 if (face_index == 2 && shape_index < n)
804 if (face_index == 3 && shape_index >= this->n_dofs_per_cell() - n)
811 const unsigned int in2 = shape_index % n2;
814 if (face_index == 0 && (shape_index % n) == 0)
817 if (face_index == 1 && (shape_index % n) == n - 1)
820 if (face_index == 2 && in2 < n)
823 if (face_index == 3 && in2 >= n2 - n)
826 if (face_index == 4 && shape_index < n2)
829 if (face_index == 5 && shape_index >= this->n_dofs_per_cell() - n2)
842template <
int dim,
int spacedim>
843std::pair<Table<2, bool>, std::vector<unsigned int>>
847 constant_modes.fill(
true);
848 return std::pair<Table<2, bool>, std::vector<unsigned int>>(
849 constant_modes, std::vector<unsigned int>(1, 0));
856template <
int dim,
int spacedim>
860 Polynomials::generate_complete_Lagrange_basis(points.get_points()))
869template <
int dim,
int spacedim>
875 std::ostringstream namebuf;
876 bool equidistant =
true;
877 std::vector<double> points(this->degree + 1);
879 const auto *
const polynomial_space =
882 std::vector<unsigned int> lexicographic =
883 polynomial_space->get_numbering_inverse();
884 for (
unsigned int j = 0; j <= this->degree; ++j)
885 points[j] = this->unit_support_points[lexicographic[j]][0];
888 for (
unsigned int j = 0; j <= this->degree; ++j)
889 if (
std::abs(points[j] -
static_cast<double>(j) / this->degree) > 1e-15)
894 if (this->degree == 0 &&
std::abs(points[0] - 0.5) < 1e-15)
897 if (equidistant ==
true)
899 if (this->degree > 2)
900 namebuf <<
"FE_DGQArbitraryNodes<"
902 <<
">(QIterated(QTrapezoid()," << this->degree <<
"))";
905 << this->degree <<
")";
906 return namebuf.str();
911 bool gauss_lobatto =
true;
912 for (
unsigned int j = 0; j <= this->degree; ++j)
913 if (points[j] != points_gl.point(j)[0])
915 gauss_lobatto =
false;
919 if (gauss_lobatto ==
true)
922 << this->degree <<
")";
923 return namebuf.str();
927 const QGauss<1> points_g(this->degree + 1);
929 for (
unsigned int j = 0; j <= this->degree; ++j)
930 if (points[j] != points_g.point(j)[0])
939 <<
">(QGauss(" << this->degree + 1 <<
"))";
940 return namebuf.str();
945 bool gauss_log =
true;
946 for (
unsigned int j = 0; j <= this->degree; ++j)
947 if (points[j] != points_glog.point(j)[0])
953 if (gauss_log ==
true)
956 <<
">(QGaussLog(" << this->degree + 1 <<
"))";
957 return namebuf.str();
962 <<
">(QUnknownNodes(" << this->degree + 1 <<
"))";
963 return namebuf.str();
968template <
int dim,
int spacedim>
973 std::vector<double> &nodal_values)
const
976 this->get_unit_support_points().size());
980 for (
unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
984 nodal_values[i] = support_point_values[i](0);
990template <
int dim,
int spacedim>
991std::unique_ptr<FiniteElement<dim, spacedim>>
995 std::vector<Point<1>> qpoints(this->degree + 1);
996 const auto *
const polynomial_space =
999 std::vector<unsigned int> lexicographic =
1000 polynomial_space->get_numbering_inverse();
1001 for (
unsigned int i = 0; i <= this->degree; ++i)
1002 qpoints[i] =
Point<1>(this->unit_support_points[lexicographic[i]][0]);
1005 return std::make_unique<FE_DGQArbitraryNodes<dim, spacedim>>(pquadrature);
1012template <
int dim,
int spacedim>
1015 Polynomials::Legendre::generate_complete_basis(degree))
1020template <
int dim,
int spacedim>
1021std::pair<Table<2, bool>, std::vector<unsigned int>>
1027 constant_modes(0, 0) =
true;
1028 return std::pair<Table<2, bool>, std::vector<unsigned int>>(
1029 constant_modes, std::vector<unsigned int>(1, 0));
1034template <
int dim,
int spacedim>
1044template <
int dim,
int spacedim>
1045std::unique_ptr<FiniteElement<dim, spacedim>>
1048 return std::make_unique<FE_DGQLegendre<dim, spacedim>>(this->degree);
1055template <
int dim,
int spacedim>
1058 Polynomials::HermiteLikeInterpolation::generate_complete_basis(degree))
1063template <
int dim,
int spacedim>
1073template <
int dim,
int spacedim>
1074std::unique_ptr<FiniteElement<dim, spacedim>>
1077 return std::make_unique<FE_DGQHermite<dim, spacedim>>(this->degree);
1083#include "fe/fe_dgq.inst"
virtual std::string get_name() const override
FE_DGQArbitraryNodes(const Quadrature< 1 > &points)
virtual void convert_generalized_support_point_values_to_dof_values(const std::vector< Vector< double > > &support_point_values, std::vector< double > &nodal_values) const override
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
virtual std::string get_name() const override
FE_DGQHermite(const unsigned int degree)
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const override
FE_DGQLegendre(const unsigned int degree)
virtual std::string get_name() const override
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const override
virtual bool hp_constraints_are_implemented() const override
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
static std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const override
virtual const FullMatrix< double > & get_restriction_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const override
virtual void get_face_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix, const unsigned int face_no=0) const override
virtual FiniteElementDomination::Domination compare_for_domination(const FiniteElement< dim, spacedim > &fe_other, const unsigned int codim=0) const override final
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_vertex_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
virtual std::string get_name() const override
void rotate_indices(std::vector< unsigned int > &indices, const char direction) const
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
virtual void convert_generalized_support_point_values_to_dof_values(const std::vector< Vector< double > > &support_point_values, std::vector< double > &nodal_values) const override
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() 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 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_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix) const override
const unsigned int degree
unsigned int n_dofs_per_cell() const
std::vector< std::vector< FullMatrix< double > > > restriction
void reinit_restriction_and_prolongation_matrices(const bool isotropic_restriction_only=false, const bool isotropic_prolongation_only=false)
std::vector< Point< dim > > unit_support_points
std::vector< std::vector< FullMatrix< double > > > prolongation
const std::vector< Point< dim > > & get_points() const
unsigned int size() const
#define DEAL_II_NAMESPACE_OPEN
constexpr bool running_in_debug_mode()
#define DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NOT_IMPLEMENTED()
#define Assert(cond, exc)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
@ either_element_can_dominate
@ other_element_dominates
@ neither_element_dominates
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
std::string dim_string(const int dim, const int spacedim)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)