52 get_QGaussLobatto_points(
const unsigned int degree)
57 return std::vector<Point<1>>(1,
Point<1>(0.5));
65 template <
int dim,
int spacedim>
99 template <
int dim,
int spacedim>
106 polynomials.size() - 1,
109 polynomials.size() - 1),
111 polynomials.size() - 1)
117 polynomials.size() - 1)
131 template <
int dim,
int spacedim>
139 std::ostringstream namebuf;
141 << this->degree <<
")";
142 return namebuf.str();
147 template <
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);
168 template <
int dim,
int spacedim>
169 std::unique_ptr<FiniteElement<dim, spacedim>>
172 return std::make_unique<FE_DGQ<dim, spacedim>>(*this);
181 template <
int dim,
int spacedim>
182 std::vector<unsigned int>
185 std::vector<unsigned int> dpo(dim + 1, 0
U);
187 for (
unsigned int i = 1; i < dim; ++i)
194 template <
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;
271 template <
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.;
339 for (
unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
343 sum += interpolation_matrix(i, j);
353 template <
int dim,
int spacedim>
358 const unsigned int)
const
366 (void)interpolation_matrix;
370 typename FE::ExcInterpolationNotImplemented());
372 Assert(interpolation_matrix.
m() == 0,
374 Assert(interpolation_matrix.
n() == 0,
380 template <
int dim,
int spacedim>
386 const unsigned int)
const
394 (void)interpolation_matrix;
398 typename FE::ExcInterpolationNotImplemented());
400 Assert(interpolation_matrix.
m() == 0,
402 Assert(interpolation_matrix.
n() == 0,
408 template <
int dim,
int spacedim>
411 const unsigned int child,
418 "Prolongation matrices are only available for refined cells!"));
422 if (this->prolongation[refinement_case - 1][child].n() == 0)
424 std::lock_guard<std::mutex> lock(this->mutex);
427 if (this->prolongation[refinement_case - 1][child].n() ==
428 this->n_dofs_per_cell())
429 return this->prolongation[refinement_case - 1][child];
437 std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
439 isotropic_matrices.back().resize(
442 this->n_dofs_per_cell()));
452 isotropic_matrices.back());
479 return this->prolongation[refinement_case - 1][child];
484 template <
int dim,
int spacedim>
487 const unsigned int child,
494 "Restriction matrices are only available for refined cells!"));
498 if (this->restriction[refinement_case - 1][child].n() == 0)
500 std::lock_guard<std::mutex> lock(this->mutex);
503 if (this->restriction[refinement_case - 1][child].n() ==
504 this->n_dofs_per_cell())
505 return this->restriction[refinement_case - 1][child];
513 std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
515 isotropic_matrices.back().resize(
518 this->n_dofs_per_cell()));
527 this_nonconst.
restriction[refinement_case - 1].swap(
528 isotropic_matrices.back());
555 return this->restriction[refinement_case - 1][child];
560 template <
int dim,
int spacedim>
569 template <
int dim,
int spacedim>
570 std::vector<std::pair<unsigned int, unsigned int>>
578 return std::vector<std::pair<unsigned int, unsigned int>>();
583 template <
int dim,
int spacedim>
584 std::vector<std::pair<unsigned int, unsigned int>>
592 return std::vector<std::pair<unsigned int, unsigned int>>();
597 template <
int dim,
int spacedim>
598 std::vector<std::pair<unsigned int, unsigned int>>
601 const unsigned int)
const
607 return std::vector<std::pair<unsigned int, unsigned int>>();
612 template <
int dim,
int spacedim>
616 const unsigned int codim)
const
636 if (this->degree < fe_dgq_other->degree)
638 else if (this->degree == fe_dgq_other->degree)
646 if (this->degree < fe_q_other->degree)
648 else if (this->degree == fe_q_other->degree)
656 if (this->degree < fe_bernstein_other->degree)
658 else if (this->degree == fe_bernstein_other->degree)
666 if (this->degree < fe_bubbles_other->degree)
668 else if (this->degree == fe_bubbles_other->degree)
676 if (this->degree < fe_dg0_other->degree)
678 else if (this->degree == fe_dg0_other->degree)
686 if (this->degree < fe_q_iso_q1_other->degree)
688 else if (this->degree == fe_q_iso_q1_other->degree)
696 if (this->degree < fe_hierarchical_other->degree)
698 else if (this->degree == fe_hierarchical_other->degree)
706 if (this->degree < fe_dgp_other->degree)
708 else if (this->degree == fe_dgp_other->degree)
716 if (fe_nothing->is_dominating())
731 template <
int dim,
int spacedim>
734 const unsigned int face_index)
const
739 unsigned int n = this->degree + 1;
750 bool support_points_on_boundary =
true;
751 for (
unsigned int d = 0;
d < dim; ++
d)
753 support_points_on_boundary =
false;
754 for (
unsigned int d = 0;
d < dim; ++
d)
756 support_points_on_boundary =
false;
757 if (support_points_on_boundary ==
false)
760 unsigned int n2 = n * n;
772 return (((shape_index == 0) && (face_index == 0)) ||
773 ((shape_index == this->degree) && (face_index == 1)));
778 if (face_index == 0 && (shape_index % n) == 0)
780 if (face_index == 1 && (shape_index % n) == this->degree)
782 if (face_index == 2 && shape_index < n)
784 if (face_index == 3 && shape_index >= this->n_dofs_per_cell() - n)
791 const unsigned int in2 = shape_index % n2;
794 if (face_index == 0 && (shape_index % n) == 0)
797 if (face_index == 1 && (shape_index % n) == n - 1)
800 if (face_index == 2 && in2 < n)
803 if (face_index == 3 && in2 >= n2 - n)
806 if (face_index == 4 && shape_index < n2)
809 if (face_index == 5 && shape_index >= this->n_dofs_per_cell() - n2)
822 template <
int dim,
int spacedim>
823 std::pair<Table<2, bool>, std::vector<unsigned int>>
827 constant_modes.fill(
true);
828 return std::pair<Table<2, bool>, std::vector<unsigned int>>(
829 constant_modes, std::vector<unsigned int>(1, 0));
836 template <
int dim,
int spacedim>
849 template <
int dim,
int spacedim>
855 std::ostringstream namebuf;
856 bool equidistant =
true;
857 std::vector<double> points(this->degree + 1);
859 auto *
const polynomial_space =
862 std::vector<unsigned int> lexicographic =
863 polynomial_space->get_numbering_inverse();
864 for (
unsigned int j = 0; j <= this->degree; ++j)
868 for (
unsigned int j = 0; j <= this->degree; ++j)
869 if (std::abs(points[j] -
static_cast<double>(j) / this->degree) > 1
e-15)
874 if (this->degree == 0 && std::abs(points[0] - 0.5) < 1
e-15)
877 if (equidistant ==
true)
879 if (this->degree > 2)
880 namebuf <<
"FE_DGQArbitraryNodes<"
882 <<
">(QIterated(QTrapezoid()," << this->degree <<
"))";
885 << this->degree <<
")";
886 return namebuf.str();
891 bool gauss_lobatto =
true;
892 for (
unsigned int j = 0; j <= this->degree; ++j)
893 if (points[j] != points_gl.
point(j)(0))
895 gauss_lobatto =
false;
899 if (gauss_lobatto ==
true)
902 << this->degree <<
")";
903 return namebuf.str();
907 const QGauss<1> points_g(this->degree + 1);
909 for (
unsigned int j = 0; j <= this->degree; ++j)
910 if (points[j] != points_g.
point(j)(0))
919 <<
">(QGauss(" << this->degree + 1 <<
"))";
920 return namebuf.str();
925 bool gauss_log =
true;
926 for (
unsigned int j = 0; j <= this->degree; ++j)
927 if (points[j] != points_glog.
point(j)(0))
933 if (gauss_log ==
true)
936 <<
">(QGaussLog(" << this->degree + 1 <<
"))";
937 return namebuf.str();
942 <<
">(QUnknownNodes(" << this->degree + 1 <<
"))";
943 return namebuf.str();
948 template <
int dim,
int spacedim>
953 std::vector<double> &nodal_values)
const
956 this->get_unit_support_points().size());
960 for (
unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
964 nodal_values[i] = support_point_values[i](0);
970 template <
int dim,
int spacedim>
971 std::unique_ptr<FiniteElement<dim, spacedim>>
975 std::vector<Point<1>> qpoints(this->degree + 1);
976 auto *
const polynomial_space =
979 std::vector<unsigned int> lexicographic =
980 polynomial_space->get_numbering_inverse();
981 for (
unsigned int i = 0; i <= this->degree; ++i)
985 return std::make_unique<FE_DGQArbitraryNodes<dim, spacedim>>(pquadrature);
992 template <
int dim,
int spacedim>
995 Polynomials::Legendre::generate_complete_basis(degree))
1000 template <
int dim,
int spacedim>
1001 std::pair<Table<2, bool>, std::vector<unsigned int>>
1007 constant_modes(0, 0) =
true;
1008 return std::pair<Table<2, bool>, std::vector<unsigned int>>(
1009 constant_modes, std::vector<unsigned int>(1, 0));
1014 template <
int dim,
int spacedim>
1024 template <
int dim,
int spacedim>
1025 std::unique_ptr<FiniteElement<dim, spacedim>>
1028 return std::make_unique<FE_DGQLegendre<dim, spacedim>>(this->degree);
1035 template <
int dim,
int spacedim>
1038 Polynomials::HermiteLikeInterpolation::generate_complete_basis(degree))
1043 template <
int dim,
int spacedim>
1053 template <
int dim,
int spacedim>
1054 std::unique_ptr<FiniteElement<dim, spacedim>>
1057 return std::make_unique<FE_DGQHermite<dim, spacedim>>(this->degree);
1063 #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 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::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 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 std::vector< Point< dim > > & get_points() const
const Point< dim > & point(const unsigned int i) 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 & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
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 > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
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)
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)