52 get_QGaussLobatto_points(
const unsigned int degree)
57 return std::vector<Point<1>>(1,
Point<1>(0.5));
65template <
int dim,
int spacedim>
99template <
int dim,
int spacedim>
106 polynomials.size() - 1,
109 polynomials.size() - 1),
111 polynomials.size() - 1)
117 polynomials.size() - 1)
131template <
int dim,
int spacedim>
139 std::ostringstream namebuf;
141 << this->degree <<
")";
142 return namebuf.str();
147template <
int dim,
int spacedim>
151 std::vector<double> & nodal_values)
const
154 this->get_unit_support_points().size());
158 for (
unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
162 nodal_values[i] = support_point_values[i](0);
168template <
int dim,
int spacedim>
169std::unique_ptr<FiniteElement<dim, spacedim>>
172 return std::make_unique<FE_DGQ<dim, spacedim>>(*this);
181template <
int dim,
int spacedim>
182std::vector<unsigned int>
185 std::vector<unsigned int> dpo(dim + 1, 0
U);
187 for (
unsigned int i = 1; i < dim; ++i)
194template <
int dim,
int spacedim>
197 const char direction)
const
199 const unsigned int n = this->degree + 1;
201 for (
unsigned int i = 1; i < dim; ++i)
210 for (
unsigned int i = n; i > 0;)
220 for (
unsigned int iz = 0; iz < ((dim > 2) ? n : 1); ++iz)
221 for (
unsigned int j = 0; j < n; ++j)
222 for (
unsigned int i = 0; i < n; ++i)
224 unsigned int k = n * i - j + n - 1 + n * n * iz;
231 for (
unsigned int iz = 0; iz < ((dim > 2) ? n : 1); ++iz)
232 for (
unsigned int iy = 0; iy < n; ++iy)
233 for (
unsigned int ix = 0; ix < n; ++ix)
235 unsigned int k = n * ix - iy + n - 1 + n * n * iz;
243 for (
unsigned int iz = 0; iz < n; ++iz)
244 for (
unsigned int iy = 0; iy < n; ++iy)
245 for (
unsigned int ix = 0; ix < n; ++ix)
247 unsigned int k = n * (n * iy - iz + n - 1) + ix;
255 for (
unsigned int iz = 0; iz < n; ++iz)
256 for (
unsigned int iy = 0; iy < n; ++iy)
257 for (
unsigned int ix = 0; ix < n; ++ix)
259 unsigned int k = n * (n * iy - iz + n - 1) + ix;
271template <
int dim,
int spacedim>
283 typename FE::ExcInterpolationNotImplemented());
290 Assert(interpolation_matrix.
m() == this->n_dofs_per_cell(),
292 this->n_dofs_per_cell()));
303 this->n_dofs_per_cell());
307 for (
unsigned int j = 0; j < this->n_dofs_per_cell(); ++j)
313 for (
unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
314 cell_interpolation(j, i) = this->poly_space->compute_value(i, p);
317 source_interpolation(j, i) = source_fe.
poly_space->compute_value(i, p);
324 cell_interpolation.gauss_jordan();
325 cell_interpolation.mmult(interpolation_matrix, source_interpolation);
328 for (
unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
330 if (
std::fabs(interpolation_matrix(i, j)) < 1
e-15)
331 interpolation_matrix(i, j) = 0.;
338 for (
unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
342 sum += interpolation_matrix(i, j);
351template <
int dim,
int spacedim>
356 const unsigned int)
const
364 (void)interpolation_matrix;
368 typename FE::ExcInterpolationNotImplemented());
370 Assert(interpolation_matrix.
m() == 0,
372 Assert(interpolation_matrix.
n() == 0,
378template <
int dim,
int spacedim>
384 const unsigned int)
const
392 (void)interpolation_matrix;
396 typename FE::ExcInterpolationNotImplemented());
398 Assert(interpolation_matrix.
m() == 0,
400 Assert(interpolation_matrix.
n() == 0,
406template <
int dim,
int spacedim>
409 const unsigned int child,
416 "Prolongation matrices are only available for refined cells!"));
420 if (this->prolongation[refinement_case - 1][child].n() == 0)
422 std::lock_guard<std::mutex> lock(this->mutex);
425 if (this->prolongation[refinement_case - 1][child].n() ==
426 this->n_dofs_per_cell())
427 return this->prolongation[refinement_case - 1][child];
435 std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
437 isotropic_matrices.back().resize(
440 this->n_dofs_per_cell()));
450 isotropic_matrices.back());
477 return this->prolongation[refinement_case - 1][child];
482template <
int dim,
int spacedim>
485 const unsigned int child,
492 "Restriction matrices are only available for refined cells!"));
496 if (this->restriction[refinement_case - 1][child].n() == 0)
498 std::lock_guard<std::mutex> lock(this->mutex);
501 if (this->restriction[refinement_case - 1][child].n() ==
502 this->n_dofs_per_cell())
503 return this->restriction[refinement_case - 1][child];
511 std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
513 isotropic_matrices.back().resize(
516 this->n_dofs_per_cell()));
525 this_nonconst.
restriction[refinement_case - 1].swap(
526 isotropic_matrices.back());
553 return this->restriction[refinement_case - 1][child];
558template <
int dim,
int spacedim>
567template <
int dim,
int spacedim>
568std::vector<std::pair<unsigned int, unsigned int>>
576 return std::vector<std::pair<unsigned int, unsigned int>>();
581template <
int dim,
int spacedim>
582std::vector<std::pair<unsigned int, unsigned int>>
590 return std::vector<std::pair<unsigned int, unsigned int>>();
595template <
int dim,
int spacedim>
596std::vector<std::pair<unsigned int, unsigned int>>
599 const unsigned int)
const
605 return std::vector<std::pair<unsigned int, unsigned int>>();
610template <
int dim,
int spacedim>
614 const unsigned int codim)
const
634 if (this->degree < fe_dgq_other->degree)
636 else if (this->degree == fe_dgq_other->degree)
644 if (this->degree < fe_q_other->degree)
646 else if (this->degree == fe_q_other->degree)
654 if (this->degree < fe_bernstein_other->degree)
656 else if (this->degree == fe_bernstein_other->degree)
664 if (this->degree < fe_bubbles_other->degree)
666 else if (this->degree == fe_bubbles_other->degree)
674 if (this->degree < fe_dg0_other->degree)
676 else if (this->degree == fe_dg0_other->degree)
684 if (this->degree < fe_q_iso_q1_other->degree)
686 else if (this->degree == fe_q_iso_q1_other->degree)
694 if (this->degree < fe_hierarchical_other->degree)
696 else if (this->degree == fe_hierarchical_other->degree)
704 if (this->degree < fe_dgp_other->degree)
706 else if (this->degree == fe_dgp_other->degree)
714 if (fe_nothing->is_dominating())
729template <
int dim,
int spacedim>
732 const unsigned int face_index)
const
737 unsigned int n = this->degree + 1;
748 bool support_points_on_boundary =
true;
749 for (
unsigned int d = 0;
d < dim; ++
d)
751 support_points_on_boundary =
false;
752 for (
unsigned int d = 0;
d < dim; ++
d)
754 support_points_on_boundary =
false;
755 if (support_points_on_boundary ==
false)
758 unsigned int n2 = n * n;
770 return (((shape_index == 0) && (face_index == 0)) ||
771 ((shape_index == this->degree) && (face_index == 1)));
776 if (face_index == 0 && (shape_index % n) == 0)
778 if (face_index == 1 && (shape_index % n) == this->degree)
780 if (face_index == 2 && shape_index < n)
782 if (face_index == 3 && shape_index >= this->n_dofs_per_cell() - n)
789 const unsigned int in2 = shape_index % n2;
792 if (face_index == 0 && (shape_index % n) == 0)
795 if (face_index == 1 && (shape_index % n) == n - 1)
798 if (face_index == 2 && in2 < n)
801 if (face_index == 3 && in2 >= n2 - n)
804 if (face_index == 4 && shape_index < n2)
807 if (face_index == 5 && shape_index >= this->n_dofs_per_cell() - n2)
820template <
int dim,
int spacedim>
821std::pair<Table<2, bool>, std::vector<unsigned int>>
825 constant_modes.fill(
true);
826 return std::pair<Table<2, bool>, std::vector<unsigned int>>(
827 constant_modes, std::vector<unsigned int>(1, 0));
834template <
int dim,
int spacedim>
847template <
int dim,
int spacedim>
853 std::ostringstream namebuf;
854 bool equidistant =
true;
855 std::vector<double> points(this->degree + 1);
857 auto *
const polynomial_space =
860 std::vector<unsigned int> lexicographic =
861 polynomial_space->get_numbering_inverse();
862 for (
unsigned int j = 0; j <= this->degree; j++)
866 for (
unsigned int j = 0; j <= this->degree; j++)
867 if (
std::abs(points[j] -
static_cast<double>(j) / this->degree) > 1
e-15)
872 if (this->degree == 0 &&
std::abs(points[0] - 0.5) < 1
e-15)
875 if (equidistant ==
true)
877 if (this->degree > 2)
878 namebuf <<
"FE_DGQArbitraryNodes<"
880 <<
">(QIterated(QTrapezoid()," << this->degree <<
"))";
883 << this->degree <<
")";
884 return namebuf.str();
889 bool gauss_lobatto =
true;
890 for (
unsigned int j = 0; j <= this->degree; j++)
891 if (points[j] != points_gl.
point(j)(0))
893 gauss_lobatto =
false;
897 if (gauss_lobatto ==
true)
900 << this->degree <<
")";
901 return namebuf.str();
905 const QGauss<1> points_g(this->degree + 1);
907 for (
unsigned int j = 0; j <= this->degree; j++)
908 if (points[j] != points_g.
point(j)(0))
917 <<
">(QGauss(" << this->degree + 1 <<
"))";
918 return namebuf.str();
923 bool gauss_log =
true;
924 for (
unsigned int j = 0; j <= this->degree; j++)
925 if (points[j] != points_glog.
point(j)(0))
931 if (gauss_log ==
true)
934 <<
">(QGaussLog(" << this->degree + 1 <<
"))";
935 return namebuf.str();
940 <<
">(QUnknownNodes(" << this->degree + 1 <<
"))";
941 return namebuf.str();
946template <
int dim,
int spacedim>
951 std::vector<double> & nodal_values)
const
954 this->get_unit_support_points().size());
958 for (
unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
962 nodal_values[i] = support_point_values[i](0);
968template <
int dim,
int spacedim>
969std::unique_ptr<FiniteElement<dim, spacedim>>
973 std::vector<Point<1>> qpoints(this->degree + 1);
974 auto *
const polynomial_space =
977 std::vector<unsigned int> lexicographic =
978 polynomial_space->get_numbering_inverse();
979 for (
unsigned int i = 0; i <= this->degree; ++i)
983 return std::make_unique<FE_DGQArbitraryNodes<dim, spacedim>>(pquadrature);
990template <
int dim,
int spacedim>
993 Polynomials::Legendre::generate_complete_basis(degree))
998template <
int dim,
int spacedim>
999std::pair<Table<2, bool>, std::vector<unsigned int>>
1005 constant_modes(0, 0) =
true;
1006 return std::pair<Table<2, bool>, std::vector<unsigned int>>(
1007 constant_modes, std::vector<unsigned int>(1, 0));
1012template <
int dim,
int spacedim>
1022template <
int dim,
int spacedim>
1023std::unique_ptr<FiniteElement<dim, spacedim>>
1026 return std::make_unique<FE_DGQLegendre<dim, spacedim>>(this->degree);
1033template <
int dim,
int spacedim>
1036 Polynomials::HermiteLikeInterpolation::generate_complete_basis(degree))
1041template <
int dim,
int spacedim>
1051template <
int dim,
int spacedim>
1052std::unique_ptr<FiniteElement<dim, spacedim>>
1055 return std::make_unique<FE_DGQHermite<dim, spacedim>>(this->degree);
1061#include "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 std::unique_ptr< ScalarPolynomialsBase< dim > > poly_space
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 Point< dim > & point(const unsigned int i) const
const std::vector< Point< dim > > & get_points() const
unsigned int size() const
SymmetricTensor< rank, dim, Number > sum(const SymmetricTensor< rank, dim, Number > &local, const MPI_Comm &mpi_communicator)
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcNotImplemented()
#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)
Expression fabs(const Expression &x)
std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
@ either_element_can_dominate
@ other_element_dominates
@ neither_element_dominates
void L2(Vector< number > &result, const FEValuesBase< dim > &fe, const std::vector< double > &input, const double factor=1.)
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)
std::vector< Polynomial< double > > generate_complete_Lagrange_basis(const std::vector< Point< 1 > > &points)
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)
std::vector< Point< dim > > unit_support_points(const std::vector< Point< 1 > > &line_support_points, const std::vector< unsigned int > &renumbering)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)