<ol>
+ <li> New: FE_Q_Bubbles describes a FiniteElement based on FE_Q
+ enriched by bubble functions.
+ <br>
+ (Daniel Arndt, 2015/08/12)
+ </li>
<li> New: MultithreadInfo::set_thread_limit() can now be called more than
once and the environment variable DEAL_II_NUM_THREADS will be respected
#include <deal.II/base/polynomial_space.h>
#include <deal.II/base/tensor_product_polynomials.h>
#include <deal.II/base/tensor_product_polynomials_const.h>
+#include <deal.II/base/tensor_product_polynomials_bubbles.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/fe_poly.h>
/*@{*/
/**
- * This class collects the basic methods used in FE_Q and FE_Q_DG0. There is
- * no public constructor for this class as it is not functional as a stand-
+ * This class collects the basic methods used in FE_Q, FE_Q_DG0 and FE_Q_Bubbles.
+ * There is no public constructor for this class as it is not functional as a stand-
* alone. The completion of definitions is left to the derived classes.
*
* @author Wolfgang Bangerth, 1998, 2003; Guido Kanschat, 2001; Ralf Hartmann,
* Mutex for protecting initialization of restriction and embedding matrix.
*/
mutable Threads::Mutex mutex;
+
+ /*
+ * The highest polynomial degree of the underlying tensor product space
+ * without any enrichment. For FE_Q*(p) this is p. Note that enrichments
+ * may lead to a difference to degree.
+ */
+ const unsigned int q_degree;
};
--- /dev/null
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2012 - 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+#ifndef __deal2__fe_q_bubbles_h
+#define __deal2__fe_q_bubbles_h
+
+#include <deal.II/base/config.h>
+#include <deal.II/base/tensor_product_polynomials_bubbles.h>
+#include <deal.II/fe/fe_q_base.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+
+/*!@addtogroup fe */
+/*@{*/
+
+/**
+ * Implementation of a scalar Lagrange finite element @p Q_p^+ that yields the
+ * finite element space of continuous, piecewise polynomials of degree @p p in
+ * each coordinate direction plus some bubble enrichment space spanned by
+ * $(2x_j-1)^{p-1}\prod_{i=0}^{dim-1}(x_i(1-x_i))$. Therefore the highest polynomial
+ * degree is $p+1$.
+ * This class is realized using tensor product polynomials based on equidistant
+ * or given support points.
+ *
+ * The standard constructor of this class takes the degree @p p of this finite
+ * element. Alternatively, it can take a quadrature formula @p points defining
+ * the support points of the Lagrange interpolation in one coordinate direction.
+ *
+ * For more information about the <tt>spacedim</tt> template parameter check
+ * the documentation of FiniteElement or the one of Triangulation.
+ *
+ * Due to the fact that the enrichments are small almost everywhere
+ * for large p, the condition number for the mass and stiffness matrix fastly
+ * increaseses with increasing p.
+ * Below you see a comparison with FE_Q(QGaussLobatto(p+1)) for dim=1.
+ *
+ * <p ALIGN="center">
+ * @image html fe_q_bubbles_conditioning.png
+ * </p>
+ *
+ * Therefore, this element should be used with care for $p>3$.
+ *
+ * <h3>Implementation</h3>
+ *
+ * The constructor creates a TensorProductPolynomials object that includes the
+ * tensor product of @p LagrangeEquidistant polynomials of degree @p p plus the
+ * bubble enrichments. This @p TensorProductPolynomialsBubbles object
+ * provides all values and derivatives of the shape functions. In case a
+ * quadrature rule is given, the constructor creates a
+ * TensorProductPolynomialsBubbles object that includes the tensor product of @p
+ * Lagrange polynomials with the support points from @p points and the bubble enrichments
+ * as defined above.
+ *
+ * Furthermore the constructor fills the @p interface_constrains, the @p
+ * prolongation (embedding) and the @p restriction matrices.
+ *
+ * <h3>Numbering of the degrees of freedom (DoFs)</h3>
+ *
+ * The original ordering of the shape functions represented by the
+ * TensorProductPolynomialsBubbles is a tensor product
+ * numbering. However, the shape functions on a cell are renumbered
+ * beginning with the shape functions whose support points are at the
+ * vertices, then on the line, on the quads, and finally (for 3d) on
+ * the hexes. Finally, there are support points for the bubble enrichments
+ * in the middle of the cell.
+ *
+ */
+template <int dim, int spacedim=dim>
+class FE_Q_Bubbles : public FE_Q_Base<TensorProductPolynomialsBubbles<dim>,dim,spacedim>
+{
+public:
+ /**
+ * Constructor for tensor product polynomials of degree @p p plus bubble enrichments
+ *
+ */
+ FE_Q_Bubbles (const unsigned int p);
+
+ /**
+ * Constructor for tensor product polynomials with support points @p points
+ * plus bubble enrichments based on a one-dimensional quadrature
+ * formula. The degree of the finite element is <tt>points.size()</tt>.
+ * Note that the first point has to be 0 and the last one 1.
+ */
+ FE_Q_Bubbles (const Quadrature<1> &points);
+
+ /**
+ * Return a string that uniquely identifies a finite element. This class
+ * returns <tt>FE_Q_Bubbles<dim>(degree)</tt>, with @p dim and @p degree
+ * replaced by appropriate values.
+ */
+ virtual std::string get_name () const;
+
+ /**
+ * Interpolate a set of scalar values, computed in the generalized support
+ * points.
+ */
+ virtual void interpolate(std::vector<double> &local_dofs,
+ const std::vector<double> &values) const;
+
+ /**
+ * Interpolate a set of vector values, computed in the generalized support
+ * points.
+ *
+ * Since a finite element often only interpolates part of a vector,
+ * <tt>offset</tt> is used to determine the first component of the vector to
+ * be interpolated. Maybe consider changing your data structures to use the
+ * next function.
+ */
+ virtual void interpolate(std::vector<double> &local_dofs,
+ const std::vector<Vector<double> > &values,
+ unsigned int offset = 0) const;
+
+ /**
+ * Interpolate a set of vector values, computed in the generalized support
+ * points.
+ */
+ virtual void interpolate(
+ std::vector<double> &local_dofs,
+ const VectorSlice<const std::vector<std::vector<double> > > &values) const;
+
+ /**
+ * Return the matrix interpolating from the given finite element to the
+ * present one. The size of the matrix is then @p dofs_per_cell times
+ * <tt>source.dofs_per_cell</tt>.
+ *
+ * These matrices are only available if the source element is also a @p
+ * FE_Q_Bubbles element. Otherwise, an exception of type
+ * FiniteElement<dim,spacedim>::ExcInterpolationNotImplemented is thrown.
+ */
+ virtual void
+ get_interpolation_matrix (const FiniteElement<dim,spacedim> &source,
+ FullMatrix<double> &matrix) const;
+
+ virtual const FullMatrix<double> &
+ get_prolongation_matrix (const unsigned int child,
+ const RefinementCase<dim> &refinement_case) const;
+
+ virtual const FullMatrix<double> &
+ get_restriction_matrix (const unsigned int child,
+ const RefinementCase<dim> &refinement_case) const;
+
+ /**
+ * Check for non-zero values on a face.
+ *
+ * This function returns @p true, if the shape function @p shape_index has
+ * non-zero values on the face @p face_index.
+ *
+ * Implementation of the interface in FiniteElement
+ */
+ virtual bool has_support_on_face (const unsigned int shape_index,
+ const unsigned int face_index) const;
+
+protected:
+ /**
+ * @p clone function instead of a copy constructor.
+ *
+ * This function is needed by the constructors of @p FESystem.
+ */
+ virtual FiniteElement<dim,spacedim> *clone() const;
+
+private:
+
+ /**
+ * Returns the restriction_is_additive flags.
+ * Only the last components for the bubble enrichments are true.
+ */
+ static std::vector<bool> get_riaf_vector(const unsigned int degree);
+
+ /**
+ * Only for internal use. Its full name is @p get_dofs_per_object_vector
+ * function and it creates the @p dofs_per_object vector that is needed
+ * within the constructor to be passed to the constructor of @p
+ * FiniteElementData.
+ */
+ static std::vector<unsigned int> get_dpo_vector(const unsigned int degree);
+
+ /**
+ * Number of additional bubble functions
+ */
+ const unsigned int n_bubbles;
+};
+
+
+
+/*@}*/
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
fe_poly_tensor.cc
fe_q_base.cc
fe_q.cc
+ fe_q_bubbles.cc
fe_q_dg0.cc
fe_q_hierarchical.cc
fe_q_iso_q1.cc
fe_poly.inst.in
fe_poly_tensor.inst.in
fe_q_base.inst.in
+ fe_q_bubbles.inst.in
fe_q_dg0.inst.in
fe_q_hierarchical.inst.in
fe_q.inst.in
#include <deal.II/base/qprojector.h>
#include <deal.II/base/tensor_product_polynomials.h>
#include <deal.II/base/tensor_product_polynomials_const.h>
+#include <deal.II/base/tensor_product_polynomials_bubbles.h>
#include <deal.II/base/polynomials_p.h>
#include <deal.II/base/polynomials_piecewise.h>
#include <deal.II/fe/fe_poly.h>
// ---------------------------------------------------------------------
//
-// Copyright (C) 1998 - 2014 by the deal.II authors
+// Copyright (C) 1998 - 2015 by the deal.II authors
//
// This file is part of the deal.II library.
//
#if deal_II_dimension <= deal_II_space_dimension
template class FE_Poly<TensorProductPolynomials<deal_II_dimension>, deal_II_dimension, deal_II_space_dimension>;
template class FE_Poly<TensorProductPolynomialsConst<deal_II_dimension>, deal_II_dimension, deal_II_space_dimension>;
+ template class FE_Poly<TensorProductPolynomialsBubbles<deal_II_dimension>, deal_II_dimension, deal_II_space_dimension>;
template class FE_Poly<TensorProductPolynomials<deal_II_dimension,Polynomials::PiecewisePolynomial<double> >, deal_II_dimension, deal_II_space_dimension>;
template class FE_Poly<PolynomialSpace<deal_II_dimension>, deal_II_dimension, deal_II_space_dimension>;
template class FE_Poly<PolynomialsP<deal_II_dimension>, deal_II_dimension, deal_II_space_dimension>;
#include <deal.II/base/template_constraints.h>
#include <deal.II/base/tensor_product_polynomials.h>
#include <deal.II/base/tensor_product_polynomials_const.h>
+#include <deal.II/base/tensor_product_polynomials_bubbles.h>
#include <deal.II/base/polynomials_piecewise.h>
#include <deal.II/fe/fe_q_base.h>
#include <deal.II/fe/fe_nothing.h>
{
const unsigned int dim = 2;
+ unsigned int q_deg = fe.degree;
+ if (types_are_equal<POLY, TensorProductPolynomialsBubbles<dim> >::value)
+ q_deg = fe.degree-1;
+
// restricted to each face, the traces of the shape functions is an
// element of P_{k} (in 2d), or Q_{k} (in 3d), where k is the degree of
// the element. from this, we interpolate between mother and cell face.
// Add midpoint
constraint_points.push_back (Point<dim-1> (0.5));
- if (fe.degree>1)
+ if (q_deg>1)
{
- const unsigned int n=fe.degree-1;
- const double step=1./fe.degree;
+ const unsigned int n=q_deg-1;
+ const double step=1./q_deg;
// subface 0
for (unsigned int i=1; i<=n; ++i)
constraint_points.push_back (
const std::vector<unsigned int> &index_map_inverse =
fe.poly_space.get_numbering_inverse();
const std::vector<unsigned int> face_index_map =
- FE_Q_Helper::face_lexicographic_to_hierarchic_numbering<dim>(fe.degree);
+ FE_Q_Helper::face_lexicographic_to_hierarchic_numbering<dim>(q_deg);
Assert(std::abs(fe.poly_space.compute_value(index_map_inverse[0],Point<dim>())
- 1.) < 1e-14,
ExcInternalError());
for (unsigned int i=0; i<constraint_points.size(); ++i)
- for (unsigned int j=0; j<fe.degree+1; ++j)
+ for (unsigned int j=0; j<q_deg+1; ++j)
{
Point<dim> p;
p[0] = constraint_points[i](0);
{
const unsigned int dim = 3;
+ unsigned int q_deg = fe.degree;
+ if (types_are_equal<POLY,TensorProductPolynomialsBubbles<dim> >::value)
+ q_deg = fe.degree-1;
+
// For a detailed documentation of the interpolation see the
// FE_Q_Base<2>::initialize_constraints function.
constraint_points.push_back (Point<dim-1> (0.5, 0));
constraint_points.push_back (Point<dim-1> (0.5, 1));
- if (fe.degree>1)
+ if (q_deg>1)
{
- const unsigned int n=fe.degree-1;
- const double step=1./fe.degree;
+ const unsigned int n=q_deg-1;
+ const double step=1./q_deg;
std::vector<Point<dim-2> > line_support_points(n);
for (unsigned int i=0; i<n; ++i)
line_support_points[i](0)=(i+1)*step;
// Now construct relation between destination (child) and source (mother)
// dofs.
- const unsigned int pnts=(fe.degree+1)*(fe.degree+1);
+ const unsigned int pnts=(q_deg+1)*(q_deg+1);
// use that the element evaluates to 1 at index 0 and along the line at
// zero
const std::vector<unsigned int> &index_map_inverse =
fe.poly_space.get_numbering_inverse();
const std::vector<unsigned int> face_index_map =
- FE_Q_Helper::face_lexicographic_to_hierarchic_numbering<dim>(fe.degree);
+ FE_Q_Helper::face_lexicographic_to_hierarchic_numbering<dim>(q_deg);
Assert(std::abs(fe.poly_space.compute_value(index_map_inverse[0],Point<dim>())
- 1.) < 1e-14,
ExcInternalError());
for (unsigned int i=0; i<constraint_points.size(); ++i)
{
- const double interval = (double) (fe.degree * 2);
+ const double interval = (double) (q_deg * 2);
bool mirror[dim - 1];
Point<dim> constraint_point;
for (unsigned int j=0; j<pnts; ++j)
{
- unsigned int indices[2] = { j % (fe.degree+1), j / (fe.degree+1) };
+ unsigned int indices[2] = { j % (q_deg+1), j / (q_deg+1) };
for (unsigned int k = 0; k<2; ++k)
if (mirror[k])
- indices[k] = fe.degree - indices[k];
+ indices[k] = q_deg - indices[k];
const unsigned int
- new_index = indices[1] * (fe.degree + 1) + indices[0];
+ new_index = indices[1] * (q_deg + 1) + indices[0];
fe.interface_constraints(i,face_index_map[j]) =
fe.poly_space.compute_value (index_map_inverse[new_index],
const std::vector<bool> &restriction_is_additive_flags)
:
FE_Poly<POLY,dim,spacedim>(poly_space, fe_data, restriction_is_additive_flags,
- std::vector<ComponentMask>(1, std::vector<bool>(1,true)))
+ std::vector<ComponentMask>(1, std::vector<bool>(1,true))),
+ q_degree (types_are_equal<POLY, TensorProductPolynomialsBubbles<dim> >::value
+ ?this->degree-1
+ :this->degree)
{}
// distinguish q/q_dg0 case: need to be flexible enough to allow more
// degrees of freedom than there are FE_Q degrees of freedom for derived
// class FE_Q_DG0 that otherwise shares 95% of the code.
- const unsigned int q_dofs_per_cell = Utilities::fixed_power<dim>(this->degree+1);
+ const unsigned int q_dofs_per_cell = Utilities::fixed_power<dim>(q_degree+1);
Assert(q_dofs_per_cell == this->dofs_per_cell ||
- q_dofs_per_cell+1 == this->dofs_per_cell, ExcInternalError());
+ q_dofs_per_cell+1 == this->dofs_per_cell ||
+ q_dofs_per_cell+dim == this->dofs_per_cell , ExcInternalError());
{
std::vector<unsigned int> renumber(q_dofs_per_cell);
- const FiniteElementData<dim> fe(get_dpo_vector(this->degree),1,
- this->degree);
+ const FiniteElementData<dim> fe(get_dpo_vector(q_degree),1,
+ q_degree);
FETools::hierarchic_to_lexicographic_numbering (fe, renumber);
- if (this->dofs_per_cell > q_dofs_per_cell)
- renumber.push_back(q_dofs_per_cell);
+ for (unsigned int i= q_dofs_per_cell; i<this->dofs_per_cell; ++i)
+ renumber.push_back(i);
this->poly_space.set_numbering(renumber);
}
x_source_fe.dofs_per_cell));
// only evaluate Q dofs
- const unsigned int q_dofs_per_cell = Utilities::fixed_power<dim>(this->degree+1);
+ const unsigned int q_dofs_per_cell = Utilities::fixed_power<dim>(q_degree+1);
const unsigned int source_q_dofs_per_cell = Utilities::fixed_power<dim>(source_fe->degree+1);
// evaluation is simply done by evaluating the other FE's basis functions on
}
// cut off very small values
- const double eps = 2e-13*this->degree*dim;
+ const double eps = 2e-13*q_degree*dim;
for (unsigned int i=0; i<this->dofs_per_cell; ++i)
for (unsigned int j=0; j<source_fe->dofs_per_cell; ++j)
if (std::fabs(interpolation_matrix(i,j)) < eps)
// Rule of thumb for FP accuracy, that can be expected for a given
// polynomial degree. This value is used to cut off values close to
// zero.
- double eps = 2e-13*this->degree*(dim-1);
+ double eps = 2e-13*q_degree*(dim-1);
// compute the interpolation matrix by simply taking the value at the
// support points.
return;
const unsigned int codim = dim-1;
- this->unit_face_support_points.resize(Utilities::fixed_power<codim>(this->degree+1));
+ this->unit_face_support_points.resize(Utilities::fixed_power<codim>(q_degree+1));
// find renumbering of faces and assign from values of quadrature
std::vector<unsigned int> face_index_map =
- FE_Q_Helper::face_lexicographic_to_hierarchic_numbering<dim>(this->degree);
+ FE_Q_Helper::face_lexicographic_to_hierarchic_numbering<dim>(q_degree);
Quadrature<1> support_1d(points);
Quadrature<codim> support_quadrature(support_1d);
this->unit_face_support_points.resize(support_quadrature.size());
Assert (this->adjust_quad_dof_index_for_face_orientation_table.n_elements()==8*this->dofs_per_quad,
ExcInternalError());
- const unsigned int n=this->degree-1;
+ const unsigned int n=q_degree-1;
Assert(n*n==this->dofs_per_quad, ExcInternalError());
// alias for the table to fill
return this->prolongation[refinement_case-1][child];
// distinguish q/q_dg0 case: only treat Q dofs first
- const unsigned int q_dofs_per_cell = Utilities::fixed_power<dim>(this->degree+1);
+ const unsigned int q_dofs_per_cell = Utilities::fixed_power<dim>(q_degree+1);
// compute the interpolation matrices in much the same way as we do for
// the constraints. it's actually simpler here, since we don't have this
// interpolation matrix is formed by a permutation of the indices of the
// cell matrix. The value eps is used a threshold to decide when certain
// evaluations of the Lagrange polynomials are zero or one.
- const double eps = 1e-15*this->degree*dim;
+ const double eps = 1e-15*q_degree*dim;
#ifdef DEBUG
// in DEBUG mode, check that the evaluation of support points in the
// the tensor product structure of this element and only evaluate 1D
// information from the polynomial. This makes the cost of this function
// almost negligible also for high order elements
- const unsigned int dofs1d = this->degree+1;
+ const unsigned int dofs1d = q_degree+1;
std::vector<Table<2,double> >
subcell_evaluations (dim, Table<2,double>(dofs1d, dofs1d));
const std::vector<unsigned int> &index_map_inverse =
FullMatrix<double> restriction(this->dofs_per_cell, this->dofs_per_cell);
// distinguish q/q_dg0 case
- const unsigned int q_dofs_per_cell = Utilities::fixed_power<dim>(this->degree+1);
+ const unsigned int q_dofs_per_cell = Utilities::fixed_power<dim>(q_degree+1);
// for these Lagrange interpolation polynomials, construction of the
// restriction matrices is relatively simple. the reason is that the
// (compute on one child) by the same value (compute on a later child),
// so we don't have to care about this
- const double eps = 1e-15*this->degree*dim;
+ const double eps = 1e-15*q_degree*dim;
const std::vector<unsigned int> &index_map_inverse =
this->poly_space.get_numbering_inverse();
- const unsigned int dofs1d = this->degree+1;
+ const unsigned int dofs1d = q_degree+1;
std::vector<Tensor<1,dim> > evaluations1d (dofs1d);
restriction.reinit(this->dofs_per_cell, this->dofs_per_cell);
FE_Q_Base<POLY,dim,spacedim>::get_constant_modes () const
{
Table<2,bool> constant_modes(1, this->dofs_per_cell);
- // FE_Q_DG0 should not use this function...
- AssertDimension(this->dofs_per_cell, Utilities::fixed_power<dim>(this->degree+1));
- constant_modes.fill(true);
+ // We here just care for the constant mode due to the polynomial space
+ // without any enrichments
+ // As there may be more constant modes derived classes may to implement this
+ // themselves. An example for this is FE_Q_DG0.
+ for (unsigned int i=0; i<Utilities::fixed_power<dim>(q_degree+1); ++i)
+ constant_modes(0, i) = true;
return std::pair<Table<2,bool>, std::vector<unsigned int> >
(constant_modes, std::vector<unsigned int>(1, 0));
}
#if deal_II_dimension <= deal_II_space_dimension
template class FE_Q_Base<TensorProductPolynomials<deal_II_dimension>, deal_II_dimension, deal_II_space_dimension>;
template class FE_Q_Base<TensorProductPolynomialsConst<deal_II_dimension>, deal_II_dimension, deal_II_space_dimension>;
+ template class FE_Q_Base<TensorProductPolynomialsBubbles<deal_II_dimension>, deal_II_dimension, deal_II_space_dimension>;
template class FE_Q_Base<TensorProductPolynomials<deal_II_dimension,Polynomials::PiecewisePolynomial<double> >, deal_II_dimension, deal_II_space_dimension>;
#endif
}
--- /dev/null
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2012 - 2013 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+#include <deal.II/base/quadrature.h>
+#include <deal.II/base/qprojector.h>
+#include <deal.II/base/template_constraints.h>
+#include <deal.II/fe/fe_q_bubbles.h>
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_tools.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_q1.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+
+#include <vector>
+#include <sstream>
+
+DEAL_II_NAMESPACE_OPEN
+
+
+namespace FE_Q_Bubbles_Helper
+{
+ namespace
+ {
+ template <int dim, int spacedim>
+ inline
+ void
+ compute_embedding_matrices(const FE_Q_Bubbles<dim, spacedim> &fe,
+ std::vector<std::vector<FullMatrix<double> > > &matrices,
+ const bool isotropic_only)
+ {
+ const unsigned int dpc = fe.dofs_per_cell;
+ const unsigned int degree = fe.degree;
+
+ // Initialize quadrature formula on fine cells
+ MappingQ1<dim, spacedim> mapping;
+
+ Quadrature<dim> *q_fine;
+ Quadrature<1> q_dummy(std::vector<Point<1> >(1), std::vector<double> (1,1.));
+ switch (dim)
+ {
+ case 1:
+ if (spacedim==1)
+ q_fine = new QGauss<dim> (degree+1);
+ else if (spacedim==2)
+ q_fine = new QAnisotropic<dim>(QGauss<1>(degree+1), q_dummy);
+ else
+ q_fine = new QAnisotropic<dim>(QGauss<1>(degree+1), q_dummy, q_dummy);
+ break;
+ case 2:
+ if (spacedim==2)
+ q_fine = new QGauss<dim> (degree+1);
+ else
+ q_fine = new QAnisotropic<dim>(QGauss<1>(degree+1), QGauss<1>(degree+1), q_dummy);
+ break;
+ case 3:
+ q_fine = new QGauss<dim> (degree+1);
+ break;
+ default:
+ AssertThrow(false, ExcInternalError());
+ }
+
+ const unsigned int nq = q_fine->size();
+
+ // loop over all possible refinement cases
+ unsigned int ref_case = (isotropic_only)
+ ? RefinementCase<dim>::isotropic_refinement
+ : RefinementCase<dim>::cut_x;
+ for (; ref_case <= RefinementCase<dim>::isotropic_refinement; ++ref_case)
+ {
+ const unsigned int nc
+ = GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case));
+
+ for (unsigned int i=0; i<nc; ++i)
+ {
+ Assert(matrices[ref_case-1][i].n() == dpc,
+ ExcDimensionMismatch(matrices[ref_case-1][i].n(),dpc));
+ Assert(matrices[ref_case-1][i].m() == dpc,
+ ExcDimensionMismatch(matrices[ref_case-1][i].m(),dpc));
+ }
+
+ // create a respective refinement on the triangulation
+ Triangulation<dim, spacedim> tr;
+ GridGenerator::hyper_cube (tr, 0, 1);
+ tr.begin_active()->set_refine_flag(RefinementCase<dim>(ref_case));
+ tr.execute_coarsening_and_refinement();
+
+ DoFHandler<dim, spacedim> dh(tr);
+ dh.distribute_dofs(fe);
+
+ FEValues<dim, spacedim> fine (mapping, fe, *q_fine,
+ update_quadrature_points
+ | update_JxW_values | update_values);
+
+ const unsigned int n_dofs = dh.n_dofs();
+
+ FullMatrix<double> fine_mass(n_dofs);
+ FullMatrix<double> coarse_rhs_matrix(n_dofs, dpc);
+
+ std::vector<std::vector<types::global_dof_index> > child_ldi
+ (nc, std::vector<types::global_dof_index>(fe.dofs_per_cell));
+
+ //now create the mass matrix and all the right_hand sides
+ unsigned int child_no = 0;
+ typename dealii::DoFHandler<dim>::active_cell_iterator cell
+ = dh.begin_active();
+ for (; cell!=dh.end(); ++cell, ++child_no)
+ {
+ fine.reinit(cell);
+ cell->get_dof_indices(child_ldi[child_no]);
+
+ for (unsigned int q=0; q<nq; ++q)
+ for (unsigned int i=0; i<dpc; ++i)
+ for (unsigned int j=0; j<dpc; ++j)
+ {
+ const unsigned int gdi=child_ldi[child_no][i];
+ const unsigned int gdj=child_ldi[child_no][j];
+ fine_mass(gdi, gdj)+=fine.shape_value(i,q)
+ *fine.shape_value(j,q)
+ *fine.JxW(q);
+ Point<dim> quad_tmp;
+ for (unsigned int k=0; k<dim; ++k)
+ quad_tmp(k) = fine.quadrature_point(q)(k);
+ coarse_rhs_matrix(gdi, j)
+ +=fine.shape_value(i,q)
+ *fe.shape_value(j, quad_tmp)
+ *fine.JxW(q);
+ }
+ }
+
+ //now solve for all right-hand sides simultaneously
+ dealii::FullMatrix<double> solution (n_dofs, dpc);
+ fine_mass.gauss_jordan();
+ fine_mass.mmult(solution, coarse_rhs_matrix);
+
+ //and distribute to the fine cell matrices
+ for (unsigned int child_no=0; child_no<nc; ++child_no)
+ for (unsigned int i=0; i<dpc; ++i)
+ for (unsigned int j=0; j<dpc; ++j)
+ {
+ const unsigned int gdi=child_ldi[child_no][i];
+ //remove small entries
+ if (std::fabs(solution(gdi, j)) > 1.e-12)
+ matrices[ref_case-1][child_no](i,j)=solution(gdi, j);
+ }
+ }
+ }
+ }
+}
+
+
+template <int dim, int spacedim>
+FE_Q_Bubbles<dim,spacedim>::FE_Q_Bubbles (const unsigned int q_degree)
+ :
+ FE_Q_Base<TensorProductPolynomialsBubbles<dim>, dim, spacedim> (
+ TensorProductPolynomialsBubbles<dim>(Polynomials::LagrangeEquidistant::generate_complete_basis(q_degree)),
+ FiniteElementData<dim>(get_dpo_vector(q_degree),
+ 1, q_degree+1,
+ FiniteElementData<dim>::H1),
+ get_riaf_vector(q_degree)),
+ n_bubbles((q_degree<=1)?1:dim)
+{
+ Assert (q_degree > 0,
+ ExcMessage ("This element can only be used for polynomial degrees "
+ "greater than zero"));
+
+ std::vector<Point<1> > support_points_1d(q_degree+1);
+ for (unsigned int i=0; i<=q_degree; ++i)
+ support_points_1d[i][0] = static_cast<double>(i)/q_degree;
+
+ this->initialize(support_points_1d);
+
+ // adjust unit support point for discontinuous node
+ Point<dim> point;
+ for (unsigned int d=0; d<dim; ++d)
+ point[d] = 0.5;
+ for (unsigned int i=0; i<n_bubbles; ++i)
+ this->unit_support_points.push_back(point);
+ AssertDimension(this->dofs_per_cell, this->unit_support_points.size());
+
+ this->reinit_restriction_and_prolongation_matrices();
+ if (dim == spacedim)
+ {
+ FE_Q_Bubbles_Helper::compute_embedding_matrices
+ (*this, this->prolongation, false);
+ // Fill restriction matrices with L2-projection
+ FETools::compute_projection_matrices (*this, this->restriction);
+ }
+
+}
+
+
+
+template <int dim, int spacedim>
+FE_Q_Bubbles<dim,spacedim>::FE_Q_Bubbles (const Quadrature<1> &points)
+ :
+ FE_Q_Base<TensorProductPolynomialsBubbles<dim>, dim, spacedim> (
+ TensorProductPolynomialsBubbles<dim>(Polynomials::generate_complete_Lagrange_basis(points.get_points())),
+ FiniteElementData<dim>(get_dpo_vector(points.size()-1),
+ 1, points.size(),
+ FiniteElementData<dim>::H1),
+ get_riaf_vector(points.size()-1)),
+ n_bubbles((points.size()-1<=1)?1:dim)
+{
+ const unsigned int q_degree = points.size()-1;
+ (void) q_degree;
+ Assert (q_degree > 0,
+ ExcMessage ("This element can only be used for polynomial degrees "
+ "at least one"));
+
+ this->initialize(points.get_points());
+
+ // adjust unit support point for discontinuous node
+ Point<dim> point;
+ for (unsigned int d=0; d<dim; ++d)
+ point[d] = 0.5;
+ for (unsigned int i=0; i< n_bubbles; ++i)
+ this->unit_support_points.push_back(point);
+ AssertDimension(this->dofs_per_cell, this->unit_support_points.size());
+
+ this->reinit_restriction_and_prolongation_matrices();
+ if (dim == spacedim)
+ {
+ FE_Q_Bubbles_Helper::compute_embedding_matrices
+ (*this, this->prolongation, false);
+ // Fill restriction matrices with L2-projection
+ FETools::compute_projection_matrices (*this, this->restriction);
+ }
+}
+
+
+
+template <int dim, int spacedim>
+std::string
+FE_Q_Bubbles<dim,spacedim>::get_name () const
+{
+ // note that the FETools::get_fe_from_name function depends on the
+ // particular format of the string this function returns, so they have to be
+ // kept in synch
+
+ std::ostringstream namebuf;
+ bool type = true;
+ const unsigned int n_points = this->degree;
+ std::vector<double> points(n_points);
+ const unsigned int dofs_per_cell = this->dofs_per_cell;
+ const std::vector<Point<dim> > &unit_support_points = this->unit_support_points;
+ unsigned int index = 0;
+
+ // Decode the support points in one coordinate direction.
+ for (unsigned int j=0; j<dofs_per_cell; j++)
+ {
+ if ((dim>1) ? (unit_support_points[j](1)==0 &&
+ ((dim>2) ? unit_support_points[j](2)==0: true)) : true)
+ {
+ if (index == 0)
+ points[index] = unit_support_points[j](0);
+ else if (index == 1)
+ points[n_points-1] = unit_support_points[j](0);
+ else
+ points[index-1] = unit_support_points[j](0);
+
+ index++;
+ }
+ }
+ // Do not consider the discontinuous node for dimension 1
+ Assert (index == n_points || (dim==1 && index == n_points+n_bubbles),
+ ExcMessage ("Could not decode support points in one coordinate direction."));
+
+ // Check whether the support points are equidistant.
+ for (unsigned int j=0; j<n_points; j++)
+ if (std::fabs(points[j] - (double)j/(this->degree-1)) > 1e-15)
+ {
+ type = false;
+ break;
+ }
+
+ if (type == true)
+ namebuf << "FE_Q_Bubbles<" << dim << ">(" << this->degree-1 << ")";
+ else
+ {
+
+ // Check whether the support points come from QGaussLobatto.
+ const QGaussLobatto<1> points_gl(n_points);
+ type = true;
+ for (unsigned int j=0; j<n_points; j++)
+ if (points[j] != points_gl.point(j)(0))
+ {
+ type = false;
+ break;
+ }
+ if (type == true)
+ namebuf << "FE_Q_Bubbles<" << dim << ">(QGaussLobatto(" << this->degree << "))";
+ else
+ namebuf << "FE_Q_Bubbles<" << dim << ">(QUnknownNodes(" << this->degree-1 << "))";
+ }
+ return namebuf.str();
+}
+
+
+
+template <int dim, int spacedim>
+FiniteElement<dim,spacedim> *
+FE_Q_Bubbles<dim,spacedim>::clone() const
+{
+ return new FE_Q_Bubbles<dim,spacedim>(*this);
+}
+
+
+
+template <int dim, int spacedim>
+void
+FE_Q_Bubbles<dim,spacedim>::interpolate(std::vector<double> &local_dofs,
+ const std::vector<double> &values) const
+{
+ Assert (values.size() == this->unit_support_points.size(),
+ ExcDimensionMismatch(values.size(),
+ this->unit_support_points.size()));
+ Assert (local_dofs.size() == this->dofs_per_cell,
+ ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell));
+ Assert (this->n_components() == 1,
+ ExcDimensionMismatch(this->n_components(), 1));
+
+ std::copy(values.begin(), values.end(), local_dofs.begin());
+
+ // We don't use the bubble functions for local interpolation
+ for (unsigned int i = 0; i<n_bubbles; ++i)
+ local_dofs[local_dofs.size()-i-1] = 0.;
+}
+
+
+
+template <int dim, int spacedim>
+void
+FE_Q_Bubbles<dim,spacedim>::interpolate(std::vector<double> &local_dofs,
+ const std::vector<Vector<double> > &values,
+ unsigned int offset) const
+{
+ Assert (values.size() == this->unit_support_points.size(),
+ ExcDimensionMismatch(values.size(),
+ this->unit_support_points.size()));
+ Assert (local_dofs.size() == this->dofs_per_cell,
+ ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell));
+ Assert (values[0].size() >= offset+this->n_components(),
+ ExcDimensionMismatch(values[0].size(),offset+this->n_components()));
+
+ for (unsigned int i=0; i<this->dofs_per_cell-1; ++i)
+ {
+ const std::pair<unsigned int, unsigned int> index
+ = this->system_to_component_index(i);
+ local_dofs[i] = values[i](offset+index.first);
+ }
+
+ // We don't use the bubble functions for local interpolation
+ for (unsigned int i = 0; i<n_bubbles; ++i)
+ local_dofs[local_dofs.size()-i-1] = 0.;
+}
+
+
+
+template <int dim, int spacedim>
+void
+FE_Q_Bubbles<dim,spacedim>::interpolate(
+ std::vector<double> &local_dofs,
+ const VectorSlice<const std::vector<std::vector<double> > > &values) const
+{
+ Assert (values[0].size() == this->unit_support_points.size(),
+ ExcDimensionMismatch(values.size(),
+ this->unit_support_points.size()));
+ Assert (local_dofs.size() == this->dofs_per_cell,
+ ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell));
+ Assert (values.size() == this->n_components(),
+ ExcDimensionMismatch(values.size(), this->n_components()));
+
+ for (unsigned int i=0; i<this->dofs_per_cell-1; ++i)
+ {
+ const std::pair<unsigned int, unsigned int> index
+ = this->system_to_component_index(i);
+ local_dofs[i] = values[index.first][i];
+ }
+
+ // We don't use the bubble functions for local interpolation
+ for (unsigned int i = 0; i<n_bubbles; ++i)
+ local_dofs[local_dofs.size()-i-1] = 0.;
+}
+
+
+
+template <int dim, int spacedim>
+void
+FE_Q_Bubbles<dim,spacedim>::
+get_interpolation_matrix (const FiniteElement<dim,spacedim> &x_source_fe,
+ FullMatrix<double> &interpolation_matrix) const
+{
+ // We don't know how to do this properly, yet.
+ // However, for SolutionTransfer to work we need to provide an implementation
+ // for the case that the x_source_fe is identical to this FE
+ typedef FE_Q_Bubbles<dim,spacedim> FEQBUBBLES;
+ typedef FiniteElement<dim,spacedim> FEL;
+
+ AssertThrow ((x_source_fe.get_name().find ("FE_Q_Bubbles<") == 0)
+ ||
+ (dynamic_cast<const FEQBUBBLES *>(&x_source_fe) != 0)
+ ,
+ typename FEL::
+ ExcInterpolationNotImplemented());
+ Assert (interpolation_matrix.m() == this->dofs_per_cell,
+ ExcDimensionMismatch (interpolation_matrix.m(),
+ this->dofs_per_cell));
+ Assert (interpolation_matrix.n() == x_source_fe.dofs_per_cell,
+ ExcDimensionMismatch (interpolation_matrix.m(),
+ x_source_fe.dofs_per_cell));
+
+ //Provide a short cut in case we are just inquiring the identity
+ if (dynamic_cast<const FEQBUBBLES *>(&x_source_fe)->degree == this->degree)
+ for (unsigned int i=0; i<interpolation_matrix.m(); ++i)
+ interpolation_matrix.set(i,i,1.);
+ //else we need to do more...
+ else
+ Assert(false, typename FEL::ExcInterpolationNotImplemented())
+ }
+
+
+
+template <int dim, int spacedim>
+std::vector<bool>
+FE_Q_Bubbles<dim,spacedim>::get_riaf_vector(const unsigned int q_deg)
+{
+ unsigned int n_cont_dofs = Utilities::fixed_power<dim>(q_deg+1);
+ const unsigned int n_bubbles = (q_deg<=1?1:dim);
+ return std::vector<bool> (n_cont_dofs+n_bubbles,true);
+}
+
+
+
+template <int dim, int spacedim>
+std::vector<unsigned int>
+FE_Q_Bubbles<dim,spacedim>::get_dpo_vector(const unsigned int q_deg)
+{
+ std::vector<unsigned int> dpo(dim+1, 1U);
+ for (unsigned int i=1; i<dpo.size(); ++i)
+ dpo[i]=dpo[i-1]*(q_deg-1);
+
+ dpo[dim]+=(q_deg<=1?1:dim);//all the bubble functions are discontinuous
+ return dpo;
+}
+
+
+
+template <int dim, int spacedim>
+bool
+FE_Q_Bubbles<dim,spacedim>::has_support_on_face
+(const unsigned int shape_index,
+ const unsigned int face_index) const
+{
+ // discontinuous functions have no support on faces
+ if (shape_index >= this->n_dofs_per_cell()-n_bubbles)
+ return false;
+ else
+ return FE_Q_Base<TensorProductPolynomialsBubbles<dim>,dim,spacedim>::has_support_on_face(shape_index, face_index);
+}
+
+
+
+template <int dim, int spacedim>
+const FullMatrix<double> &
+FE_Q_Bubbles<dim,spacedim>::get_prolongation_matrix
+(const unsigned int child,
+ const RefinementCase<dim> &refinement_case) const
+{
+ Assert (refinement_case<RefinementCase<dim>::isotropic_refinement+1,
+ ExcIndexRange(refinement_case,0,RefinementCase<dim>::isotropic_refinement+1));
+ Assert (refinement_case!=RefinementCase<dim>::no_refinement,
+ ExcMessage("Prolongation matrices are only available for refined cells!"));
+ Assert (child<GeometryInfo<dim>::n_children(refinement_case),
+ ExcIndexRange(child,0,GeometryInfo<dim>::n_children(refinement_case)));
+
+ Assert (this->prolongation[refinement_case-1][child].n() != 0,
+ ExcMessage("This prolongation matrix has not been computed yet!"));
+ // finally return the matrix
+ return this->prolongation[refinement_case-1][child];
+}
+
+
+
+template <int dim, int spacedim>
+const FullMatrix<double> &
+FE_Q_Bubbles<dim,spacedim>::get_restriction_matrix
+(const unsigned int child,
+ const RefinementCase<dim> &refinement_case) const
+{
+ Assert (refinement_case<RefinementCase<dim>::isotropic_refinement+1,
+ ExcIndexRange(refinement_case,0,RefinementCase<dim>::isotropic_refinement+1));
+ Assert (refinement_case!=RefinementCase<dim>::no_refinement,
+ ExcMessage("Restriction matrices are only available for refined cells!"));
+ Assert (child<GeometryInfo<dim>::n_children(refinement_case),
+ ExcIndexRange(child,0,GeometryInfo<dim>::n_children(refinement_case)));
+
+ Assert(this->restriction[refinement_case-1][child].n() != 0,
+ ExcMessage("This restriction matrix has not been computed yet!"));
+
+ //finally return the matrix
+ return this->restriction[refinement_case-1][child];
+}
+
+// explicit instantiations
+#include "fe_q_bubbles.inst"
+
+DEAL_II_NAMESPACE_CLOSE
--- /dev/null
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 1998 - 2013 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
+ {
+#if deal_II_dimension <= deal_II_space_dimension
+ template class FE_Q_Bubbles<deal_II_dimension, deal_II_space_dimension>;
+#endif
+ }
projects / dealii.git / commitdiff
