17 #include <deal.II/base/qprojector.h> 18 #include <deal.II/base/quadrature.h> 19 #include <deal.II/base/quadrature_lib.h> 20 #include <deal.II/base/std_cxx14/memory.h> 21 #include <deal.II/base/template_constraints.h> 23 #include <deal.II/dofs/dof_accessor.h> 24 #include <deal.II/dofs/dof_handler.h> 26 #include <deal.II/fe/fe_nothing.h> 27 #include <deal.II/fe/fe_q_bubbles.h> 28 #include <deal.II/fe/fe_tools.h> 29 #include <deal.II/fe/fe_values.h> 30 #include <deal.II/fe/mapping_q1.h> 32 #include <deal.II/grid/grid_generator.h> 33 #include <deal.II/grid/tria.h> 39 DEAL_II_NAMESPACE_OPEN
48 template <
int dim,
int spacedim>
51 const ::FE_Q_Bubbles<dim, spacedim> & fe,
53 const bool isotropic_only)
55 const unsigned int dpc = fe.dofs_per_cell;
56 const unsigned int degree = fe.degree;
59 std::unique_ptr<Quadrature<dim>> q_fine;
61 std::vector<double>(1, 1.));
66 q_fine = std_cxx14::make_unique<QGauss<dim>>(degree + 1);
67 else if (spacedim == 2)
68 q_fine = std_cxx14::make_unique<QAnisotropic<dim>>(
71 q_fine = std_cxx14::make_unique<QAnisotropic<dim>>(
76 q_fine = std_cxx14::make_unique<QGauss<dim>>(degree + 1);
78 q_fine = std_cxx14::make_unique<QAnisotropic<dim>>(
82 q_fine = std_cxx14::make_unique<QGauss<dim>>(degree + 1);
89 const unsigned int nq = q_fine->size();
92 unsigned int ref_case = (isotropic_only) ?
95 for (; ref_case <= RefinementCase<dim>::isotropic_refinement;
98 const unsigned int nc =
101 for (
unsigned int i = 0; i < nc; ++i)
103 Assert(matrices[ref_case - 1][i].n() == dpc,
106 Assert(matrices[ref_case - 1][i].m() == dpc,
118 dh.distribute_dofs(fe);
126 const unsigned int n_dofs = dh.n_dofs();
131 std::vector<std::vector<types::global_dof_index>> child_ldi(
132 nc, std::vector<types::global_dof_index>(fe.dofs_per_cell));
135 unsigned int child_no = 0;
136 typename ::DoFHandler<dim>::active_cell_iterator cell =
138 for (; cell != dh.end(); ++cell, ++child_no)
141 cell->get_dof_indices(child_ldi[child_no]);
143 for (
unsigned int q = 0; q < nq; ++q)
144 for (
unsigned int i = 0; i < dpc; ++i)
145 for (
unsigned int j = 0; j < dpc; ++j)
147 const unsigned int gdi = child_ldi[child_no][i];
148 const unsigned int gdj = child_ldi[child_no][j];
149 fine_mass(gdi, gdj) += fine.shape_value(i, q) *
150 fine.shape_value(j, q) *
153 for (
unsigned int k = 0; k < dim; ++k)
154 quad_tmp(k) = fine.quadrature_point(q)(k);
155 coarse_rhs_matrix(gdi, j) +=
156 fine.shape_value(i, q) * fe.shape_value(j, quad_tmp) *
163 fine_mass.gauss_jordan();
164 fine_mass.mmult(solution, coarse_rhs_matrix);
167 for (
unsigned int child_no = 0; child_no < nc; ++child_no)
168 for (
unsigned int i = 0; i < dpc; ++i)
169 for (
unsigned int j = 0; j < dpc; ++j)
171 const unsigned int gdi = child_ldi[child_no][i];
173 if (std::fabs(solution(gdi, j)) > 1.e-12)
174 matrices[ref_case - 1][child_no](i, j) = solution(gdi, j);
183 template <
int dim,
int spacedim>
193 get_riaf_vector(q_degree))
194 , n_bubbles((q_degree <= 1) ? 1 : dim)
197 ExcMessage(
"This element can only be used for polynomial degrees " 198 "greater than zero"));
204 for (
unsigned int d = 0; d < dim; ++d)
206 for (
unsigned int i = 0; i <
n_bubbles; ++i)
213 internal::FE_Q_Bubbles::compute_embedding_matrices(*
this,
223 template <
int dim,
int spacedim>
227 Polynomials::generate_complete_Lagrange_basis(points.get_points())),
233 ,
n_bubbles((points.size() - 1 <= 1) ? 1 : dim)
236 ExcMessage(
"This element can only be used for polynomial degrees " 243 for (
unsigned int d = 0; d < dim; ++d)
245 for (
unsigned int i = 0; i <
n_bubbles; ++i)
252 internal::FE_Q_Bubbles::compute_embedding_matrices(*
this,
262 template <
int dim,
int spacedim>
270 std::ostringstream namebuf;
272 const unsigned int n_points = this->
degree;
273 std::vector<double> points(n_points);
277 unsigned int index = 0;
282 if ((dim > 1) ? (unit_support_points[j](1) == 0 &&
283 ((dim > 2) ? unit_support_points[j](2) == 0 :
true)) :
287 points[index] = unit_support_points[j](0);
289 points[n_points - 1] = unit_support_points[j](0);
291 points[index - 1] = unit_support_points[j](0);
297 Assert(index == n_points || (dim == 1 && index == n_points +
n_bubbles),
299 "Could not decode support points in one coordinate direction."));
302 for (
unsigned int j = 0; j < n_points; j++)
303 if (std::fabs(points[j] - (
double)j / (this->degree - 1)) > 1e-15)
311 if (this->degree > 3)
313 <<
">(QIterated(QTrapez()," << this->degree - 1 <<
"))";
316 <<
">(" << this->degree - 1 <<
")";
323 for (
unsigned int j = 0; j < n_points; j++)
324 if (points[j] != points_gl.
point(j)(0))
330 namebuf <<
"FE_Q_Bubbles<" << dim <<
">(" << this->degree - 1 <<
")";
332 namebuf <<
"FE_Q_Bubbles<" << dim <<
">(QUnknownNodes(" << this->degree
335 return namebuf.str();
340 template <
int dim,
int spacedim>
341 std::unique_ptr<FiniteElement<dim, spacedim>>
344 return std_cxx14::make_unique<FE_Q_Bubbles<dim, spacedim>>(*this);
349 template <
int dim,
int spacedim>
354 std::vector<double> & nodal_values)
const 367 const std::pair<unsigned int, unsigned int> index =
369 nodal_values[i] = support_point_values[i](index.first);
373 for (
unsigned int i = 0; i <
n_bubbles; ++i)
374 nodal_values[nodal_values.size() - i - 1] = 0.;
379 template <
int dim,
int spacedim>
391 (x_source_fe.
get_name().find(
"FE_Q_Bubbles<") == 0) ||
392 (dynamic_cast<const FEQBUBBLES *>(&x_source_fe) !=
nullptr),
401 auto casted_fe =
dynamic_cast<const FEQBUBBLES *
>(&x_source_fe);
402 if (casted_fe !=
nullptr && casted_fe->degree == this->degree)
403 for (
unsigned int i = 0; i < interpolation_matrix.
m(); ++i)
404 interpolation_matrix.
set(i, i, 1.);
415 template <
int dim,
int spacedim>
419 unsigned int n_cont_dofs = Utilities::fixed_power<dim>(q_deg + 1);
420 const unsigned int n_bubbles = (q_deg <= 1 ? 1 : dim);
421 return std::vector<bool>(n_cont_dofs +
n_bubbles,
true);
426 template <
int dim,
int spacedim>
427 std::vector<unsigned int>
430 std::vector<unsigned int> dpo(dim + 1, 1U);
431 for (
unsigned int i = 1; i < dpo.size(); ++i)
432 dpo[i] = dpo[i - 1] * (q_deg - 1);
435 (q_deg <= 1 ? 1 : dim);
441 template <
int dim,
int spacedim>
444 const unsigned int shape_index,
445 const unsigned int face_index)
const 457 template <
int dim,
int spacedim>
460 const unsigned int child,
469 "Prolongation matrices are only available for refined cells!"));
476 ExcMessage(
"This prolongation matrix has not been computed yet!"));
483 template <
int dim,
int spacedim>
486 const unsigned int child,
495 "Restriction matrices are only available for refined cells!"));
502 ExcMessage(
"This restriction matrix has not been computed yet!"));
505 return this->
restriction[refinement_case - 1][child];
509 #include "fe_q_bubbles.inst" 511 DEAL_II_NAMESPACE_CLOSE
Transformed quadrature weights.
const unsigned int n_bubbles
virtual std::string get_name() const override
virtual const FullMatrix< double > & get_prolongation_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case) const override
std::vector< std::vector< FullMatrix< double > > > restriction
#define AssertDimension(dim1, dim2)
const unsigned int degree
static std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
virtual const FullMatrix< double > & get_restriction_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case) const override
Transformed quadrature points.
#define AssertThrow(cond, exc)
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const override
static::ExceptionBase & ExcInterpolationNotImplemented()
static::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
const Point< dim > & point(const unsigned int i) const
std::vector< Point< dim > > unit_support_points
void set(const size_type i, const size_type j, const number value)
virtual void get_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix) const override
virtual void execute_coarsening_and_refinement()
const std::vector< Point< dim > > & get_points() const
unsigned int size() const
static::ExceptionBase & ExcMessage(std::string arg1)
std::vector< std::vector< FullMatrix< double > > > prolongation
#define Assert(cond, exc)
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
virtual std::string get_name() const =0
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
unsigned int n_components() const
std::string dim_string(const int dim, const int spacedim)
const unsigned int dofs_per_cell
static unsigned int n_children(const RefinementCase< dim > &refinement_case)
void initialize(const std::vector< Point< 1 >> &support_points_1d)
std::pair< unsigned int, unsigned int > system_to_component_index(const unsigned int index) const
static std::vector< bool > get_riaf_vector(const unsigned int degree)
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
unsigned int n_dofs_per_cell() const
void reinit_restriction_and_prolongation_matrices(const bool isotropic_restriction_only=false, const bool isotropic_prolongation_only=false)
FE_Q_Bubbles(const unsigned int p)
active_cell_iterator begin_active(const unsigned int level=0) const
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
static::ExceptionBase & ExcInternalError()