17 #include <deal.II/base/quadrature.h> 18 #include <deal.II/base/quadrature_lib.h> 19 #include <deal.II/fe/fe.h> 20 #include <deal.II/fe/fe_dgq.h> 21 #include <deal.II/fe/fe_tools.h> 27 DEAL_II_NAMESPACE_OPEN
32 std::vector<Point<1> >
33 get_QGaussLobatto_points (
const unsigned int degree)
38 return std::vector<Point<1> >(1,
Point<1>(0.5));
44 template <
int dim,
int spacedim>
69 template <
int dim,
int spacedim>
88 template <
int dim,
int spacedim>
96 std::ostringstream namebuf;
99 <<
">(" << this->
degree <<
")";
100 return namebuf.str();
105 template <
int dim,
int spacedim>
118 template <
int dim,
int spacedim>
119 std::vector<unsigned int>
122 std::vector<unsigned int> dpo(dim+1, 0U);
124 for (
unsigned int i=1; i<dim; ++i)
131 template <
int dim,
int spacedim>
134 const char direction)
const 136 const unsigned int n = this->
degree+1;
138 for (
unsigned int i=1; i<dim; ++i)
147 for (
unsigned int i=n; i>0;)
157 for (
unsigned int iz=0; iz<((dim>2) ? n:1); ++iz)
158 for (
unsigned int j=0; j<n; ++j)
159 for (
unsigned int i=0; i<n; ++i)
161 unsigned int k = n*i-j+n-1 + n*n*iz;
168 for (
unsigned int iz=0; iz<((dim>2) ? n:1); ++iz)
169 for (
unsigned int iy=0; iy<n; ++iy)
170 for (
unsigned int ix=0; ix<n; ++ix)
172 unsigned int k = n*ix-iy+n-1 + n*n*iz;
180 for (
unsigned int iz=0; iz<n; ++iz)
181 for (
unsigned int iy=0; iy<n; ++iy)
182 for (
unsigned int ix=0; ix<n; ++ix)
184 unsigned int k = n*(n*iy-iz+n-1) + ix;
192 for (
unsigned int iz=0; iz<n; ++iz)
193 for (
unsigned int iy=0; iy<n; ++iy)
194 for (
unsigned int ix=0; ix<n; ++ix)
196 unsigned int k = n*(n*iy-iz+n-1) + ix;
208 template <
int dim,
int spacedim>
219 typename FE::ExcInterpolationNotImplemented() );
251 cell_interpolation(j,i)
255 source_interpolation(j,i)
264 cell_interpolation.
mmult (interpolation_matrix,
265 source_interpolation);
270 if (std::fabs(interpolation_matrix(i,j)) < 1e-15)
271 interpolation_matrix(i,j) = 0.;
282 sum += interpolation_matrix(i,j);
284 Assert (std::fabs(sum-1) < 5e-14*std::max(this->
degree,1U)*dim,
291 template <
int dim,
int spacedim>
303 (void)interpolation_matrix;
306 typename FE::ExcInterpolationNotImplemented());
308 Assert (interpolation_matrix.
m() == 0,
311 Assert (interpolation_matrix.
n() == 0,
318 template <
int dim,
int spacedim>
331 (void)interpolation_matrix;
334 typename FE::ExcInterpolationNotImplemented());
336 Assert (interpolation_matrix.
m() == 0,
339 Assert (interpolation_matrix.
n() == 0,
346 template <
int dim,
int spacedim>
355 ExcMessage(
"Prolongation matrices are only available for refined cells!"));
360 if (this->
prolongation[refinement_case-1][child].n() == 0)
374 std::vector<std::vector<FullMatrix<double> > >
376 isotropic_matrices.back().
383 isotropic_matrices,
true);
384 this_nonconst.
prolongation[refinement_case-1].swap(isotropic_matrices.back());
412 template <
int dim,
int spacedim>
421 ExcMessage(
"Restriction matrices are only available for refined cells!"));
426 if (this->
restriction[refinement_case-1][child].n() == 0)
431 if (this->
restriction[refinement_case-1][child].n() ==
433 return this->
restriction[refinement_case-1][child];
440 std::vector<std::vector<FullMatrix<double> > >
442 isotropic_matrices.back().
449 isotropic_matrices,
true);
450 this_nonconst.
restriction[refinement_case-1].swap(isotropic_matrices.back());
473 return this->
restriction[refinement_case-1][child];
478 template <
int dim,
int spacedim>
487 template <
int dim,
int spacedim>
488 std::vector<std::pair<unsigned int, unsigned int> >
497 std::vector<std::pair<unsigned int, unsigned int> > ();
502 template <
int dim,
int spacedim>
503 std::vector<std::pair<unsigned int, unsigned int> >
512 std::vector<std::pair<unsigned int, unsigned int> > ();
517 template <
int dim,
int spacedim>
518 std::vector<std::pair<unsigned int, unsigned int> >
527 std::vector<std::pair<unsigned int, unsigned int> > ();
532 template <
int dim,
int spacedim>
544 template <
int dim,
int spacedim>
547 const unsigned int face_index)
const 554 unsigned int n = this->
degree+1;
565 bool support_points_on_boundary =
true;
566 for (
unsigned int d=0; d<dim; ++d)
568 support_points_on_boundary =
false;
569 for (
unsigned int d=0; d<dim; ++d)
571 support_points_on_boundary =
false;
572 if (support_points_on_boundary ==
false)
575 unsigned int n2 = n*n;
587 return (((shape_index == 0) && (face_index == 0)) ||
588 ((shape_index == this->
degree) && (face_index == 1)));
593 if (face_index==0 && (shape_index % n) == 0)
595 if (face_index==1 && (shape_index % n) == this->
degree)
597 if (face_index==2 && shape_index < n)
599 if (face_index==3 && shape_index >= this->dofs_per_cell-n)
606 const unsigned int in2 = shape_index % n2;
609 if (face_index==0 && (shape_index % n) == 0)
612 if (face_index==1 && (shape_index % n) == n-1)
615 if (face_index==2 && in2 < n )
618 if (face_index==3 && in2 >= n2-n)
621 if (face_index==4 && shape_index < n2)
624 if (face_index==5 && shape_index >= this->dofs_per_cell - n2)
637 template <
int dim,
int spacedim>
638 std::pair<Table<2,bool>, std::vector<unsigned int> >
642 constant_modes.
fill(
true);
643 return std::pair<Table<2,bool>, std::vector<unsigned int> >
644 (constant_modes, std::vector<unsigned int>(1, 0));
650 template <
int dim,
int spacedim>
662 template <
int dim,
int spacedim>
664 :
FE_DGQ<dim,spacedim>(
Polynomials::generate_complete_Lagrange_basis(points.get_points()))
674 template <
int dim,
int spacedim>
680 std::ostringstream namebuf;
681 bool equidistant =
true;
682 std::vector<double> points(this->
degree+1);
685 for (
unsigned int j=0; j<=this->
degree; j++)
689 for (
unsigned int j=0; j<=this->
degree; j++)
690 if (std::abs(points[j] - (
double)j/this->
degree) > 1e-15)
695 if (this->
degree == 0 && std::abs(points[0]-0.5) < 1e-15)
698 if (equidistant ==
true)
704 return namebuf.str();
709 bool gauss_lobatto =
true;
710 for (
unsigned int j=0; j<=this->
degree; j++)
711 if (points[j] != points_gl.
point(j)(0))
713 gauss_lobatto =
false;
717 if (gauss_lobatto ==
true)
720 return namebuf.str();
726 for (
unsigned int j=0; j<=this->
degree; j++)
727 if (points[j] != points_g.
point(j)(0))
736 return namebuf.str();
741 bool gauss_log =
true;
742 for (
unsigned int j=0; j<=this->
degree; j++)
743 if (points[j] != points_glog.
point(j)(0))
749 if (gauss_log ==
true)
752 return namebuf.str();
757 return namebuf.str();
762 template <
int dim,
int spacedim>
767 std::vector<Point<1> > qpoints(this->
degree+1);
769 for (
unsigned int i=0; i<=this->
degree; ++i)
780 template <
int dim,
int spacedim>
787 template <
int dim,
int spacedim>
788 std::pair<Table<2,bool>, std::vector<unsigned int> >
794 constant_modes(0,0) =
true;
795 return std::pair<Table<2,bool>, std::vector<unsigned int> >
796 (constant_modes, std::vector<unsigned int>(1, 0));
801 template <
int dim,
int spacedim>
811 template <
int dim,
int spacedim>
822 template <
int dim,
int spacedim>
824 :
FE_DGQ<dim,spacedim>(degree < 3 ?
825 Polynomials::generate_complete_Lagrange_basis(get_QGaussLobatto_points(degree))
827 Polynomials::HermiteInterpolation::generate_complete_basis(degree))
832 template <
int dim,
int spacedim>
833 std::pair<Table<2,bool>, std::vector<unsigned int> >
845 for (
unsigned int i=0; i<(dim>2?2:1); ++i)
846 for (
unsigned int j=0; j<(dim>1?2:1); ++j)
847 for (
unsigned int k=0; k<2; ++k)
849 return std::pair<Table<2,bool>, std::vector<unsigned int> >
850 (constant_modes, std::vector<unsigned int>(1, 0));
856 template <
int dim,
int spacedim>
866 template <
int dim,
int spacedim>
876 #include "fe_dgq.inst" 879 DEAL_II_NAMESPACE_CLOSE
static::ExceptionBase & ExcFEHasNoSupportPoints()
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const
virtual FiniteElement< dim, spacedim > * clone() const
std::vector< std::vector< FullMatrix< double > > > restriction
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const
double compute_value(const unsigned int i, const Point< dim > &p) const
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const
virtual const FullMatrix< double > & get_prolongation_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const
FE_DGQLegendre(const unsigned int degree)
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const
virtual FiniteElement< dim, spacedim > * clone() const
const unsigned int degree
virtual std::string get_name() const
#define AssertThrow(cond, exc)
static::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
virtual void get_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix) const
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 FiniteElement< dim, spacedim > * clone() const
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_vertex_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const
const std::vector< Point< dim > > & get_points() const
unsigned int size() const
static::ExceptionBase & ExcMessage(std::string arg1)
virtual void get_subface_interpolation_matrix(const FiniteElement< dim, spacedim > &source, const unsigned int subface, FullMatrix< double > &matrix) const
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const
std::vector< std::vector< FullMatrix< double > > > prolongation
virtual std::string get_name() const
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 std::size_t memory_consumption() const
virtual std::string get_name() const
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)
std::string dim_string(const int dim, const int spacedim)
const unsigned int dofs_per_cell
FE_DGQArbitraryNodes(const Quadrature< 1 > &points)
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_quad_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const
virtual void get_face_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix) const
virtual FiniteElement< dim, spacedim > * clone() const
FE_DGQHermite(const unsigned int degree)
virtual bool hp_constraints_are_implemented() const
void rotate_indices(std::vector< unsigned int > &indices, const char direction) const
void reinit_restriction_and_prolongation_matrices(const bool isotropic_restriction_only=false, const bool isotropic_prolongation_only=false)
virtual std::string get_name() const
static::ExceptionBase & ExcNotImplemented()
virtual FiniteElementDomination::Domination compare_for_face_domination(const FiniteElement< dim, spacedim > &fe_other) const
virtual const FullMatrix< double > & get_restriction_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const
TensorProductPolynomials< dim > poly_space
void fill(InputIterator entries, const bool C_style_indexing=true)
static::ExceptionBase & ExcInternalError()