17 #include <deal.II/base/quadrature.h> 18 #include <deal.II/base/quadrature_lib.h> 19 #include <deal.II/base/std_cxx14/memory.h> 21 #include <deal.II/fe/fe.h> 22 #include <deal.II/fe/fe_dgq.h> 23 #include <deal.II/fe/fe_tools.h> 25 #include <deal.II/lac/vector.h> 31 DEAL_II_NAMESPACE_OPEN
41 get_QGaussLobatto_points(
const unsigned int degree)
46 return std::vector<Point<1>>(1,
Point<1>(0.5));
54 template <
int dim,
int spacedim>
69 std::vector<bool>(1, true)))
86 template <
int dim,
int spacedim>
93 polynomials.size() - 1,
96 polynomials.size() - 1),
98 polynomials.size() - 1)
104 polynomials.size() - 1)
106 std::vector<bool>(1, true)))
118 template <
int dim,
int spacedim>
126 std::ostringstream namebuf;
129 return namebuf.str();
134 template <
int dim,
int spacedim>
138 std::vector<double> & nodal_values)
const 149 nodal_values[i] = support_point_values[i](0);
155 template <
int dim,
int spacedim>
156 std::unique_ptr<FiniteElement<dim, spacedim>>
159 return std_cxx14::make_unique<FE_DGQ<dim, spacedim>>(*this);
168 template <
int dim,
int spacedim>
169 std::vector<unsigned int>
172 std::vector<unsigned int> dpo(dim + 1, 0U);
174 for (
unsigned int i = 1; i < dim; ++i)
181 template <
int dim,
int spacedim>
184 const char direction)
const 186 const unsigned int n = this->
degree + 1;
188 for (
unsigned int i = 1; i < dim; ++i)
197 for (
unsigned int i = n; i > 0;)
207 for (
unsigned int iz = 0; iz < ((dim > 2) ? n : 1); ++iz)
208 for (
unsigned int j = 0; j < n; ++j)
209 for (
unsigned int i = 0; i < n; ++i)
211 unsigned int k = n * i - j + n - 1 + n * n * iz;
218 for (
unsigned int iz = 0; iz < ((dim > 2) ? n : 1); ++iz)
219 for (
unsigned int iy = 0; iy < n; ++iy)
220 for (
unsigned int ix = 0; ix < n; ++ix)
222 unsigned int k = n * ix - iy + n - 1 + n * n * iz;
230 for (
unsigned int iz = 0; iz < n; ++iz)
231 for (
unsigned int iy = 0; iy < n; ++iy)
232 for (
unsigned int ix = 0; ix < n; ++ix)
234 unsigned int k = n * (n * iy - iz + n - 1) + ix;
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;
258 template <
int dim,
int spacedim>
270 typename FE::ExcInterpolationNotImplemented());
311 cell_interpolation.
mmult(interpolation_matrix, source_interpolation);
316 if (std::fabs(interpolation_matrix(i, j)) < 1e-15)
317 interpolation_matrix(i, j) = 0.;
328 sum += interpolation_matrix(i, j);
330 Assert(std::fabs(sum - 1) < 5e-14 * std::max(this->
degree, 1U) * dim,
337 template <
int dim,
int spacedim>
349 (void)interpolation_matrix;
353 typename FE::ExcInterpolationNotImplemented());
355 Assert(interpolation_matrix.
m() == 0,
357 Assert(interpolation_matrix.
n() == 0,
363 template <
int dim,
int spacedim>
376 (void)interpolation_matrix;
380 typename FE::ExcInterpolationNotImplemented());
382 Assert(interpolation_matrix.
m() == 0,
384 Assert(interpolation_matrix.
n() == 0,
390 template <
int dim,
int spacedim>
393 const unsigned int child,
402 "Prolongation matrices are only available for refined cells!"));
409 if (this->
prolongation[refinement_case - 1][child].n() == 0)
414 if (this->
prolongation[refinement_case - 1][child].n() ==
424 std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
426 isotropic_matrices.back().resize(
438 isotropic_matrices.back());
470 template <
int dim,
int spacedim>
473 const unsigned int child,
482 "Restriction matrices are only available for refined cells!"));
489 if (this->
restriction[refinement_case - 1][child].n() == 0)
494 if (this->
restriction[refinement_case - 1][child].n() ==
496 return this->
restriction[refinement_case - 1][child];
504 std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
506 isotropic_matrices.back().resize(
517 this_nonconst.
restriction[refinement_case - 1].swap(
518 isotropic_matrices.back());
545 return this->
restriction[refinement_case - 1][child];
550 template <
int dim,
int spacedim>
559 template <
int dim,
int spacedim>
560 std::vector<std::pair<unsigned int, unsigned int>>
568 return std::vector<std::pair<unsigned int, unsigned int>>();
573 template <
int dim,
int spacedim>
574 std::vector<std::pair<unsigned int, unsigned int>>
582 return std::vector<std::pair<unsigned int, unsigned int>>();
587 template <
int dim,
int spacedim>
588 std::vector<std::pair<unsigned int, unsigned int>>
596 return std::vector<std::pair<unsigned int, unsigned int>>();
601 template <
int dim,
int spacedim>
614 template <
int dim,
int spacedim>
617 const unsigned int face_index)
const 624 unsigned int n = this->
degree + 1;
635 bool support_points_on_boundary =
true;
636 for (
unsigned int d = 0; d < dim; ++d)
638 support_points_on_boundary =
false;
639 for (
unsigned int d = 0; d < dim; ++d)
641 support_points_on_boundary =
false;
642 if (support_points_on_boundary ==
false)
645 unsigned int n2 = n * n;
657 return (((shape_index == 0) && (face_index == 0)) ||
658 ((shape_index == this->
degree) && (face_index == 1)));
663 if (face_index == 0 && (shape_index % n) == 0)
665 if (face_index == 1 && (shape_index % n) == this->
degree)
667 if (face_index == 2 && shape_index < n)
669 if (face_index == 3 && shape_index >= this->dofs_per_cell - n)
676 const unsigned int in2 = shape_index % n2;
679 if (face_index == 0 && (shape_index % n) == 0)
682 if (face_index == 1 && (shape_index % n) == n - 1)
685 if (face_index == 2 && in2 < n)
688 if (face_index == 3 && in2 >= n2 - n)
691 if (face_index == 4 && shape_index < n2)
694 if (face_index == 5 && shape_index >= this->dofs_per_cell - n2)
707 template <
int dim,
int spacedim>
708 std::pair<Table<2, bool>, std::vector<unsigned int>>
712 constant_modes.
fill(
true);
713 return std::pair<Table<2, bool>, std::vector<unsigned int>>(
714 constant_modes, std::vector<unsigned int>(1, 0));
719 template <
int dim,
int spacedim>
731 template <
int dim,
int spacedim>
735 Polynomials::generate_complete_Lagrange_basis(points.get_points()))
744 template <
int dim,
int spacedim>
750 std::ostringstream namebuf;
751 bool equidistant =
true;
752 std::vector<double> points(this->
degree + 1);
754 std::vector<unsigned int> lexicographic =
756 for (
unsigned int j = 0; j <= this->
degree; j++)
760 for (
unsigned int j = 0; j <= this->
degree; j++)
761 if (std::abs(points[j] - (
double)j / this->
degree) > 1e-15)
766 if (this->
degree == 0 && std::abs(points[0] - 0.5) < 1e-15)
769 if (equidistant ==
true)
772 namebuf <<
"FE_DGQArbitraryNodes<" 774 <<
">(QIterated(QTrapez()," << this->
degree <<
"))";
777 << this->degree <<
")";
778 return namebuf.str();
783 bool gauss_lobatto =
true;
784 for (
unsigned int j = 0; j <= this->
degree; j++)
785 if (points[j] != points_gl.
point(j)(0))
787 gauss_lobatto =
false;
791 if (gauss_lobatto ==
true)
795 return namebuf.str();
801 for (
unsigned int j = 0; j <= this->
degree; j++)
802 if (points[j] != points_g.
point(j)(0))
811 <<
">(QGauss(" << this->
degree + 1 <<
"))";
812 return namebuf.str();
817 bool gauss_log =
true;
818 for (
unsigned int j = 0; j <= this->
degree; j++)
819 if (points[j] != points_glog.
point(j)(0))
825 if (gauss_log ==
true)
828 <<
">(QGaussLog(" << this->
degree + 1 <<
"))";
829 return namebuf.str();
834 <<
">(QUnknownNodes(" << this->
degree + 1 <<
"))";
835 return namebuf.str();
840 template <
int dim,
int spacedim>
845 std::vector<double> & nodal_values)
const 856 nodal_values[i] = support_point_values[i](0);
862 template <
int dim,
int spacedim>
863 std::unique_ptr<FiniteElement<dim, spacedim>>
867 std::vector<Point<1>> qpoints(this->
degree + 1);
868 std::vector<unsigned int> lexicographic =
870 for (
unsigned int i = 0; i <= this->
degree; ++i)
874 return std_cxx14::make_unique<FE_DGQArbitraryNodes<dim, spacedim>>(
882 template <
int dim,
int spacedim>
885 Polynomials::Legendre::generate_complete_basis(degree))
890 template <
int dim,
int spacedim>
891 std::pair<Table<2, bool>, std::vector<unsigned int>>
897 constant_modes(0, 0) =
true;
898 return std::pair<Table<2, bool>, std::vector<unsigned int>>(
899 constant_modes, std::vector<unsigned int>(1, 0));
904 template <
int dim,
int spacedim>
914 template <
int dim,
int spacedim>
915 std::unique_ptr<FiniteElement<dim, spacedim>>
918 return std_cxx14::make_unique<FE_DGQLegendre<dim, spacedim>>(this->
degree);
925 template <
int dim,
int spacedim>
928 Polynomials::HermiteLikeInterpolation::generate_complete_basis(degree))
933 template <
int dim,
int spacedim>
943 template <
int dim,
int spacedim>
944 std::unique_ptr<FiniteElement<dim, spacedim>>
947 return std_cxx14::make_unique<FE_DGQHermite<dim, spacedim>>(this->
degree);
953 #include "fe_dgq.inst" 956 DEAL_II_NAMESPACE_CLOSE
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const override
static::ExceptionBase & ExcFEHasNoSupportPoints()
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_quad_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
virtual const FullMatrix< double > & get_prolongation_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const override
const std::vector< Point< dim > > & get_unit_support_points() const
std::vector< std::vector< FullMatrix< double > > > restriction
#define AssertDimension(dim1, dim2)
virtual FiniteElementDomination::Domination compare_for_face_domination(const FiniteElement< dim, spacedim > &fe_other) const override
double compute_value(const unsigned int i, const Point< dim > &p) const
virtual bool hp_constraints_are_implemented() const override
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_vertex_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
FE_DGQLegendre(const unsigned int degree)
const unsigned int degree
virtual std::string get_name() const override
virtual void get_face_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix) const override
#define AssertThrow(cond, exc)
static::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
const std::vector< unsigned int > & get_numbering_inverse() const
const Point< dim > & point(const unsigned int i) const
std::vector< Point< dim > > unit_support_points
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
virtual std::size_t memory_consumption() const override
unsigned int size() const
static::ExceptionBase & ExcMessage(std::string arg1)
virtual std::string get_name() const override
std::vector< std::vector< FullMatrix< double > > > prolongation
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const override
static std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
#define Assert(cond, exc)
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
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 void get_subface_interpolation_matrix(const FiniteElement< dim, spacedim > &source, const unsigned int subface, FullMatrix< double > &matrix) const override
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
void mmult(FullMatrix< number2 > &C, const FullMatrix< number2 > &B, const bool adding=false) const
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
std::string dim_string(const int dim, const int spacedim)
const unsigned int dofs_per_cell
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const override
FE_DGQArbitraryNodes(const Quadrature< 1 > &points)
virtual void get_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix) const override
FE_DGQHermite(const unsigned int degree)
void rotate_indices(std::vector< unsigned int > &indices, const char direction) const
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
void reinit_restriction_and_prolongation_matrices(const bool isotropic_restriction_only=false, const bool isotropic_prolongation_only=false)
virtual std::string get_name() const override
virtual std::string get_name() const override
virtual const FullMatrix< double > & get_restriction_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const override
static::ExceptionBase & ExcNotImplemented()
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
TensorProductPolynomials< dim > poly_space
void fill(InputIterator entries, const bool C_style_indexing=true)
static::ExceptionBase & ExcInternalError()