]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implement FE_Q_Bubbles that is based on TensorProductPolynomialsBubbles
authorDaniel Arndt <d.arndt@math.uni-goettingen.de>
Wed, 12 Aug 2015 09:50:50 +0000 (04:50 -0500)
committerDaniel Arndt <d.arndt@math.uni-goettingen.de>
Wed, 12 Aug 2015 11:38:57 +0000 (06:38 -0500)
12 files changed:
doc/doxygen/images/fe_q_bubbles_conditioning.png [new file with mode: 0644]
doc/news/changes.h
include/deal.II/fe/fe_poly.templates.h
include/deal.II/fe/fe_q_base.h
include/deal.II/fe/fe_q_bubbles.h [new file with mode: 0644]
source/fe/CMakeLists.txt
source/fe/fe_poly.cc
source/fe/fe_poly.inst.in
source/fe/fe_q_base.cc
source/fe/fe_q_base.inst.in
source/fe/fe_q_bubbles.cc [new file with mode: 0644]
source/fe/fe_q_bubbles.inst.in [new file with mode: 0644]

diff --git a/doc/doxygen/images/fe_q_bubbles_conditioning.png b/doc/doxygen/images/fe_q_bubbles_conditioning.png
new file mode 100644 (file)
index 0000000..a061872
Binary files /dev/null and b/doc/doxygen/images/fe_q_bubbles_conditioning.png differ
index 2db0318d828f0ae38c0869720c58d888d31b83eb..225378735879bf8cce1a622ce1162a71de80d380 100644 (file)
@@ -66,6 +66,11 @@ inconvenience this causes.
 
 
 <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
index d93c690e5867f1986b4cf0609a87cc40ee747c74..ce2a7bf4d66849f6cace81ab9f673d804d7f7b81 100644 (file)
@@ -18,6 +18,7 @@
 #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>
 
index b89de91721f6da911a4d86becd18500d914554a7..a5a9f8337709ebe8760ec329a7770cf047605d55 100644 (file)
@@ -27,8 +27,8 @@ DEAL_II_NAMESPACE_OPEN
 /*@{*/
 
 /**
- * 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,
@@ -324,6 +324,13 @@ private:
    * 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;
 };
 
 
diff --git a/include/deal.II/fe/fe_q_bubbles.h b/include/deal.II/fe/fe_q_bubbles.h
new file mode 100644 (file)
index 0000000..1e6b9cd
--- /dev/null
@@ -0,0 +1,205 @@
+// ---------------------------------------------------------------------
+// $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
index b4130f5407ce99b15086bedfeff1d8feb2b2aeb9..73fee37d1e3a7b0065367cafbb184bdeae5e23c0 100644 (file)
@@ -35,6 +35,7 @@ SET(_src
   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
@@ -72,6 +73,7 @@ SET(_inst
   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
index 147b4eb984aac45cfcb9c245d9adb65d8332d635..ebd25e79063c0d4411c2710ae769b69402d73632 100644 (file)
@@ -17,6 +17,7 @@
 #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>
index 6331952dd9bc3cf849eaf5832ea6c287445dc536..dcd33ed098ef6011c91da0a64ee403a57914e940 100644 (file)
@@ -1,6 +1,6 @@
 // ---------------------------------------------------------------------
 //
-// 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.
 //
@@ -20,6 +20,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS
 #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>;
index 866290fac412b646bb17db5568ec38ea51ba7f74..f6eb9fdbdc3042695cc3284ac57e24e88d91ee54 100644 (file)
@@ -19,6 +19,7 @@
 #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>
@@ -124,6 +125,10 @@ struct FE_Q_Base<POLY,xdim,xspacedim>::Implementation
   {
     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.
@@ -172,10 +177,10 @@ struct FE_Q_Base<POLY,xdim,xspacedim>::Implementation
     // 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 (
@@ -197,13 +202,13 @@ struct FE_Q_Base<POLY,xdim,xspacedim>::Implementation
     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);
@@ -227,6 +232,10 @@ struct FE_Q_Base<POLY,xdim,xspacedim>::Implementation
   {
     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.
 
@@ -253,10 +262,10 @@ struct FE_Q_Base<POLY,xdim,xspacedim>::Implementation
     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;
@@ -310,14 +319,14 @@ struct FE_Q_Base<POLY,xdim,xspacedim>::Implementation
 
     // 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());
@@ -327,7 +336,7 @@ struct FE_Q_Base<POLY,xdim,xspacedim>::Implementation
 
     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;
 
@@ -373,14 +382,14 @@ struct FE_Q_Base<POLY,xdim,xspacedim>::Implementation
 
         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],
@@ -405,7 +414,10 @@ FE_Q_Base<POLY,dim,spacedim>::FE_Q_Base (const POLY &poly_space,
                                          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)
 {}
 
 
@@ -422,17 +434,18 @@ FE_Q_Base<POLY,dim,spacedim>::initialize (const std::vector<Point<1> > &points)
   // 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);
   }
 
@@ -471,7 +484,7 @@ get_interpolation_matrix (const FiniteElement<dim,spacedim> &x_source_fe,
                                     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
@@ -504,7 +517,7 @@ get_interpolation_matrix (const FiniteElement<dim,spacedim> &x_source_fe,
         }
 
       // 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)
@@ -596,7 +609,7 @@ get_subface_interpolation_matrix (const FiniteElement<dim,spacedim> &x_source_fe
       // 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.
@@ -888,11 +901,11 @@ void FE_Q_Base<POLY,dim,spacedim>::initialize_unit_face_support_points (const st
     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());
@@ -914,7 +927,7 @@ FE_Q_Base<POLY,dim,spacedim>::initialize_quad_dof_index_permutation ()
   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
@@ -1137,7 +1150,7 @@ FE_Q_Base<POLY,dim,spacedim>
         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
@@ -1145,7 +1158,7 @@ FE_Q_Base<POLY,dim,spacedim>
       // 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
@@ -1167,7 +1180,7 @@ FE_Q_Base<POLY,dim,spacedim>
       // 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 =
@@ -1315,7 +1328,7 @@ FE_Q_Base<POLY,dim,spacedim>
 
       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
@@ -1339,11 +1352,11 @@ FE_Q_Base<POLY,dim,spacedim>
       // (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);
@@ -1531,9 +1544,12 @@ std::pair<Table<2,bool>, std::vector<unsigned int> >
 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));
 }
index f0b13ce5c1bac62c8064a75cc51edec9a6f04f09..1e36e6bee3ad64ca19feb76c372d3992fa2b5bcc 100644 (file)
@@ -20,6 +20,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS
 #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
   }
diff --git a/source/fe/fe_q_bubbles.cc b/source/fe/fe_q_bubbles.cc
new file mode 100644 (file)
index 0000000..7616ac5
--- /dev/null
@@ -0,0 +1,527 @@
+// ---------------------------------------------------------------------
+// $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
diff --git a/source/fe/fe_q_bubbles.inst.in b/source/fe/fe_q_bubbles.inst.in
new file mode 100644 (file)
index 0000000..7cf26de
--- /dev/null
@@ -0,0 +1,24 @@
+// ---------------------------------------------------------------------
+// $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
+  }

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.