Split up error_estimator into two files, improve the generic face estimate part so...

git-svn-id: https://svn.dealii.org/trunk@30400 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/include/deal.II/dofs/dof_accessor.templates.h b/deal.II/include/deal.II/dofs/dof_accessor.templates.h
index 8006dba0d046b583338f78646f66c8a8bd3fa802..989591f843a11f8edac6caf2b03c71df820c4b9e 100644 (file)
--- a/deal.II/include/deal.II/dofs/dof_accessor.templates.h
+++ b/deal.II/include/deal.II/dofs/dof_accessor.templates.h
@@ -3492,6 +3492,7 @@ namespace internal
 
 
 template <class DH, bool lda>
+inline
 TriaIterator<DoFAccessor<DH::dimension-1,DH,lda> >
 DoFCellAccessor<DH,lda>::face (const unsigned int i) const
 {
@@ -3519,6 +3520,7 @@ DoFCellAccessor<DH,lda>::get_dof_indices (std::vector<types::global_dof_index> &
 }
 
 
+
 template<class DH, bool lda>
 inline
 void DoFCellAccessor<DH,lda>::get_mg_dof_indices (std::vector<types::global_dof_index> &dof_indices) const
@@ -3526,6 +3528,8 @@ void DoFCellAccessor<DH,lda>::get_mg_dof_indices (std::vector<types::global_dof_
   DoFAccessor<dim, DH,lda>::get_mg_dof_indices (this->level (), dof_indices);
 }
 
+
+
 template<class DH, bool lda>
 inline
 void DoFCellAccessor<DH,lda>::set_mg_dof_indices (const std::vector<types::global_dof_index> &dof_indices)
@@ -3534,6 +3538,7 @@ void DoFCellAccessor<DH,lda>::set_mg_dof_indices (const std::vector<types::globa
 }
 
 
+
 template<class DH, bool lda>
 inline
 void DoFCellAccessor<DH,lda>::get_active_or_mg_dof_indices (std::vector<types::global_dof_index> &dof_indices) const
@@ -3548,6 +3553,7 @@ void DoFCellAccessor<DH,lda>::get_active_or_mg_dof_indices (std::vector<types::g
 
 template <class DH, bool lda>
 template <class InputVector, typename number>
+inline
 void
 DoFCellAccessor<DH,lda>::get_dof_values (const InputVector &values,
                                          Vector<number>    &local_values) const
@@ -3564,6 +3570,7 @@ DoFCellAccessor<DH,lda>::get_dof_values (const InputVector &values,
 
 template <class DH, bool lda>
 template <class InputVector, typename ForwardIterator>
+inline
 void
 DoFCellAccessor<DH,lda>::get_dof_values (const InputVector &values,
                                          ForwardIterator    local_values_begin,
@@ -3579,6 +3586,7 @@ DoFCellAccessor<DH,lda>::get_dof_values (const InputVector &values,
 
 template <class DH, bool lda>
 template <class InputVector, typename ForwardIterator>
+inline
 void
 DoFCellAccessor<DH,lda>::get_dof_values (const ConstraintMatrix &constraints,
                                          const InputVector      &values,
@@ -3596,6 +3604,7 @@ DoFCellAccessor<DH,lda>::get_dof_values (const ConstraintMatrix &constraints,
 
 template <class DH, bool lda>
 template <class OutputVector, typename number>
+inline
 void
 DoFCellAccessor<DH,lda>::set_dof_values (const Vector<number> &local_values,
                                          OutputVector         &values) const

diff --git a/deal.II/include/deal.II/numerics/error_estimator.h b/deal.II/include/deal.II/numerics/error_estimator.h
index f47c98b30b0d3ce07c08c40ecfcd1fe49c3a9d88..6f15d97c2e759c8af141e3982cc2cfd6f5c9b5ac 100644 (file)
--- a/deal.II/include/deal.II/numerics/error_estimator.h
+++ b/deal.II/include/deal.II/numerics/error_estimator.h
@@ -234,114 +234,59 @@ class KellyErrorEstimator
 {
 public:
   /**
-   * Implementation of the error
-   * estimator described above. You
-   * may give a coefficient, but
-   * there is a default value which
-   * denotes the constant
-   * coefficient with value
-   * one. The coefficient function
-   * may either be a scalar one, in
-   * which case it is used for all
-   * components of the finite
-   * element, or a vector-valued
-   * one with as many components as
-   * there are in the finite
-   * element; in the latter case,
-   * each component is weighted by
-   * the respective component in
-   * the coefficient.
+   * Implementation of the error estimator described above. You may give a
+   * coefficient, but there is a default value which denotes the constant
+   * coefficient with value one. The coefficient function may either be a
+   * scalar one, in which case it is used for all components of the finite
+   * element, or a vector-valued one with as many components as there are in
+   * the finite element; in the latter case, each component is weighted by the
+   * respective component in the coefficient.
    *
-   * You might give a list of
-   * components you want to
-   * evaluate, in case the finite
-   * element used by the
-   * DoFHandler object is
-   * vector-valued. You then have
-   * to set those entries to true
-   * in the bit-vector
-   * @p component_mask
-   * (see @ref GlossComponentMask)
-   * for which the
-   * respective component is to be
-   * used in the error
-   * estimator. The default is to
-   * use all components, which is
-   * done by either providing a
-   * bit-vector with all-set
-   * entries, or an empty
-   * bit-vector.
+   * You might give a list of components you want to evaluate, in case the
+   * finite element used by the DoFHandler object is vector-valued. You then
+   * have to set those entries to true in the bit-vector @p component_mask
+   * (see @ref GlossComponentMask) for which the respective component is to be
+   * used in the error estimator. The default is to use all components, which
+   * is done by either providing a bit-vector with all-set entries, or an
+   * empty bit-vector.
    *
-   * The @p subdomain_id parameter
-   * indicates whether we shall compute
-   * indicators for all cells (in case its
-   * value is the default,
-   * <tt>numbers::invalid_unsigned_int</tt>),
-   * or only for the cells belonging to a
-   * certain subdomain with the given
-   * indicator. The latter case is used for
-   * parallel computations where all
-   * processor nodes have the global grid
-   * stored, and could well compute all the
-   * indicators for all cells themselves,
-   * but where it is more efficient to have
-   * each process compute only indicators
-   * for the cells it owns, and have them
-   * exchange the resulting information
-   * afterwards. This is in particular true
-   * for the case where meshes are very
-   * large and computing indicators for @em
-   * every cell is too expensive, while
-   * computing indicators for only local
-   * cells is acceptable. Note that if you
-   * only ask for the indicators of a
-   * certain subdomain to be computed, you
-   * must nevertheless make sure that this
-   * function has access to the correct
-   * node values of @em all degrees of
-   * freedom. This is since the function
-   * needs to access DoF values on
-   * neighboring cells as well, even if
-   * they belong to a different subdomain.
+   * The @p subdomain_id parameter indicates whether we shall compute
+   * indicators for all cells (in case its value is the default,
+   * <tt>numbers::invalid_unsigned_int</tt>), or only for the cells belonging
+   * to a certain subdomain with the given indicator. The latter case is used
+   * for parallel computations where all processor nodes have the global grid
+   * stored, and could well compute all the indicators for all cells
+   * themselves, but where it is more efficient to have each process compute
+   * only indicators for the cells it owns, and have them exchange the
+   * resulting information afterwards. This is in particular true for the case
+   * where meshes are very large and computing indicators for @em every cell
+   * is too expensive, while computing indicators for only local cells is
+   * acceptable. Note that if you only ask for the indicators of a certain
+   * subdomain to be computed, you must nevertheless make sure that this
+   * function has access to the correct node values of @em all degrees of
+   * freedom. This is since the function needs to access DoF values on
+   * neighboring cells as well, even if they belong to a different subdomain.
    *
-   * The @p material_id parameter has a similar
-   * meaning: if not set to its default value
-   * (which is numbers::invalid_material_id),
-   * it means that indicators will only be
-   * computed for cells with this particular
-   * material id.
+   * The @p material_id parameter has a similar meaning: if not set to its
+   * default value (which is numbers::invalid_material_id), it means that
+   * indicators will only be computed for cells with this particular material
+   * id.
    *
-   * The @p n_threads parameter used to
-   * indicate the number of threads to be
-   * used to compute the error
-   * estimator. This parameter is now
-   * ignored, with the number of threads
-   * determined automatically. The
-   * parameter is retained for
-   * compatibility with old versions of the
-   * library.
+   * The @p n_threads parameter used to indicate the number of threads to be
+   * used to compute the error estimator. This parameter is now ignored, with
+   * the number of threads determined automatically. The parameter is retained
+   * for compatibility with old versions of the library.
    *
-   * @note If the DoFHandler object
-   * given as an argument to this
-   * function builds on a
-   * parallel::distributed::Triangulation,
-   * this function skips
-   * computations on all cells that
-   * are not locally owned. In that
-   * case, the only valid value for
-   * the subdomain_id argument
-   * (besides the invalid value) is
-   * the subdomain id that is
-   * associated with the currently
+   * @note If the DoFHandler object given as an argument to this function
+   * builds on a parallel::distributed::Triangulation, this function skips
+   * computations on all cells that are not locally owned. In that case, the
+   * only valid value for the subdomain_id argument (besides the invalid
+   * value) is the subdomain id that is associated with the currently
    * processor, as reported by
    * parallel::distributed::Triangulation::locally_owned_subdomain(). Even
-   * though nothing is computed on
-   * cells that we don't locally
-   * own, the error indicator
-   * vector must still have a
-   * length equal to the number of
-   * active cell in the mesh as
-   * reported by
+   * though nothing is computed on cells that we don't locally own, the error
+   * indicator vector must still have a length equal to the number of active
+   * cell in the mesh as reported by
    * parallel::distributed::Triangulation::n_locally_owned_active_cells().
    */
   template <typename InputVector, class DH>
@@ -358,8 +303,7 @@ public:
                         const types::material_id       material_id = numbers::invalid_material_id);
 
   /**
-   * Calls the @p estimate
-   * function, see above, with
+   * Calls the @p estimate function, see above, with
    * <tt>mapping=MappingQ1@<dim@>()</tt>.
    */
   template <typename InputVector, class DH>
@@ -375,31 +319,17 @@ public:
                         const types::material_id       material_id = numbers::invalid_material_id);
 
   /**
-   * Same function as above, but
-   * accepts more than one solution
-   * vector and returns one error
-   * vector for each solution
-   * vector. For the reason of
-   * existence of this function,
-   * see the general documentation
-   * of this class.
+   * Same function as above, but accepts more than one solution vector and
+   * returns one error vector for each solution vector. For the reason of
+   * existence of this function, see the general documentation of this class.
    *
-   * Since we do not want to force
-   * the user of this function to
-   * copy around their solution
-   * vectors, the vector of
-   * solution vectors takes
-   * pointers to the solutions,
-   * rather than being a vector of
-   * vectors. This makes it simpler
-   * to have the solution vectors
-   * somewhere in memory, rather
-   * than to have them collected
-   * somewhere special. (Note that
-   * it is not possible to
-   * construct of vector of
-   * references, so we had to use a
-   * vector of pointers.)
+   * Since we do not want to force the user of this function to copy around
+   * their solution vectors, the vector of solution vectors takes pointers to
+   * the solutions, rather than being a vector of vectors. This makes it
+   * simpler to have the solution vectors somewhere in memory, rather than to
+   * have them collected somewhere special. (Note that it is not possible to
+   * construct of vector of references, so we had to use a vector of
+   * pointers.)
    */
   template <typename InputVector, class DH>
   static void estimate (const Mapping<dim, spacedim>          &mapping,
@@ -415,8 +345,7 @@ public:
                         const types::material_id           material_id = numbers::invalid_material_id);
 
   /**
-   * Calls the @p estimate
-   * function, see above, with
+   * Calls the @p estimate function, see above, with
    * <tt>mapping=MappingQ1@<dim@>()</tt>.
    */
   template <typename InputVector, class DH>
@@ -433,10 +362,8 @@ public:
 
 
   /**
-   * Equivalent to the set of functions
-   * above, except that this one takes a
-   * quadrature collection for hp finite
-   * element dof handlers.
+   * Equivalent to the set of functions above, except that this one takes a
+   * quadrature collection for hp finite element dof handlers.
    */
   template <typename InputVector, class DH>
   static void estimate (const Mapping<dim, spacedim>      &mapping,
@@ -453,10 +380,8 @@ public:
 
 
   /**
-   * Equivalent to the set of functions
-   * above, except that this one takes a
-   * quadrature collection for hp finite
-   * element dof handlers.
+   * Equivalent to the set of functions above, except that this one takes a
+   * quadrature collection for hp finite element dof handlers.
    */
   template <typename InputVector, class DH>
   static void estimate (const DH                &dof,
@@ -472,10 +397,8 @@ public:
 
 
   /**
-   * Equivalent to the set of functions
-   * above, except that this one takes a
-   * quadrature collection for hp finite
-   * element dof handlers.
+   * Equivalent to the set of functions above, except that this one takes a
+   * quadrature collection for hp finite element dof handlers.
    */
   template <typename InputVector, class DH>
   static void estimate (const Mapping<dim, spacedim>          &mapping,
@@ -492,10 +415,8 @@ public:
 
 
   /**
-   * Equivalent to the set of functions
-   * above, except that this one takes a
-   * quadrature collection for hp finite
-   * element dof handlers.
+   * Equivalent to the set of functions above, except that this one takes a
+   * quadrature collection for hp finite element dof handlers.
    */
   template <typename InputVector, class DH>
   static void estimate (const DH                    &dof,
@@ -561,51 +482,28 @@ class KellyErrorEstimator<1,spacedim>
 {
 public:
   /**
-   * Implementation of the error
-   * estimator described above. You
-   * may give a coefficient, but
-   * there is a default value which
-   * denotes the constant
-   * coefficient with value
-   * one. The coefficient function
-   * may either be a scalar one, in
-   * which case it is used for all
-   * components of the finite
-   * element, or a vector-valued
-   * one with as many components as
-   * there are in the finite
-   * element; in the latter case,
-   * each component is weighted by
-   * the respective component in
-   * the coefficient.
+   * Implementation of the error estimator described above. You may give a
+   * coefficient, but there is a default value which denotes the constant
+   * coefficient with value one. The coefficient function may either be a
+   * scalar one, in which case it is used for all components of the finite
+   * element, or a vector-valued one with as many components as there are in
+   * the finite element; in the latter case, each component is weighted by the
+   * respective component in the coefficient.
    *
-   * You might give a list of components
-   * you want to evaluate, in case the
-   * finite element used by the DoFHandler
-   * object is vector-valued. You then have
-   * to set those entries to true in the
-   * bit-vector @p component_mask for which
-   * the respective component is to be used
-   * in the error estimator. The default is
-   * to use all components, which is done
-   * by either providing a bit-vector with
-   * all-set entries, or an empty
-   * bit-vector. All the other parameters
-   * are as in the general case used for 2d
-   * and higher.
+   * You might give a list of components you want to evaluate, in case the
+   * finite element used by the DoFHandler object is vector-valued. You then
+   * have to set those entries to true in the bit-vector @p component_mask for
+   * which the respective component is to be used in the error estimator. The
+   * default is to use all components, which is done by either providing a
+   * bit-vector with all-set entries, or an empty bit-vector. All the other
+   * parameters are as in the general case used for 2d and higher.
    *
-   * The estimator supports multithreading
-   * and splits the cells to
-   * <tt>multithread_info.n_default_threads</tt>
-   * (default) threads. The number of
-   * threads to be used in multithreaded
-   * mode can be set with the last
-   * parameter of the error estimator.
-   * Multithreading is not presently
-   * implemented for 1d, but we retain the
-   * respective parameter for compatibility
-   * with the function signature in the
-   * general case.
+   * The estimator supports multithreading and splits the cells to
+   * <tt>multithread_info.n_default_threads</tt> (default) threads. The number
+   * of threads to be used in multithreaded mode can be set with the last
+   * parameter of the error estimator.  Multithreading is not presently
+   * implemented for 1d, but we retain the respective parameter for
+   * compatibility with the function signature in the general case.
    */
   template <typename InputVector, class DH>
   static void estimate (const Mapping<1,spacedim>  &mapping,
@@ -621,8 +519,7 @@ public:
                         const types::material_id       material_id = numbers::invalid_material_id);
 
   /**
-   * Calls the @p estimate
-   * function, see above, with
+   * Calls the @p estimate function, see above, with
    * <tt>mapping=MappingQ1<1>()</tt>.
    */
   template <typename InputVector, class DH>
@@ -638,31 +535,17 @@ public:
                         const types::material_id       material_id = numbers::invalid_material_id);
 
   /**
-   * Same function as above, but
-   * accepts more than one solution
-   * vectors and returns one error
-   * vector for each solution
-   * vector. For the reason of
-   * existence of this function,
-   * see the general documentation
-   * of this class.
+   * Same function as above, but accepts more than one solution vectors and
+   * returns one error vector for each solution vector. For the reason of
+   * existence of this function, see the general documentation of this class.
    *
-   * Since we do not want to force
-   * the user of this function to
-   * copy around their solution
-   * vectors, the vector of
-   * solution vectors takes
-   * pointers to the solutions,
-   * rather than being a vector of
-   * vectors. This makes it simpler
-   * to have the solution vectors
-   * somewhere in memory, rather
-   * than to have them collected
-   * somewhere special. (Note that
-   * it is not possible to
-   * construct of vector of
-   * references, so we had to use a
-   * vector of pointers.)
+   * Since we do not want to force the user of this function to copy around
+   * their solution vectors, the vector of solution vectors takes pointers to
+   * the solutions, rather than being a vector of vectors. This makes it
+   * simpler to have the solution vectors somewhere in memory, rather than to
+   * have them collected somewhere special. (Note that it is not possible to
+   * construct of vector of references, so we had to use a vector of
+   * pointers.)
    */
   template <typename InputVector, class DH>
   static void estimate (const Mapping<1,spacedim>          &mapping,
@@ -678,8 +561,7 @@ public:
                         const types::material_id           material_id = numbers::invalid_material_id);
 
   /**
-   * Calls the @p estimate
-   * function, see above, with
+   * Calls the @p estimate function, see above, with
    * <tt>mapping=MappingQ1<1>()</tt>.
    */
   template <typename InputVector, class DH>
@@ -696,10 +578,8 @@ public:
 
 
   /**
-   * Equivalent to the set of functions
-   * above, except that this one takes a
-   * quadrature collection for hp finite
-   * element dof handlers.
+   * Equivalent to the set of functions above, except that this one takes a
+   * quadrature collection for hp finite element dof handlers.
    */
   template <typename InputVector, class DH>
   static void estimate (const Mapping<1,spacedim>      &mapping,
@@ -716,10 +596,8 @@ public:
 
 
   /**
-   * Equivalent to the set of functions
-   * above, except that this one takes a
-   * quadrature collection for hp finite
-   * element dof handlers.
+   * Equivalent to the set of functions above, except that this one takes a
+   * quadrature collection for hp finite element dof handlers.
    */
   template <typename InputVector, class DH>
   static void estimate (const DH                &dof,
@@ -735,10 +613,8 @@ public:
 
 
   /**
-   * Equivalent to the set of functions
-   * above, except that this one takes a
-   * quadrature collection for hp finite
-   * element dof handlers.
+   * Equivalent to the set of functions above, except that this one takes a
+   * quadrature collection for hp finite element dof handlers.
    */
   template <typename InputVector, class DH>
   static void estimate (const Mapping<1,spacedim>          &mapping,
@@ -755,10 +631,8 @@ public:
 
 
   /**
-   * Equivalent to the set of functions
-   * above, except that this one takes a
-   * quadrature collection for hp finite
-   * element dof handlers.
+   * Equivalent to the set of functions above, except that this one takes a
+   * quadrature collection for hp finite element dof handlers.
    */
   template <typename InputVector, class DH>
   static void estimate (const DH                    &dof,

diff --git a/deal.II/source/dofs/dof_accessor.cc b/deal.II/source/dofs/dof_accessor.cc
index 59b8c4833404a969bec960f284849ed4f09998a1..53b67ac4d1997cb95e09b88cd505ee9f63f7434e 100644 (file)
--- a/deal.II/source/dofs/dof_accessor.cc
+++ b/deal.II/source/dofs/dof_accessor.cc
@@ -117,8 +117,8 @@ get_interpolated_dof_values (const InputVector &values,
           typename BaseClass::ExcVectorDoesNotMatch());
 
   if (!this->has_children())
-    // if this cell has no children: simply
-    // return the exact values on this cell
+    // if this cell has no children: simply return the exact values on this
+    // cell
     this->get_dof_values (values, interpolated_values);
   else
     // otherwise clobber them from the children
@@ -128,80 +128,45 @@ get_interpolated_dof_values (const InputVector &values,
 
       interpolated_values = 0;
 
-      // later on we will have to
-      // push the values interpolated
-      // from the child to the mother
-      // cell into the output
-      // vector. unfortunately, there
-      // are two types of elements:
-      // ones where you add up the
-      // contributions from the
-      // different child cells, and
-      // ones where you overwrite.
+      // later on we will have to push the values interpolated from the child
+      // to the mother cell into the output vector. unfortunately, there are
+      // two types of elements: ones where you add up the contributions from
+      // the different child cells, and ones where you overwrite.
       //
-      // an example for the first is
-      // piecewise constant (and
-      // discontinuous) elements,
-      // where we build the value on
-      // the coarse cell by averaging
-      // the values from the cell
-      // (i.e. by adding up a
-      // fraction of the values of
-      // their values)
+      // an example for the first is piecewise constant (and discontinuous)
+      // elements, where we build the value on the coarse cell by averaging
+      // the values from the cell (i.e. by adding up a fraction of the values
+      // of their values)
       //
-      // an example for the latter
-      // are the usual continuous
-      // elements. the value on a
-      // vertex of a coarse cell must
-      // there be the same,
-      // irrespective of the adjacent
-      // cell we are presently on. so
-      // we always overwrite. in
-      // fact, we must, since we
-      // cannot know in advance how
-      // many neighbors there will
-      // be, so there is no way to
-      // compute the average with
-      // fixed factors
+      // an example for the latter are the usual continuous elements. the
+      // value on a vertex of a coarse cell must there be the same,
+      // irrespective of the adjacent cell we are presently on. so we always
+      // overwrite. in fact, we must, since we cannot know in advance how many
+      // neighbors there will be, so there is no way to compute the average
+      // with fixed factors
       //
-      // so we have to find out to
-      // which type this element
-      // belongs. the difficulty is:
-      // the finite element may be a
-      // composed one, so we can only
-      // hope to do this for each
-      // shape function
-      // individually. in fact, there
-      // are even weird finite
-      // elements (for example the
-      // Raviart-Thomas element)
-      // which have shape functions
-      // that are additive (interior
-      // ones) and others that are
-      // overwriting (face degrees of
-      // freedom that need to be
-      // continuous across the
-      // face). to avoid checking
-      // this over and over again, we
-      // do this once now and cache
-      // the results
+      // so we have to find out to which type this element belongs. the
+      // difficulty is: the finite element may be a composed one, so we can
+      // only hope to do this for each shape function individually. in fact,
+      // there are even weird finite elements (for example the Raviart-Thomas
+      // element) which have shape functions that are additive (interior ones)
+      // and others that are overwriting (face degrees of freedom that need to
+      // be continuous across the face). to avoid checking this over and over
+      // again, we do this once now and cache the results
       std::vector<bool> restriction_is_additive (dofs_per_cell);
       for (unsigned int i=0; i<dofs_per_cell; ++i)
         restriction_is_additive[i] = fe.restriction_is_additive(i);
 
       for (unsigned int child=0; child<this->n_children(); ++child)
         {
-          // get the values from the present
-          // child, if necessary by
+          // get the values from the present child, if necessary by
           // interpolation itself
           this->child(child)->get_interpolated_dof_values (values,
                                                            tmp1);
-          // interpolate these to the mother
-          // cell
+          // interpolate these to the mother cell
           fe.get_restriction_matrix(child, this->refinement_case()).vmult (tmp2, tmp1);
 
-          // and add up or set them
-          // in the output vector
+          // and add up or set them in the output vector
           for (unsigned int i=0; i<dofs_per_cell; ++i)
             if (restriction_is_additive[i])
               interpolated_values(i) += tmp2(i);
@@ -232,8 +197,7 @@ set_dof_values_by_interpolation (const Vector<number> &local_values,
           typename BaseClass::ExcVectorDoesNotMatch());
 
   if (!this->has_children())
-    // if this cell has no children: simply
-    // set the values on this cell
+    // if this cell has no children: simply set the values on this cell
     this->set_dof_values (local_values, values);
   else
     // otherwise distribute them to the children
@@ -245,12 +209,8 @@ set_dof_values_by_interpolation (const Vector<number> &local_values,
           Assert (this->child(child)->get_fe().dofs_per_cell == dofs_per_cell,
                   ExcNotImplemented());
 
-          // prolong the given data
-          // to the present
-          // cell. FullMatrix only
-          // wants us to call vmult
-          // if the matrix size is
-          // actually non-zero, so
+          // prolong the given data to the present cell. FullMatrix only wants
+          // us to call vmult if the matrix size is actually non-zero, so
           // check that case
           if (tmp.size() > 0)
             {

diff --git a/deal.II/source/numerics/CMakeLists.txt b/deal.II/source/numerics/CMakeLists.txt
index db675fcd7a0d0733fd8f300a4435c1880bf127c2..7ab970c44704b08c2d51d3d83fd19e5093d38a93 100644 (file)
--- a/deal.II/source/numerics/CMakeLists.txt
+++ b/deal.II/source/numerics/CMakeLists.txt
@@ -29,6 +29,7 @@ SET(_src
   derivative_approximation.cc
   dof_output_operator.cc
   error_estimator.cc
+  error_estimator_1d.cc
   fe_field_function.cc
   matrix_tools.cc
   histogram.cc
@@ -55,6 +56,7 @@ SET(_inst
   derivative_approximation.inst.in
   dof_output_operator.inst.in
   error_estimator.inst.in
+  error_estimator_1d.inst.in
   fe_field_function.inst.in
   matrix_tools.inst.in
   point_value_history.inst.in

diff --git a/deal.II/source/numerics/error_estimator.cc b/deal.II/source/numerics/error_estimator.cc
index 46771ecc59e28fb605701816613fa113ff534736..3fb633c4a5e269bdbe40554145ea8fc14e545ea3 100644 (file)
--- a/deal.II/source/numerics/error_estimator.cc
+++ b/deal.II/source/numerics/error_estimator.cc
@@ -51,24 +51,14 @@ DEAL_II_NAMESPACE_OPEN
 
 namespace
 {
-  inline
-  double sqr (const double x)
-  {
-    return x*x;
-  }
-
-
   template <typename CellIterator>
   inline
   void advance_by_n (CellIterator &cell,
                      const unsigned int n)
   {
-    // store a pointer to the end
-    // iterator, since we can't get at
-    // it any more once cell is already
-    // the end iterator (in that case
-    // dereferencing cell-> triggers an
-    // assertion)
+    // store a pointer to the end iterator, since we can't get at it any more
+    // once cell is already the end iterator (in that case dereferencing
+    // cell-> triggers an assertion)
     const CellIterator endc = cell->get_dof_handler().end();
     for (unsigned int t=0; ((t<n) && (cell!=endc)); ++t, ++cell)
       ;
@@ -81,60 +71,31 @@ namespace internal
   namespace
   {
     /**
-     * All small temporary data
-     * objects that are needed once
-     * per thread by the several
-     * functions of the error
-     * estimator are gathered in this
-     * struct. The reason for this
-     * structure is mainly that we
-     * have a number of functions
-     * that operate on cells or faces
-     * and need a number of small
-     * temporary data objects. Since
-     * these functions may run in
-     * parallel, we cannot make these
-     * objects member variables of
-     * the enclosing class. On the
-     * other hand, declaring them
-     * locally in each of these
-     * functions would require their
-     * reallocating every time we
-     * visit the next cell or face,
-     * which we found can take a
-     * significant amount of time if
-     * it happens often even in the
-     * single threaded case (10-20
-     * per cent in our measurements);
-     * however, most importantly,
-     * memory allocation requires
-     * synchronisation in
-     * multithreaded mode. While that
-     * is done by the C++ library and
-     * has not to be handcoded, it
-     * nevertheless seriously damages
-     * the ability to efficiently run
-     * the functions of this class in
-     * parallel, since they are quite
-     * often blocked by these
-     * synchronisation points,
-     * slowing everything down by a
-     * factor of two or three.
+     * All small temporary data objects that are needed once per thread by the
+     * several functions of the error estimator are gathered in this
+     * struct. The reason for this structure is mainly that we have a number
+     * of functions that operate on cells or faces and need a number of small
+     * temporary data objects. Since these functions may run in parallel, we
+     * cannot make these objects member variables of the enclosing class. On
+     * the other hand, declaring them locally in each of these functions would
+     * require their reallocating every time we visit the next cell or face,
+     * which we found can take a significant amount of time if it happens
+     * often even in the single threaded case (10-20 per cent in our
+     * measurements); however, most importantly, memory allocation requires
+     * synchronisation in multithreaded mode. While that is done by the C++
+     * library and has not to be handcoded, it nevertheless seriously damages
+     * the ability to efficiently run the functions of this class in parallel,
+     * since they are quite often blocked by these synchronisation points,
+     * slowing everything down by a factor of two or three.
      *
-     * Thus, every thread gets an
-     * instance of this class to work
-     * with and needs not allocate
-     * memory itself, or synchronise
-     * with other threads.
+     * Thus, every thread gets an instance of this class to work with and
+     * needs not allocate memory itself, or synchronise with other threads.
      *
-     * The sizes of the arrays are
-     * initialized with the maximal number of
-     * entries necessary for the hp
-     * case. Within the loop over individual
-     * cells, we then resize the arrays as
-     * necessary. Since for std::vector
-     * resizing to a smaller size doesn't
-     * imply memory allocation, this is fast.
+     * The sizes of the arrays are initialized with the maximal number of
+     * entries necessary for the hp case. Within the loop over individual
+     * cells, we then resize the arrays as necessary. Since for std::vector
+     * resizing to a smaller size doesn't imply memory allocation, this is
+     * fast.
      */
     template <class DH>
     struct ParallelData
@@ -148,14 +109,12 @@ namespace internal
       const dealii::hp::FECollection<dim,spacedim> finite_element;
 
       /**
-       * The quadrature formulas to be used for
-       * the faces.
+       * The quadrature formulas to be used for the faces.
        */
       const dealii::hp::QCollection<dim-1> face_quadratures;
 
       /**
-       * FEFaceValues objects to integrate over
-       * the faces of the current and
+       * FEFaceValues objects to integrate over the faces of the current and
        * potentially of neighbor cells.
        */
       dealii::hp::FEFaceValues<dim,spacedim>    fe_face_values_cell;
@@ -163,38 +122,22 @@ namespace internal
       dealii::hp::FESubfaceValues<dim,spacedim> fe_subface_values;
 
       /**
-       * A vector to store the jump
-       * of the normal vectors in
-       * the quadrature points for
-       * each of the solution
-       * vectors (i.e. a temporary
-       * value). This vector is not
-       * allocated inside the
-       * functions that use it, but
-       * rather globally, since
-       * memory allocation is slow,
-       * in particular in presence
-       * of multiple threads where
-       * synchronisation makes
-       * things even slower.
+       * A vector to store the jump of the normal vectors in the quadrature
+       * points for each of the solution vectors (i.e. a temporary
+       * value). This vector is not allocated inside the functions that use
+       * it, but rather globally, since memory allocation is slow, in
+       * particular in presence of multiple threads where synchronisation
+       * makes things even slower.
        */
       std::vector<std::vector<std::vector<double> > > phi;
 
       /**
-       * A vector for the gradients of
-       * the finite element function
-       * on one cell
+       * A vector for the gradients of the finite element function on one cell
        *
-       * Let psi be a short name
-       * for <tt>a grad u_h</tt>, where
-       * the third index be the
-       * component of the finite
-       * element, and the second
-       * index the number of the
-       * quadrature point. The
-       * first index denotes the
-       * index of the solution
-       * vector.
+       * Let psi be a short name for <tt>a grad u_h</tt>, where the third
+       * index be the component of the finite element, and the second index
+       * the number of the quadrature point. The first index denotes the index
+       * of the solution vector.
        */
       std::vector<std::vector<std::vector<Tensor<1,spacedim> > > > psi;
 
@@ -204,43 +147,35 @@ namespace internal
       std::vector<std::vector<std::vector<Tensor<1,spacedim> > > > neighbor_psi;
 
       /**
-       * The normal vectors of the finite
-       * element function on one face
+       * The normal vectors of the finite element function on one face
        */
       std::vector<Point<spacedim> > normal_vectors;
 
       /**
-       * Two arrays needed for the
-       * values of coefficients in
-       * the jumps, if they are
-       * given.
+       * Two arrays needed for the values of coefficients in the jumps, if
+       * they are given.
        */
       std::vector<double>                  coefficient_values1;
       std::vector<dealii::Vector<double> > coefficient_values;
 
       /**
-       * Array for the products of
-       * Jacobian determinants and
-       * weights of quadraturs
-       * points.
+       * Array for the products of Jacobian determinants and weights of
+       * quadraturs points.
        */
       std::vector<double>          JxW_values;
 
       /**
-       * The subdomain id we are to care
-       * for.
+       * The subdomain id we are to care for.
        */
       const types::subdomain_id subdomain_id;
       /**
-       * The material id we are to care
-       * for.
+       * The material id we are to care for.
        */
       const types::material_id material_id;
 
       /**
-       * Some more references to input data to
-       * the KellyErrorEstimator::estimate()
-       * function.
+       * Some more references to input data to the
+       * KellyErrorEstimator::estimate() function.
        */
       const typename FunctionMap<spacedim>::type *neumann_bc;
       const ComponentMask                component_mask;
@@ -262,10 +197,9 @@ namespace internal
                     const Function<spacedim>                   *coefficients);
 
       /**
-       * Resize the arrays so that they fit the
-       * number of quadrature points associated
-       * with the given finite element index
-       * into the hp collections.
+       * Resize the arrays so that they fit the number of quadrature points
+       * associated with the given finite element index into the hp
+       * collections.
        */
       void resize (const unsigned int active_fe_index);
     };
@@ -363,11 +297,9 @@ namespace internal
 
 
     /**
-     * Copy data from the
-     * local_face_integrals map of a single
-     * ParallelData object into a global such
-     * map. This is the copier stage of a
-     * WorkStream pipeline.
+     * Copy data from the local_face_integrals map of a single ParallelData
+     * object into a global such map. This is the copier stage of a WorkStream
+     * pipeline.
      */
     template <class DH>
     void
@@ -375,16 +307,14 @@ namespace internal
                           std::map<typename DH::face_iterator,std::vector<double> > &face_integrals)
     {
 
-      // now copy locally computed elements
-      // into the global map
+      // now copy locally computed elements into the global map
       for (typename std::map<typename DH::face_iterator,std::vector<double> >::const_iterator
            p=local_face_integrals.begin();
            p!=local_face_integrals.end();
            ++p)
         {
-          // double check that the
-          // element does not already
-          // exists in the global map
+          // double check that the element does not already exists in the
+          // global map
           Assert (face_integrals.find (p->first) == face_integrals.end(),
                   ExcInternalError());
 
@@ -400,105 +330,40 @@ namespace internal
 
 
     /**
-     * Actually do the computation on
-     * a face which has no hanging
-     * nodes (it is regular), i.e.
-     * either on the other side there
-     * is nirvana (face is at
-     * boundary), or the other side's
-     * refinement level is the same
-     * as that of this side, then
-     * handle the integration of
-     * these both cases together.
+     * Actually do the computation based on the evaluated gradients in
+     * ParallelData.
      */
-    template <typename InputVector, class DH>
-    void
-    integrate_over_regular_face (const std::vector<const InputVector *>   &solutions,
-                                 ParallelData<DH>                        &parallel_data,
-                                 std::map<typename DH::face_iterator,std::vector<double> > &local_face_integrals,
-                                 const typename DH::active_cell_iterator &cell,
-                                 const unsigned int                       face_no,
-                                 dealii::hp::FEFaceValues<DH::dimension, DH::space_dimension> &fe_face_values_cell,
-                                 dealii::hp::FEFaceValues<DH::dimension, DH::space_dimension> &fe_face_values_neighbor)
+    template <class DH>
+    std::vector<double>
+    integrate_over_face (ParallelData<DH>                        &parallel_data,
+                         const typename DH::face_iterator        &face,
+                         dealii::hp::FEFaceValues<DH::dimension, DH::space_dimension> &fe_face_values_cell)
     {
       const unsigned int dim = DH::dimension;
 
-      const typename DH::face_iterator face = cell->face(face_no);
-      const unsigned int n_q_points         = parallel_data.face_quadratures[cell->active_fe_index()].size(),
+      const unsigned int n_q_points         = parallel_data.psi[0].size(),
                          n_components       = parallel_data.finite_element.n_components(),
-                         n_solution_vectors = solutions.size();
-
-
-      // initialize data of the restriction
-      // of this cell to the present face
-      fe_face_values_cell.reinit (cell, face_no,
-                                  cell->active_fe_index());
+                         n_solution_vectors = parallel_data.psi.size();
 
-      // get gradients of the finite element
-      // function on this cell
-      for (unsigned int n=0; n<n_solution_vectors; ++n)
-        fe_face_values_cell.get_present_fe_values()
-        .get_function_gradients (*solutions[n], parallel_data.psi[n]);
-
-      // now compute over the other side of
-      // the face
       if (face->at_boundary() == false)
-        // internal face; integrate jump
-        // of gradient across this face
         {
-          Assert (cell->neighbor(face_no).state() == IteratorState::valid,
-                  ExcInternalError());
-
-          const typename DH::active_cell_iterator neighbor = cell->neighbor(face_no);
-
-          // find which number the
-          // current face has relative to
-          // the neighboring cell
-          const unsigned int neighbor_neighbor
-            = cell->neighbor_of_neighbor (face_no);
-          Assert (neighbor_neighbor<GeometryInfo<dim>::faces_per_cell,
-                  ExcInternalError());
-
-          // get restriction of finite element
-          // function of @p{neighbor} to the
-          // common face. in the hp case, use the
-          // quadrature formula that matches the
-          // one we would use for the present
-          // cell
-          fe_face_values_neighbor.reinit (neighbor, neighbor_neighbor,
-                                          cell->active_fe_index());
-
-          // get gradients on neighbor cell
+          // compute the jump in the gradients
           for (unsigned int n=0; n<n_solution_vectors; ++n)
-            {
-              fe_face_values_neighbor.get_present_fe_values()
-              .get_function_gradients (*solutions[n],
-                                       parallel_data.neighbor_psi[n]);
-
-              // compute the jump in the gradients
-              for (unsigned int component=0; component<n_components; ++component)
-                for (unsigned int p=0; p<n_q_points; ++p)
-                  parallel_data.psi[n][p][component] -= parallel_data.neighbor_psi[n][p][component];
-            }
+            for (unsigned int component=0; component<n_components; ++component)
+              for (unsigned int p=0; p<n_q_points; ++p)
+                parallel_data.psi[n][p][component] -= parallel_data.neighbor_psi[n][p][component];
         }
 
-
       // now psi contains the following:
       // - for an internal face, psi=[grad u]
-      // - for a neumann boundary face,
-      //   psi=grad u
-      // each component being the
-      // mentioned value at one of the
-      // quadrature points
-
-      // next we have to multiply this with
-      // the normal vector. Since we have
-      // taken the difference of gradients
-      // for internal faces, we may chose
-      // the normal vector of one cell,
-      // taking that of the neighbor
-      // would only change the sign. We take
-      // the outward normal.
+      // - for a neumann boundary face, psi=grad u
+      // each component being the mentioned value at one of the quadrature
+      // points
+
+      // next we have to multiply this with the normal vector. Since we have
+      // taken the difference of gradients for internal faces, we may chose
+      // the normal vector of one cell, taking that of the neighbor would only
+      // change the sign. We take the outward normal.
 
       parallel_data.normal_vectors =
         fe_face_values_cell.get_present_fe_values().get_normal_vectors();
@@ -510,14 +375,13 @@ namespace internal
               = (parallel_data.psi[n][point][component] *
                  parallel_data.normal_vectors[point]);
 
-      // if a coefficient was given: use that
-      // to scale the jump in the gradient
+      // if a coefficient was given: use that to scale the jump in the
+      // gradient
       if (parallel_data.coefficients != 0)
         {
           // scalar coefficient
           if (parallel_data.coefficients->n_components == 1)
             {
-
               parallel_data.coefficients
               ->value_list (fe_face_values_cell.get_present_fe_values()
                             .get_quadrature_points(),
@@ -545,18 +409,15 @@ namespace internal
 
 
       if (face->at_boundary() == true)
-        // neumann boundary face. compute
-        // difference between normal
-        // derivative and boundary function
+        // neumann boundary face. compute difference between normal derivative
+        // and boundary function
         {
           const types::boundary_id boundary_indicator = face->boundary_indicator();
 
           Assert (parallel_data.neumann_bc->find(boundary_indicator) !=
                   parallel_data.neumann_bc->end(),
                   ExcInternalError ());
-          // get the values of the boundary
-          // function at the quadrature
-          // points
+          // get the values of the boundary function at the quadrature points
           if (n_components == 1)
             {
               std::vector<double> g(n_q_points);
@@ -587,26 +448,96 @@ namespace internal
 
       // now phi contains the following:
       // - for an internal face, phi=[a du/dn]
-      // - for a neumann boundary face,
-      //   phi=a du/dn-g
-      // each component being the
-      // mentioned value at one of the
-      // quadrature points
+      // - for a neumann boundary face, phi=a du/dn-g
+      // each component being the mentioned value at one of the quadrature
+      // points
 
       parallel_data.JxW_values
         = fe_face_values_cell.get_present_fe_values().get_JxW_values();
 
-      // take the square of the phi[i]
-      // for integration, and sum up
+      // take the square of the phi[i] for integration, and sum up
       std::vector<double> face_integral (n_solution_vectors, 0);
       for (unsigned int n=0; n<n_solution_vectors; ++n)
         for (unsigned int component=0; component<n_components; ++component)
           if (parallel_data.component_mask[component] == true)
             for (unsigned int p=0; p<n_q_points; ++p)
-              face_integral[n] += dealii::sqr(parallel_data.phi[n][p][component]) *
+              face_integral[n] += Utilities::fixed_power<2>(parallel_data.phi[n][p][component]) *
                                   parallel_data.JxW_values[p];
 
-      local_face_integrals[face] = face_integral;
+      return face_integral;
+    }
+
+
+
+    /**
+     * Actually do the computation on a face which has no hanging nodes (it is
+     * regular), i.e. either on the other side there is nirvana (face is at
+     * boundary), or the other side's refinement level is the same as that of
+     * this side, then handle the integration of these both cases together.
+     */
+    template <typename InputVector, class DH>
+    void
+    integrate_over_regular_face (const std::vector<const InputVector *>   &solutions,
+                                 ParallelData<DH>                        &parallel_data,
+                                 std::map<typename DH::face_iterator,std::vector<double> > &local_face_integrals,
+                                 const typename DH::active_cell_iterator &cell,
+                                 const unsigned int                       face_no,
+                                 dealii::hp::FEFaceValues<DH::dimension, DH::space_dimension> &fe_face_values_cell,
+                                 dealii::hp::FEFaceValues<DH::dimension, DH::space_dimension> &fe_face_values_neighbor)
+    {
+      const unsigned int dim = DH::dimension;
+
+      const typename DH::face_iterator face = cell->face(face_no);
+      const unsigned int n_q_points         = parallel_data.face_quadratures[cell->active_fe_index()].size(),
+                         n_components       = parallel_data.finite_element.n_components(),
+                         n_solution_vectors = solutions.size();
+
+
+      // initialize data of the restriction
+      // of this cell to the present face
+      fe_face_values_cell.reinit (cell, face_no,
+                                  cell->active_fe_index());
+
+      // get gradients of the finite element
+      // function on this cell
+      for (unsigned int n=0; n<n_solution_vectors; ++n)
+        fe_face_values_cell.get_present_fe_values()
+        .get_function_gradients (*solutions[n], parallel_data.psi[n]);
+
+      // now compute over the other side of the face
+      if (face->at_boundary() == false)
+        // internal face; integrate jump of gradient across this face
+        {
+          Assert (cell->neighbor(face_no).state() == IteratorState::valid,
+                  ExcInternalError());
+
+          const typename DH::active_cell_iterator neighbor = cell->neighbor(face_no);
+
+          // find which number the current face has relative to the
+          // neighboring cell
+          const unsigned int neighbor_neighbor
+            = cell->neighbor_of_neighbor (face_no);
+          Assert (neighbor_neighbor<GeometryInfo<dim>::faces_per_cell,
+                  ExcInternalError());
+
+          // get restriction of finite element function of @p{neighbor} to the
+          // common face. in the hp case, use the quadrature formula that
+          // matches the one we would use for the present cell
+          fe_face_values_neighbor.reinit (neighbor, neighbor_neighbor,
+                                          cell->active_fe_index());
+
+          // get gradients on neighbor cell
+          for (unsigned int n=0; n<n_solution_vectors; ++n)
+            {
+              fe_face_values_neighbor.get_present_fe_values()
+              .get_function_gradients (*solutions[n],
+                                       parallel_data.neighbor_psi[n]);
+            }
+        }
+
+      // now go to the generic function that does all the other things
+      local_face_integrals[face] = 
+        integrate_over_face (parallel_data, face, fe_face_values_cell);
     }
 
 
@@ -624,8 +555,8 @@ namespace internal
      */
     template <typename InputVector, class DH>
     void
-    integrate_over_irregular_face (const std::vector<const InputVector *>      &solutions,
-                                   ParallelData<DH>                           &parallel_data,
+    integrate_over_irregular_face (const std::vector<const InputVector *>   &solutions,
+                                   ParallelData<DH>                         &parallel_data,
                                    std::map<typename DH::face_iterator,std::vector<double> > &local_face_integrals,
                                    const typename DH::active_cell_iterator    &cell,
                                    const unsigned int                          face_no,
@@ -643,20 +574,15 @@ namespace internal
 
       Assert (neighbor.state() == IteratorState::valid, ExcInternalError());
       Assert (face->has_children(), ExcInternalError());
-      // set up a vector of the gradients
-      // of the finite element function
-      // on this cell at the quadrature
-      // points
+
+      // set up a vector of the gradients of the finite element function on
+      // this cell at the quadrature points
       //
-      // let psi be a short name for
-      // [a grad u_h], where the second
-      // index be the component of the
-      // finite element, and the first
-      // index the number of the
-      // quadrature point
-
-      // store which number @p{cell} has
-      // in the list of neighbors of
+      // let psi be a short name for [a grad u_h], where the second index be
+      // the component of the finite element, and the first index the number
+      // of the quadrature point
+
+      // store which number @p{cell} has in the list of neighbors of
       // @p{neighbor}
       const unsigned int neighbor_neighbor
         = cell->neighbor_of_neighbor (face_no);
@@ -666,138 +592,38 @@ namespace internal
       // loop over all subfaces
       for (unsigned int subface_no=0; subface_no<face->n_children(); ++subface_no)
         {
-          // get an iterator pointing to the
-          // cell behind the present subface
+          // get an iterator pointing to the cell behind the present subface
           const typename DH::active_cell_iterator neighbor_child
             = cell->neighbor_child_on_subface (face_no, subface_no);
           Assert (!neighbor_child->has_children(),
                   ExcInternalError());
 
-          // restrict the finite element
-          // on the present cell to the
-          // subface
+          // restrict the finite element on the present cell to the subface
           fe_subface_values.reinit (cell, face_no, subface_no,
                                     cell->active_fe_index());
 
-          // restrict the finite element
-          // on the neighbor cell to the
-          // common @p{subface}.
+          // restrict the finite element on the neighbor cell to the common
+          // @p{subface}.
           fe_face_values.reinit (neighbor_child, neighbor_neighbor,
                                  cell->active_fe_index());
 
-          // store the gradient of the
-          // solution in psi
+          // store the gradient of the solution in psi
           for (unsigned int n=0; n<n_solution_vectors; ++n)
             fe_subface_values.get_present_fe_values()
             .get_function_gradients (*solutions[n], parallel_data.psi[n]);
 
-          // store the gradient from the
-          // neighbor's side in
-          // @p{neighbor_psi}
+          // store the gradient from the neighbor's side in @p{neighbor_psi}
           for (unsigned int n=0; n<n_solution_vectors; ++n)
             fe_face_values.get_present_fe_values()
             .get_function_gradients (*solutions[n], parallel_data.neighbor_psi[n]);
 
-          // compute the jump in the gradients
-          for (unsigned int n=0; n<n_solution_vectors; ++n)
-            for (unsigned int component=0; component<n_components; ++component)
-              for (unsigned int p=0; p<n_q_points; ++p)
-                parallel_data.psi[n][p][component] -=
-                  parallel_data.neighbor_psi[n][p][component];
-
-          // note that unlike for the
-          // case of regular faces
-          // (treated in the other
-          // function of this class), we
-          // have not to take care of
-          // boundary faces here, since
-          // they always are regular.
-
-          // next we have to multiply this with
-          // the normal vector. Since we have
-          // taken the difference of gradients
-          // for internal faces, we may chose
-          // the normal vector of one cell,
-          // taking that of the neighbor
-          // would only change the sign. We take
-          // the outward normal.
-          //
-          // let phi be the name of the integrand
-
-          parallel_data.normal_vectors
-            = fe_face_values.get_present_fe_values().get_normal_vectors();
-
-
-          for (unsigned int n=0; n<n_solution_vectors; ++n)
-            for (unsigned int component=0; component<n_components; ++component)
-              for (unsigned int point=0; point<n_q_points; ++point)
-                parallel_data.phi[n][point][component] = (parallel_data.psi[n][point][component]*
-                                                          parallel_data.normal_vectors[point]);
-
-          // if a coefficient was given: use that
-          // to scale the jump in the gradient
-          if (parallel_data.coefficients != 0)
-            {
-              // scalar coefficient
-              if (parallel_data.coefficients->n_components == 1)
-                {
-                  parallel_data.coefficients
-                  ->value_list (fe_face_values.get_present_fe_values()
-                                .get_quadrature_points(),
-                                parallel_data.coefficient_values1);
-                  for (unsigned int n=0; n<n_solution_vectors; ++n)
-                    for (unsigned int component=0; component<n_components; ++component)
-                      for (unsigned int point=0; point<n_q_points; ++point)
-                        parallel_data.phi[n][point][component] *=
-                          parallel_data.coefficient_values1[point];
-                }
-              else
-                // vector-valued coefficient
-                {
-                  parallel_data.coefficients
-                  ->vector_value_list (fe_face_values.get_present_fe_values()
-                                       .get_quadrature_points(),
-                                       parallel_data.coefficient_values);
-                  for (unsigned int n=0; n<n_solution_vectors; ++n)
-                    for (unsigned int component=0; component<n_components; ++component)
-                      for (unsigned int point=0; point<n_q_points; ++point)
-                        parallel_data.phi[n][point][component] *=
-                          parallel_data.coefficient_values[point](component);
-                }
-            }
-
-          // get the weights for the
-          // integration. note that it
-          // does not matter whether we
-          // take the JxW values from the
-          // fe_face_values or the
-          // fe_subface_values, as the
-          // first is on the small
-          // neighbor cell, while the
-          // latter is on the refined
-          // face of the big cell here
-          parallel_data.JxW_values
-            = fe_face_values.get_present_fe_values().get_JxW_values();
-
-          // take the square of the phi[i]
-          // for integration, and sum up
-          std::vector<double> face_integral (n_solution_vectors, 0);
-          for (unsigned int n=0; n<n_solution_vectors; ++n)
-            for (unsigned int component=0; component<n_components; ++component)
-              if (parallel_data.component_mask[component] == true)
-                for (unsigned int p=0; p<n_q_points; ++p)
-                  face_integral[n] += dealii::sqr(parallel_data.phi[n][p][component]) *
-                                      parallel_data.JxW_values[p];
-
-          local_face_integrals[neighbor_child->face(neighbor_neighbor)]
-            = face_integral;
+          // call generic evaluate function
+          local_face_integrals[neighbor_child->face(neighbor_neighbor)] =
+            integrate_over_face (parallel_data, face, fe_face_values);
         }
 
-
-      // finally loop over all subfaces to
-      // collect the contributions of the
-      // subfaces and store them with the
-      // mother face
+      // finally loop over all subfaces to collect the contributions of the
+      // subfaces and store them with the mother face
       std::vector<double> sum (n_solution_vectors, 0);
       for (unsigned int subface_no=0; subface_no<face->n_children(); ++subface_no)
         {
@@ -816,29 +642,26 @@ namespace internal
 
 
     /**
-     * Computate the error on the faces of a
-     * single cell.
+     * Computate the error on the faces of a single cell.
      *
-     * This function is only needed
-     * in two or three dimensions.
-     * The error estimator in one
-     * dimension is implemented
-     * seperatly.
+     * This function is only needed in two or three dimensions.  The error
+     * estimator in one dimension is implemented separately.
      */
-    template <int dim, int spacedim, typename InputVector, class DH>
+    template <typename InputVector, class DH>
     void
     estimate_one_cell (const typename DH::active_cell_iterator &cell,
                        ParallelData<DH>                    &parallel_data,
                        std::map<typename DH::face_iterator,std::vector<double> > &local_face_integrals,
                        const std::vector<const InputVector *> &solutions)
     {
+      const unsigned int dim = DH::dimension;
+      const unsigned int spacedim = DH::dimension;
       const unsigned int n_solution_vectors = solutions.size();
 
       const types::subdomain_id subdomain_id = parallel_data.subdomain_id;
       const unsigned int material_id  = parallel_data.material_id;
 
-      // empty our own copy of the local face
-      // integrals
+      // empty our own copy of the local face integrals
       local_face_integrals.clear();
 
       // loop over all faces of this cell
@@ -848,20 +671,11 @@ namespace internal
           const typename DH::face_iterator
           face=cell->face(face_no);
 
-          // make sure we do work
-          // only once: this face
-          // may either be regular
-          // or irregular. if it is
-          // regular and has a
-          // neighbor, then we
-          // visit the face twice,
-          // once from every
-          // side. let the one with
-          // the lower index do the
-          // work. if it is at the
-          // boundary, or if the
-          // face is irregular,
-          // then do the work below
+          // make sure we do work only once: this face may either be regular
+          // or irregular. if it is regular and has a neighbor, then we visit
+          // the face twice, once from every side. let the one with the lower
+          // index do the work. if it is at the boundary, or if the face is
+          // irregular, then do the work below
           if ((face->has_children() == false) &&
               !cell->at_boundary(face_no) &&
               (!cell->neighbor_is_coarser(face_no) &&
@@ -870,24 +684,17 @@ namespace internal
                  cell->neighbor(face_no)->level() < cell->level()))))
             continue;
 
-          // if the neighboring cell is less
-          // refined than the present one,
-          // then do nothing since we
-          // integrate over the subfaces when
-          // we visit the coarse cells.
+          // if the neighboring cell is less refined than the present one,
+          // then do nothing since we integrate over the subfaces when we
+          // visit the coarse cells.
           if (face->at_boundary() == false)
             if (cell->neighbor_is_coarser(face_no))
               continue;
 
-          // if this face is part of the
-          // boundary but not of the neumann
-          // boundary -> nothing to
-          // do. However, to make things
-          // easier when summing up the
-          // contributions of the faces of
-          // cells, we enter this face into
-          // the list of faces with
-          // contribution zero.
+          // if this face is part of the boundary but not of the neumann
+          // boundary -> nothing to do. However, to make things easier when
+          // summing up the contributions of the faces of cells, we enter this
+          // face into the list of faces with contribution zero.
           if (face->at_boundary()
               &&
               (parallel_data.neumann_bc->find(face->boundary_indicator()) ==
@@ -898,12 +705,9 @@ namespace internal
               continue;
             }
 
-          // finally: note that we only have
-          // to do something if either the
-          // present cell is on the subdomain
-          // we care for (and the same for
-          // material_id), or if one of the
-          // neighbors behind the face is on
+          // finally: note that we only have to do something if either the
+          // present cell is on the subdomain we care for (and the same for
+          // material_id), or if one of the neighbors behind the face is on
           // the subdomain we care for
           if ( ! ( ((subdomain_id == numbers::invalid_subdomain_id)
                     ||
@@ -913,11 +717,8 @@ namespace internal
                     ||
                     (cell->material_id() == material_id))) )
             {
-              // ok, cell is unwanted, but
-              // maybe its neighbor behind
-              // the face we presently work
-              // on? oh is there a face at
-              // all?
+              // ok, cell is unwanted, but maybe its neighbor behind the face
+              // we presently work on? oh is there a face at all?
               if (face->at_boundary())
                 continue;
 
@@ -950,32 +751,24 @@ namespace internal
                       }
                 }
 
-              // so if none of the neighbors
-              // cares for this subdomain or
-              // material either, then try
-              // next face
+              // so if none of the neighbors cares for this subdomain or
+              // material either, then try next face
               if (care_for_cell == false)
                 continue;
             }
 
-          // so now we know that we care for
-          // this face, let's do something
-          // about it. first re-size the
-          // arrays we may use to the correct
+          // so now we know that we care for this face, let's do something
+          // about it. first re-size the arrays we may use to the correct
           // size:
           parallel_data.resize (cell->active_fe_index());
 
 
           // then do the actual integration
           if (face->has_children() == false)
-            // if the face is a regular one,
-            // i.e.  either on the other side
-            // there is nirvana (face is at
-            // boundary), or the other side's
-            // refinement level is the same
-            // as that of this side, then
-            // handle the integration of
-            // these both cases together
+            // if the face is a regular one, i.e.  either on the other side
+            // there is nirvana (face is at boundary), or the other side's
+            // refinement level is the same as that of this side, then handle
+            // the integration of these both cases together
             integrate_over_regular_face (solutions,
                                          parallel_data,
                                          local_face_integrals,
@@ -984,10 +777,8 @@ namespace internal
                                          parallel_data.fe_face_values_neighbor);
 
           else
-            // otherwise we need to do some
-            // special computations which do
-            // not fit into the framework of
-            // the above function
+            // otherwise we need to do some special computations which do not
+            // fit into the framework of the above function
             integrate_over_irregular_face (solutions,
                                            parallel_data,
                                            local_face_integrals,
@@ -1003,418 +794,6 @@ namespace internal
 
 
 
-
-
-template <int spacedim>
-template <typename InputVector, class DH>
-void
-KellyErrorEstimator<1,spacedim>::
-estimate (const Mapping<1,spacedim>      &mapping,
-          const DH   &dof_handler,
-          const Quadrature<0> &quadrature,
-          const typename FunctionMap<spacedim>::type &neumann_bc,
-          const InputVector       &solution,
-          Vector<float>           &error,
-          const ComponentMask &component_mask,
-          const Function<spacedim>     *coefficients,
-          const unsigned int       n_threads,
-          const types::subdomain_id subdomain_id,
-          const types::material_id       material_id)
-{
-  // just pass on to the other function
-  const std::vector<const InputVector *> solutions (1, &solution);
-  std::vector<Vector<float>*>              errors (1, &error);
-  estimate (mapping, dof_handler, quadrature, neumann_bc, solutions, errors,
-            component_mask, coefficients, n_threads, subdomain_id, material_id);
-}
-
-
-
-template <int spacedim>
-template <typename InputVector, class DH>
-void
-KellyErrorEstimator<1,spacedim>::
-estimate (const DH   &dof_handler,
-          const Quadrature<0> &quadrature,
-          const typename FunctionMap<spacedim>::type &neumann_bc,
-          const InputVector       &solution,
-          Vector<float>           &error,
-          const ComponentMask &component_mask,
-          const Function<spacedim>     *coefficients,
-          const unsigned int       n_threads,
-          const types::subdomain_id subdomain_id,
-          const types::material_id       material_id)
-{
-  estimate(StaticMappingQ1<1,spacedim>::mapping, dof_handler, quadrature, neumann_bc, solution,
-           error, component_mask, coefficients, n_threads, subdomain_id, material_id);
-}
-
-
-
-template <int spacedim>
-template <typename InputVector, class DH>
-void
-KellyErrorEstimator<1,spacedim>::
-estimate (const DH   &dof_handler,
-          const Quadrature<0> &quadrature,
-          const typename FunctionMap<spacedim>::type &neumann_bc,
-          const std::vector<const InputVector *> &solutions,
-          std::vector<Vector<float>*> &errors,
-          const ComponentMask &component_mask,
-          const Function<spacedim>     *coefficients,
-          const unsigned int       n_threads,
-          const types::subdomain_id subdomain_id,
-          const types::material_id       material_id)
-{
-  estimate(StaticMappingQ1<1,spacedim>::mapping, dof_handler, quadrature, neumann_bc, solutions,
-           errors, component_mask, coefficients, n_threads, subdomain_id, material_id);
-}
-
-
-
-template <int spacedim>
-template <typename InputVector, class DH>
-void
-KellyErrorEstimator<1,spacedim>::
-estimate (const Mapping<1,spacedim>      &mapping,
-          const DH   &dof_handler,
-          const hp::QCollection<0> &quadrature,
-          const typename FunctionMap<spacedim>::type &neumann_bc,
-          const InputVector       &solution,
-          Vector<float>           &error,
-          const ComponentMask &component_mask,
-          const Function<spacedim>     *coefficients,
-          const unsigned int       n_threads,
-          const types::subdomain_id subdomain_id,
-          const types::material_id       material_id)
-{
-  // just pass on to the other function
-  const std::vector<const InputVector *> solutions (1, &solution);
-  std::vector<Vector<float>*>              errors (1, &error);
-  estimate (mapping, dof_handler, quadrature, neumann_bc, solutions, errors,
-            component_mask, coefficients, n_threads, subdomain_id, material_id);
-}
-
-
-template <int spacedim>
-template <typename InputVector, class DH>
-void
-KellyErrorEstimator<1,spacedim>::
-estimate (const DH   &dof_handler,
-          const hp::QCollection<0> &quadrature,
-          const typename FunctionMap<spacedim>::type &neumann_bc,
-          const InputVector       &solution,
-          Vector<float>           &error,
-          const ComponentMask &component_mask,
-          const Function<spacedim>     *coefficients,
-          const unsigned int       n_threads,
-          const types::subdomain_id subdomain_id,
-          const types::material_id       material_id)
-{
-  estimate(StaticMappingQ1<1,spacedim>::mapping, dof_handler, quadrature, neumann_bc, solution,
-           error, component_mask, coefficients, n_threads, subdomain_id, material_id);
-}
-
-
-
-template <int spacedim>
-template <typename InputVector, class DH>
-void
-KellyErrorEstimator<1,spacedim>::
-estimate (const DH   &dof_handler,
-          const hp::QCollection<0> &quadrature,
-          const typename FunctionMap<spacedim>::type &neumann_bc,
-          const std::vector<const InputVector *> &solutions,
-          std::vector<Vector<float>*> &errors,
-          const ComponentMask &component_mask,
-          const Function<spacedim>     *coefficients,
-          const unsigned int       n_threads,
-          const types::subdomain_id subdomain_id,
-          const types::material_id       material_id)
-{
-  estimate(StaticMappingQ1<1,spacedim>::mapping, dof_handler, quadrature, neumann_bc, solutions,
-           errors, component_mask, coefficients, n_threads, subdomain_id, material_id);
-}
-
-
-
-
-template <int spacedim>
-template <typename InputVector, class DH>
-void KellyErrorEstimator<1,spacedim>::
-estimate (const Mapping<1,spacedim>                    &/*mapping*/,
-          const DH                            &/*dof_handler*/,
-          const hp::QCollection<0> &,
-          const typename FunctionMap<spacedim>::type          &/*neumann_bc*/,
-          const std::vector<const InputVector *> &/*solutions*/,
-          std::vector<Vector<float>*>            &/*errors*/,
-          const ComponentMask                &/*component_mask_*/,
-          const Function<spacedim>                   */*coefficient*/,
-          const unsigned int,
-          const types::subdomain_id          /*subdomain_id*/,
-          const types::material_id                   /*material_id*/)
-{
-  Assert (false, ExcInternalError());
-}
-
-
-
-template <int spacedim>
-template <typename InputVector, class DH>
-void KellyErrorEstimator<1,spacedim>::
-estimate (const Mapping<1,spacedim>                    &mapping,
-          const DH                 &dof_handler,
-          const Quadrature<0> &,
-          const typename FunctionMap<spacedim>::type          &neumann_bc,
-          const std::vector<const InputVector *> &solutions,
-          std::vector<Vector<float>*>              &errors,
-          const ComponentMask                  &component_mask,
-          const Function<spacedim>                   *coefficient,
-          const unsigned int,
-          const types::subdomain_id         subdomain_id_,
-          const types::material_id                  material_id)
-{
-#ifdef DEAL_II_WITH_P4EST
-  if (dynamic_cast<const parallel::distributed::Triangulation<1,spacedim>*>
-      (&dof_handler.get_tria())
-      != 0)
-    Assert ((subdomain_id_ == numbers::invalid_subdomain_id)
-            ||
-            (subdomain_id_ ==
-             dynamic_cast<const parallel::distributed::Triangulation<1,spacedim>&>
-             (dof_handler.get_tria()).locally_owned_subdomain()),
-            ExcMessage ("For parallel distributed triangulations, the only "
-                        "valid subdomain_id that can be passed here is the "
-                        "one that corresponds to the locally owned subdomain id."));
-
-  const types::subdomain_id subdomain_id
-    = ((dynamic_cast<const parallel::distributed::Triangulation<1,spacedim>*>
-        (&dof_handler.get_tria())
-        != 0)
-       ?
-       dynamic_cast<const parallel::distributed::Triangulation<1,spacedim>&>
-       (dof_handler.get_tria()).locally_owned_subdomain()
-       :
-       subdomain_id_);
-#else
-  const types::subdomain_id subdomain_id
-    = subdomain_id_;
-#endif
-
-  const unsigned int n_components       = dof_handler.get_fe().n_components();
-  const unsigned int n_solution_vectors = solutions.size();
-
-  // sanity checks
-  Assert (neumann_bc.find(numbers::internal_face_boundary_id) == neumann_bc.end(),
-          ExcInvalidBoundaryIndicator());
-
-  for (typename FunctionMap<spacedim>::type::const_iterator i=neumann_bc.begin();
-       i!=neumann_bc.end(); ++i)
-    Assert (i->second->n_components == n_components, ExcInvalidBoundaryFunction());
-
-  Assert (component_mask.represents_n_components(n_components),
-          ExcInvalidComponentMask());
-  Assert (component_mask.n_selected_components(n_components) > 0,
-          ExcInvalidComponentMask());
-
-  Assert ((coefficient == 0) ||
-          (coefficient->n_components == n_components) ||
-          (coefficient->n_components == 1),
-          ExcInvalidCoefficient());
-
-  Assert (solutions.size() > 0,
-          ExcNoSolutions());
-  Assert (solutions.size() == errors.size(),
-          ExcIncompatibleNumberOfElements(solutions.size(), errors.size()));
-  for (unsigned int n=0; n<solutions.size(); ++n)
-    Assert (solutions[n]->size() == dof_handler.n_dofs(),
-            ExcInvalidSolutionVector());
-
-  Assert ((coefficient == 0) ||
-          (coefficient->n_components == n_components) ||
-          (coefficient->n_components == 1),
-          ExcInvalidCoefficient());
-
-  for (typename FunctionMap<spacedim>::type::const_iterator i=neumann_bc.begin();
-       i!=neumann_bc.end(); ++i)
-    Assert (i->second->n_components == n_components,
-            ExcInvalidBoundaryFunction());
-
-  // reserve one slot for each cell and set
-  // it to zero
-  for (unsigned int n=0; n<n_solution_vectors; ++n)
-    (*errors[n]).reinit (dof_handler.get_tria().n_active_cells());
-
-  // fields to get the gradients on
-  // the present and the neighbor cell.
-  //
-  // for the neighbor gradient, we
-  // need several auxiliary fields,
-  // depending on the way we get it
-  // (see below)
-  std::vector<std::vector<std::vector<Tensor<1,spacedim> > > >
-  gradients_here (n_solution_vectors,
-                  std::vector<std::vector<Tensor<1,spacedim> > >(2, std::vector<Tensor<1,spacedim> >(n_components)));
-  std::vector<std::vector<std::vector<Tensor<1,spacedim> > > >
-  gradients_neighbor (gradients_here);
-  std::vector<Vector<double> >
-  grad_neighbor (n_solution_vectors, Vector<double>(n_components));
-
-  // reserve some space for
-  // coefficient values at one point.
-  // if there is no coefficient, then
-  // we fill it by unity once and for
-  // all and don't set it any more
-  Vector<double> coefficient_values (n_components);
-  if (coefficient == 0)
-    for (unsigned int c=0; c<n_components; ++c)
-      coefficient_values(c) = 1;
-
-  const QTrapez<1> quadrature;
-  const hp::QCollection<1> q_collection(quadrature);
-
-  const hp::FECollection<1,spacedim> fe (dof_handler.get_fe());
-
-  hp::MappingCollection<1,spacedim> mapping_collection;
-  mapping_collection.push_back (mapping);
-
-  hp::FEValues<1,spacedim> fe_values (mapping_collection, fe, q_collection,
-                                      update_gradients);
-
-  // loop over all cells and do something on
-  // the cells which we're told to work
-  // on. note that the error indicator is
-  // only a sum over the two contributions
-  // from the two vertices of each cell.
-  typename DH::active_cell_iterator cell = dof_handler.begin_active();
-  for (unsigned int cell_index=0; cell != dof_handler.end();
-       ++cell, ++cell_index)
-    if (((subdomain_id == numbers::invalid_subdomain_id)
-         ||
-         (cell->subdomain_id() == subdomain_id))
-        &&
-        ((material_id == numbers::invalid_material_id)
-         ||
-         (cell->material_id() == material_id)))
-      {
-        for (unsigned int n=0; n<n_solution_vectors; ++n)
-          (*errors[n])(cell_index) = 0;
-
-        // loop over the two points bounding
-        // this line. n==0 is left point,
-        // n==1 is right point
-        for (unsigned int n=0; n<2; ++n)
-          {
-            // find left or right active
-            // neighbor
-            typename DH::cell_iterator neighbor = cell->neighbor(n);
-            if (neighbor.state() == IteratorState::valid)
-              while (neighbor->has_children())
-                neighbor = neighbor->child(n==0 ? 1 : 0);
-
-            // now get the gradients on the
-            // both sides of the point
-            fe_values.reinit (cell);
-
-            for (unsigned int s=0; s<n_solution_vectors; ++s)
-              fe_values.get_present_fe_values()
-              .get_function_gradients (*solutions[s], gradients_here[s]);
-
-            if (neighbor.state() == IteratorState::valid)
-              {
-                fe_values.reinit (neighbor);
-
-                for (unsigned int s=0; s<n_solution_vectors; ++s)
-                  fe_values.get_present_fe_values()
-                  .get_function_gradients (*solutions[s],
-                                           gradients_neighbor[s]);
-
-                // extract the
-                // gradients of all the
-                // components. [0]
-                // means: x-derivative,
-                // which is the only
-                // one here
-                for (unsigned int s=0; s<n_solution_vectors; ++s)
-                  for (unsigned int c=0; c<n_components; ++c)
-                    grad_neighbor[s](c)
-                      = gradients_neighbor[s][n==0 ? 1 : 0][c][0];
-              }
-            else if (neumann_bc.find(n) != neumann_bc.end())
-              // if Neumann b.c., then fill
-              // the gradients field which
-              // will be used later on.
-              {
-                if (n_components==1)
-                  {
-                    const double
-                    v = neumann_bc.find(n)->second->value(cell->vertex(0));
-
-                    for (unsigned int s=0; s<n_solution_vectors; ++s)
-                      grad_neighbor[s](0) = v;
-                  }
-                else
-                  {
-                    Vector<double> v(n_components);
-                    neumann_bc.find(n)->second->vector_value(cell->vertex(0), v);
-
-                    for (unsigned int s=0; s<n_solution_vectors; ++s)
-                      grad_neighbor[s] = v;
-                  }
-              }
-            else
-              // fill with zeroes.
-              for (unsigned int s=0; s<n_solution_vectors; ++s)
-                grad_neighbor[s] = 0;
-
-            // if there is a
-            // coefficient, then
-            // evaluate it at the
-            // present position. if
-            // there is none, reuse the
-            // preset values.
-            if (coefficient != 0)
-              {
-                if (coefficient->n_components == 1)
-                  {
-                    const double c_value = coefficient->value (cell->vertex(n));
-                    for (unsigned int c=0; c<n_components; ++c)
-                      coefficient_values(c) = c_value;
-                  }
-                else
-                  coefficient->vector_value(cell->vertex(n),
-                                            coefficient_values);
-              }
-
-
-            for (unsigned int s=0; s<n_solution_vectors; ++s)
-              for (unsigned int component=0; component<n_components; ++component)
-                if (component_mask[component] == true)
-                  {
-                    // get gradient
-                    // here. [0] means
-                    // x-derivative
-                    // (there is no
-                    // other component
-                    // in 1d)
-                    const double grad_here = gradients_here[s][n][component][0];
-
-                    const double jump = ((grad_here - grad_neighbor[s](component)) *
-                                         coefficient_values(component));
-                    (*errors[s])(cell_index) += jump*jump * cell->diameter();
-                  }
-          }
-
-        for (unsigned int s=0; s<n_solution_vectors; ++s)
-          (*errors[s])(cell_index) = std::sqrt((*errors[s])(cell_index));
-      }
-}
-
-
-
-
-
-
 // the following function is still independent of dimension, but it
 // calls dimension dependent functions
 template <int dim, int spacedim>
@@ -1581,20 +960,14 @@ estimate (const Mapping<dim, spacedim>                  &mapping,
 
   const unsigned int n_solution_vectors = solutions.size();
 
-  // Map of integrals indexed by
-  // the corresponding face. In this map
-  // we store the integrated jump of the
-  // gradient for each face.
-  // At the end of the function, we again
-  // loop over the cells and collect the
-  // contributions of the different faces
-  // of the cell.
+  // Map of integrals indexed by the corresponding face. In this map we store
+  // the integrated jump of the gradient for each face.  At the end of the
+  // function, we again loop over the cells and collect the contributions of
+  // the different faces of the cell.
   std::map<typename DH::face_iterator,std::vector<double> > face_integrals;
 
-  // all the data needed in the error
-  // estimator by each of the threads
-  // is gathered in the following
-  // stuctures
+  // all the data needed in the error estimator by each of the threads is
+  // gathered in the following structures
   const hp::MappingCollection<dim,spacedim> mapping_collection(mapping);
   const internal::ParallelData<DH>
   parallel_data (dof_handler.get_fe(),
@@ -1612,18 +985,16 @@ estimate (const Mapping<dim, spacedim>                  &mapping,
   // now let's work on all those cells:
   WorkStream::run (dof_handler.begin_active(),
                    static_cast<typename DH::active_cell_iterator>(dof_handler.end()),
-                   std_cxx1x::bind (&internal::estimate_one_cell<dim,spacedim,InputVector,DH>,
+                   std_cxx1x::bind (&internal::estimate_one_cell<InputVector,DH>,
                                     std_cxx1x::_1, std_cxx1x::_2, std_cxx1x::_3, std_cxx1x::ref(solutions)),
                    std_cxx1x::bind (&internal::copy_local_to_global<DH>,
                                     std_cxx1x::_1, std_cxx1x::ref(face_integrals)),
                    parallel_data,
                    sample_local_face_integrals);
 
-  // finally add up the contributions of the
-  // faces for each cell
+  // finally add up the contributions of the faces for each cell
 
-  // reserve one slot for each cell and set
-  // it to zero
+  // reserve one slot for each cell and set it to zero
   for (unsigned int n=0; n<n_solution_vectors; ++n)
     {
       (*errors[n]).reinit (dof_handler.get_tria().n_active_cells());
@@ -1631,10 +1002,8 @@ estimate (const Mapping<dim, spacedim>                  &mapping,
         (*errors[n])(i)=0;
     }
 
-  // now walk over all cells and collect
-  // information from the faces. only do
-  // something if this is a cell we care for
-  // based on the subdomain id
+  // now walk over all cells and collect information from the faces. only do
+  // something if this is a cell we care for based on the subdomain id
   unsigned int present_cell=0;
   for (typename DH::active_cell_iterator cell=dof_handler.begin_active();
        cell!=dof_handler.end();
@@ -1658,9 +1027,8 @@ estimate (const Mapping<dim, spacedim>                  &mapping,
 
             for (unsigned int n=0; n<n_solution_vectors; ++n)
               {
-                // make sure that we have
-                // written a meaningful value
-                // into this slot
+                // make sure that we have written a meaningful value into this
+                // slot
                 Assert (face_integrals[cell->face(face_no)][n] >= 0,
                         ExcInternalError());
 

diff --git a/deal.II/source/numerics/error_estimator.inst.in b/deal.II/source/numerics/error_estimator.inst.in
index 34b06e672b499cdd56cd51d6af0c79fa844a1cc7..8ca601fbcfdb39c5babe7ba6d3b089f3df78bd39 100644 (file)
--- a/deal.II/source/numerics/error_estimator.inst.in
+++ b/deal.II/source/numerics/error_estimator.inst.in
@@ -24,7 +24,7 @@ template class KellyErrorEstimator<deal_II_dimension, deal_II_space_dimension>;
 
 for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS; DH : DOFHANDLER_TEMPLATES )
 {
-#if deal_II_dimension <= deal_II_space_dimension
+#if deal_II_dimension != 1 && deal_II_dimension <= deal_II_space_dimension
 
 template
 void

diff --git a/deal.II/source/numerics/error_estimator_1d.cc b/deal.II/source/numerics/error_estimator_1d.cc
new file mode 100644 (file)
index 0000000..15b65d0
--- /dev/null
+++ b/deal.II/source/numerics/error_estimator_1d.cc
@@ -0,0 +1,439 @@
+// ---------------------------------------------------------------------
+// $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.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/thread_management.h>
+#include <deal.II/base/quadrature.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/work_stream.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/parallel_vector.h>
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/parallel_block_vector.h>
+#include <deal.II/lac/petsc_vector.h>
+#include <deal.II/lac/petsc_block_vector.h>
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/lac/trilinos_block_vector.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/base/geometry_info.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/fe/fe.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/hp/fe_values.h>
+#include <deal.II/fe/fe_update_flags.h>
+#include <deal.II/fe/mapping_q1.h>
+#include <deal.II/hp/q_collection.h>
+#include <deal.II/hp/mapping_collection.h>
+#include <deal.II/numerics/error_estimator.h>
+
+#include <deal.II/base/std_cxx1x/bind.h>
+
+#include <numeric>
+#include <algorithm>
+#include <cmath>
+#include <vector>
+
+DEAL_II_NAMESPACE_OPEN
+
+
+
+template <int spacedim>
+template <typename InputVector, class DH>
+void
+KellyErrorEstimator<1,spacedim>::
+estimate (const Mapping<1,spacedim>      &mapping,
+          const DH   &dof_handler,
+          const Quadrature<0> &quadrature,
+          const typename FunctionMap<spacedim>::type &neumann_bc,
+          const InputVector       &solution,
+          Vector<float>           &error,
+          const ComponentMask &component_mask,
+          const Function<spacedim>     *coefficients,
+          const unsigned int       n_threads,
+          const types::subdomain_id subdomain_id,
+          const types::material_id       material_id)
+{
+  // just pass on to the other function
+  const std::vector<const InputVector *> solutions (1, &solution);
+  std::vector<Vector<float>*>              errors (1, &error);
+  estimate (mapping, dof_handler, quadrature, neumann_bc, solutions, errors,
+            component_mask, coefficients, n_threads, subdomain_id, material_id);
+}
+
+
+
+template <int spacedim>
+template <typename InputVector, class DH>
+void
+KellyErrorEstimator<1,spacedim>::
+estimate (const DH   &dof_handler,
+          const Quadrature<0> &quadrature,
+          const typename FunctionMap<spacedim>::type &neumann_bc,
+          const InputVector       &solution,
+          Vector<float>           &error,
+          const ComponentMask &component_mask,
+          const Function<spacedim>     *coefficients,
+          const unsigned int       n_threads,
+          const types::subdomain_id subdomain_id,
+          const types::material_id       material_id)
+{
+  estimate(StaticMappingQ1<1,spacedim>::mapping, dof_handler, quadrature, neumann_bc, solution,
+           error, component_mask, coefficients, n_threads, subdomain_id, material_id);
+}
+
+
+
+template <int spacedim>
+template <typename InputVector, class DH>
+void
+KellyErrorEstimator<1,spacedim>::
+estimate (const DH   &dof_handler,
+          const Quadrature<0> &quadrature,
+          const typename FunctionMap<spacedim>::type &neumann_bc,
+          const std::vector<const InputVector *> &solutions,
+          std::vector<Vector<float>*> &errors,
+          const ComponentMask &component_mask,
+          const Function<spacedim>     *coefficients,
+          const unsigned int       n_threads,
+          const types::subdomain_id subdomain_id,
+          const types::material_id       material_id)
+{
+  estimate(StaticMappingQ1<1,spacedim>::mapping, dof_handler, quadrature, neumann_bc, solutions,
+           errors, component_mask, coefficients, n_threads, subdomain_id, material_id);
+}
+
+
+
+template <int spacedim>
+template <typename InputVector, class DH>
+void
+KellyErrorEstimator<1,spacedim>::
+estimate (const Mapping<1,spacedim>      &mapping,
+          const DH   &dof_handler,
+          const hp::QCollection<0> &quadrature,
+          const typename FunctionMap<spacedim>::type &neumann_bc,
+          const InputVector       &solution,
+          Vector<float>           &error,
+          const ComponentMask &component_mask,
+          const Function<spacedim>     *coefficients,
+          const unsigned int       n_threads,
+          const types::subdomain_id subdomain_id,
+          const types::material_id       material_id)
+{
+  // just pass on to the other function
+  const std::vector<const InputVector *> solutions (1, &solution);
+  std::vector<Vector<float>*>              errors (1, &error);
+  estimate (mapping, dof_handler, quadrature, neumann_bc, solutions, errors,
+            component_mask, coefficients, n_threads, subdomain_id, material_id);
+}
+
+
+template <int spacedim>
+template <typename InputVector, class DH>
+void
+KellyErrorEstimator<1,spacedim>::
+estimate (const DH   &dof_handler,
+          const hp::QCollection<0> &quadrature,
+          const typename FunctionMap<spacedim>::type &neumann_bc,
+          const InputVector       &solution,
+          Vector<float>           &error,
+          const ComponentMask &component_mask,
+          const Function<spacedim>     *coefficients,
+          const unsigned int       n_threads,
+          const types::subdomain_id subdomain_id,
+          const types::material_id       material_id)
+{
+  estimate(StaticMappingQ1<1,spacedim>::mapping, dof_handler, quadrature, neumann_bc, solution,
+           error, component_mask, coefficients, n_threads, subdomain_id, material_id);
+}
+
+
+
+template <int spacedim>
+template <typename InputVector, class DH>
+void
+KellyErrorEstimator<1,spacedim>::
+estimate (const DH   &dof_handler,
+          const hp::QCollection<0> &quadrature,
+          const typename FunctionMap<spacedim>::type &neumann_bc,
+          const std::vector<const InputVector *> &solutions,
+          std::vector<Vector<float>*> &errors,
+          const ComponentMask &component_mask,
+          const Function<spacedim>     *coefficients,
+          const unsigned int       n_threads,
+          const types::subdomain_id subdomain_id,
+          const types::material_id       material_id)
+{
+  estimate(StaticMappingQ1<1,spacedim>::mapping, dof_handler, quadrature, neumann_bc, solutions,
+           errors, component_mask, coefficients, n_threads, subdomain_id, material_id);
+}
+
+
+
+
+template <int spacedim>
+template <typename InputVector, class DH>
+void KellyErrorEstimator<1,spacedim>::
+estimate (const Mapping<1,spacedim>                    &/*mapping*/,
+          const DH                            &/*dof_handler*/,
+          const hp::QCollection<0> &,
+          const typename FunctionMap<spacedim>::type          &/*neumann_bc*/,
+          const std::vector<const InputVector *> &/*solutions*/,
+          std::vector<Vector<float>*>            &/*errors*/,
+          const ComponentMask                &/*component_mask_*/,
+          const Function<spacedim>                   */*coefficient*/,
+          const unsigned int,
+          const types::subdomain_id          /*subdomain_id*/,
+          const types::material_id                   /*material_id*/)
+{
+  Assert (false, ExcInternalError());
+}
+
+
+
+template <int spacedim>
+template <typename InputVector, class DH>
+void KellyErrorEstimator<1,spacedim>::
+estimate (const Mapping<1,spacedim>                    &mapping,
+          const DH                 &dof_handler,
+          const Quadrature<0> &,
+          const typename FunctionMap<spacedim>::type          &neumann_bc,
+          const std::vector<const InputVector *> &solutions,
+          std::vector<Vector<float>*>              &errors,
+          const ComponentMask                  &component_mask,
+          const Function<spacedim>                   *coefficient,
+          const unsigned int,
+          const types::subdomain_id         subdomain_id_,
+          const types::material_id                  material_id)
+{
+#ifdef DEAL_II_WITH_P4EST
+  if (dynamic_cast<const parallel::distributed::Triangulation<1,spacedim>*>
+      (&dof_handler.get_tria())
+      != 0)
+    Assert ((subdomain_id_ == numbers::invalid_subdomain_id)
+            ||
+            (subdomain_id_ ==
+             dynamic_cast<const parallel::distributed::Triangulation<1,spacedim>&>
+             (dof_handler.get_tria()).locally_owned_subdomain()),
+            ExcMessage ("For parallel distributed triangulations, the only "
+                        "valid subdomain_id that can be passed here is the "
+                        "one that corresponds to the locally owned subdomain id."));
+
+  const types::subdomain_id subdomain_id
+    = ((dynamic_cast<const parallel::distributed::Triangulation<1,spacedim>*>
+        (&dof_handler.get_tria())
+        != 0)
+       ?
+       dynamic_cast<const parallel::distributed::Triangulation<1,spacedim>&>
+       (dof_handler.get_tria()).locally_owned_subdomain()
+       :
+       subdomain_id_);
+#else
+  const types::subdomain_id subdomain_id
+    = subdomain_id_;
+#endif
+
+  const unsigned int n_components       = dof_handler.get_fe().n_components();
+  const unsigned int n_solution_vectors = solutions.size();
+
+  // sanity checks
+  Assert (neumann_bc.find(numbers::internal_face_boundary_id) == neumann_bc.end(),
+          ExcInvalidBoundaryIndicator());
+
+  for (typename FunctionMap<spacedim>::type::const_iterator i=neumann_bc.begin();
+       i!=neumann_bc.end(); ++i)
+    Assert (i->second->n_components == n_components, ExcInvalidBoundaryFunction());
+
+  Assert (component_mask.represents_n_components(n_components),
+          ExcInvalidComponentMask());
+  Assert (component_mask.n_selected_components(n_components) > 0,
+          ExcInvalidComponentMask());
+
+  Assert ((coefficient == 0) ||
+          (coefficient->n_components == n_components) ||
+          (coefficient->n_components == 1),
+          ExcInvalidCoefficient());
+
+  Assert (solutions.size() > 0,
+          ExcNoSolutions());
+  Assert (solutions.size() == errors.size(),
+          ExcIncompatibleNumberOfElements(solutions.size(), errors.size()));
+  for (unsigned int n=0; n<solutions.size(); ++n)
+    Assert (solutions[n]->size() == dof_handler.n_dofs(),
+            ExcInvalidSolutionVector());
+
+  Assert ((coefficient == 0) ||
+          (coefficient->n_components == n_components) ||
+          (coefficient->n_components == 1),
+          ExcInvalidCoefficient());
+
+  for (typename FunctionMap<spacedim>::type::const_iterator i=neumann_bc.begin();
+       i!=neumann_bc.end(); ++i)
+    Assert (i->second->n_components == n_components,
+            ExcInvalidBoundaryFunction());
+
+  // reserve one slot for each cell and set it to zero
+  for (unsigned int n=0; n<n_solution_vectors; ++n)
+    (*errors[n]).reinit (dof_handler.get_tria().n_active_cells());
+
+  // fields to get the gradients on the present and the neighbor cell.
+  //
+  // for the neighbor gradient, we need several auxiliary fields, depending on
+  // the way we get it (see below)
+  std::vector<std::vector<std::vector<Tensor<1,spacedim> > > >
+  gradients_here (n_solution_vectors,
+                  std::vector<std::vector<Tensor<1,spacedim> > >(2, std::vector<Tensor<1,spacedim> >(n_components)));
+  std::vector<std::vector<std::vector<Tensor<1,spacedim> > > >
+  gradients_neighbor (gradients_here);
+  std::vector<Vector<double> >
+  grad_neighbor (n_solution_vectors, Vector<double>(n_components));
+
+  // reserve some space for coefficient values at one point.  if there is no
+  // coefficient, then we fill it by unity once and for all and don't set it
+  // any more
+  Vector<double> coefficient_values (n_components);
+  if (coefficient == 0)
+    for (unsigned int c=0; c<n_components; ++c)
+      coefficient_values(c) = 1;
+
+  const QTrapez<1> quadrature;
+  const hp::QCollection<1> q_collection(quadrature);
+
+  const hp::FECollection<1,spacedim> fe (dof_handler.get_fe());
+
+  hp::MappingCollection<1,spacedim> mapping_collection;
+  mapping_collection.push_back (mapping);
+
+  hp::FEValues<1,spacedim> fe_values (mapping_collection, fe, q_collection,
+                                      update_gradients);
+
+  // loop over all cells and do something on the cells which we're told to
+  // work on. note that the error indicator is only a sum over the two
+  // contributions from the two vertices of each cell.
+  typename DH::active_cell_iterator cell = dof_handler.begin_active();
+  for (unsigned int cell_index=0; cell != dof_handler.end();
+       ++cell, ++cell_index)
+    if (((subdomain_id == numbers::invalid_subdomain_id)
+         ||
+         (cell->subdomain_id() == subdomain_id))
+        &&
+        ((material_id == numbers::invalid_material_id)
+         ||
+         (cell->material_id() == material_id)))
+      {
+        for (unsigned int n=0; n<n_solution_vectors; ++n)
+          (*errors[n])(cell_index) = 0;
+
+        // loop over the two points bounding this line. n==0 is left point,
+        // n==1 is right point
+        for (unsigned int n=0; n<2; ++n)
+          {
+            // find left or right active neighbor
+            typename DH::cell_iterator neighbor = cell->neighbor(n);
+            if (neighbor.state() == IteratorState::valid)
+              while (neighbor->has_children())
+                neighbor = neighbor->child(n==0 ? 1 : 0);
+
+            // now get the gradients on the both sides of the point
+            fe_values.reinit (cell);
+
+            for (unsigned int s=0; s<n_solution_vectors; ++s)
+              fe_values.get_present_fe_values()
+              .get_function_gradients (*solutions[s], gradients_here[s]);
+
+            if (neighbor.state() == IteratorState::valid)
+              {
+                fe_values.reinit (neighbor);
+
+                for (unsigned int s=0; s<n_solution_vectors; ++s)
+                  fe_values.get_present_fe_values()
+                  .get_function_gradients (*solutions[s],
+                                           gradients_neighbor[s]);
+
+                // extract the gradients of all the components. [0] means:
+                // x-derivative, which is the only one here
+                for (unsigned int s=0; s<n_solution_vectors; ++s)
+                  for (unsigned int c=0; c<n_components; ++c)
+                    grad_neighbor[s](c)
+                      = gradients_neighbor[s][n==0 ? 1 : 0][c][0];
+              }
+            else if (neumann_bc.find(n) != neumann_bc.end())
+              // if Neumann b.c., then fill the gradients field which will be
+              // used later on.
+              {
+                if (n_components==1)
+                  {
+                    const double
+                    v = neumann_bc.find(n)->second->value(cell->vertex(0));
+
+                    for (unsigned int s=0; s<n_solution_vectors; ++s)
+                      grad_neighbor[s](0) = v;
+                  }
+                else
+                  {
+                    Vector<double> v(n_components);
+                    neumann_bc.find(n)->second->vector_value(cell->vertex(0), v);
+
+                    for (unsigned int s=0; s<n_solution_vectors; ++s)
+                      grad_neighbor[s] = v;
+                  }
+              }
+            else
+              // fill with zeroes.
+              for (unsigned int s=0; s<n_solution_vectors; ++s)
+                grad_neighbor[s] = 0;
+
+            // if there is a coefficient, then evaluate it at the present
+            // position. if there is none, reuse the preset values.
+            if (coefficient != 0)
+              {
+                if (coefficient->n_components == 1)
+                  {
+                    const double c_value = coefficient->value (cell->vertex(n));
+                    for (unsigned int c=0; c<n_components; ++c)
+                      coefficient_values(c) = c_value;
+                  }
+                else
+                  coefficient->vector_value(cell->vertex(n),
+                                            coefficient_values);
+              }
+
+
+            for (unsigned int s=0; s<n_solution_vectors; ++s)
+              for (unsigned int component=0; component<n_components; ++component)
+                if (component_mask[component] == true)
+                  {
+                    // get gradient here. [0] means x-derivative (there is no
+                    // other component in 1d)
+                    const double grad_here = gradients_here[s][n][component][0];
+
+                    const double jump = ((grad_here - grad_neighbor[s](component)) *
+                                         coefficient_values(component));
+                    (*errors[s])(cell_index) += jump*jump * cell->diameter();
+                  }
+          }
+
+        for (unsigned int s=0; s<n_solution_vectors; ++s)
+          (*errors[s])(cell_index) = std::sqrt((*errors[s])(cell_index));
+      }
+}
+
+
+// explicit instantiations
+#include "error_estimator_1d.inst"
+
+
+DEAL_II_NAMESPACE_CLOSE

diff --git a/deal.II/source/numerics/error_estimator_1d.inst.in b/deal.II/source/numerics/error_estimator_1d.inst.in
new file mode 100644 (file)
index 0000000..66b0c78
--- /dev/null
+++ b/deal.II/source/numerics/error_estimator_1d.inst.in
@@ -0,0 +1,143 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2010 - 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 (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS; DH : DOFHANDLER_TEMPLATES )
+{
+#if deal_II_dimension == 1 && deal_II_dimension <= deal_II_space_dimension
+
+template
+void
+KellyErrorEstimator<deal_II_dimension, deal_II_space_dimension>::
+estimate<VEC,DH<deal_II_dimension,deal_II_space_dimension> > (const Mapping<deal_II_dimension, deal_II_space_dimension>      &,
+          const DH<deal_II_dimension,deal_II_space_dimension>   &,
+          const Quadrature<deal_II_dimension-1> &,
+          const FunctionMap<deal_II_space_dimension>::type &,
+          const VEC       &,
+          Vector<float>           &,
+          const ComponentMask &,
+          const Function<deal_II_space_dimension>     *,
+          const unsigned int       ,
+          const unsigned int       ,
+          const types::material_id);
+
+template
+void
+KellyErrorEstimator<deal_II_dimension, deal_II_space_dimension>::
+estimate<VEC,DH<deal_II_dimension,deal_II_space_dimension> > (
+          const DH<deal_II_dimension,deal_II_space_dimension>   &,
+          const Quadrature<deal_II_dimension-1> &,
+          const FunctionMap<deal_II_space_dimension>::type &,
+          const VEC       &,
+          Vector<float>           &,
+          const ComponentMask &,
+          const Function<deal_II_space_dimension>     *,
+          const unsigned int       ,
+          const unsigned int       ,
+          const types::material_id);
+
+template
+void
+KellyErrorEstimator<deal_II_dimension, deal_II_space_dimension>::
+estimate<VEC,DH<deal_II_dimension,deal_II_space_dimension> > (const Mapping<deal_II_dimension, deal_II_space_dimension>      &,
+          const DH<deal_II_dimension,deal_II_space_dimension>   &,
+          const Quadrature<deal_II_dimension-1> &,
+          const FunctionMap<deal_II_space_dimension>::type &,
+          const std::vector<const VEC *>       &,
+          std::vector<Vector<float> *>         &,
+          const ComponentMask &,
+          const Function<deal_II_space_dimension>     *,
+          const unsigned int       ,
+          const unsigned int       ,
+          const types::material_id);
+
+template
+void
+KellyErrorEstimator<deal_II_dimension, deal_II_space_dimension>::
+estimate<VEC,DH<deal_II_dimension,deal_II_space_dimension> > (
+          const DH<deal_II_dimension,deal_II_space_dimension>   &,
+          const Quadrature<deal_II_dimension-1> &,
+          const FunctionMap<deal_II_space_dimension>::type &,
+          const std::vector<const VEC *>       &,
+          std::vector<Vector<float> *>         &,
+          const ComponentMask &,
+          const Function<deal_II_space_dimension>     *,
+          const unsigned int       ,
+          const unsigned int       ,
+          const types::material_id);
+
+template
+void
+KellyErrorEstimator<deal_II_dimension, deal_II_space_dimension>::
+estimate<VEC,DH<deal_II_dimension,deal_II_space_dimension> > (const Mapping<deal_II_dimension, deal_II_space_dimension>      &,
+          const DH<deal_II_dimension,deal_II_space_dimension>   &,
+          const hp::QCollection<deal_II_dimension-1> &,
+          const FunctionMap<deal_II_space_dimension>::type &,
+          const VEC       &,
+          Vector<float>           &,
+          const ComponentMask &,
+          const Function<deal_II_space_dimension>     *,
+          const unsigned int       ,
+          const unsigned int       ,
+          const types::material_id);
+
+template
+void
+KellyErrorEstimator<deal_II_dimension, deal_II_space_dimension>::
+estimate<VEC,DH<deal_II_dimension,deal_II_space_dimension> > (
+          const DH<deal_II_dimension,deal_II_space_dimension>   &,
+          const hp::QCollection<deal_II_dimension-1> &,
+          const FunctionMap<deal_II_space_dimension>::type &,
+          const VEC       &,
+          Vector<float>           &,
+          const ComponentMask &,
+          const Function<deal_II_space_dimension>     *,
+          const unsigned int       ,
+          const unsigned int       ,
+          const types::material_id);
+
+template
+void
+KellyErrorEstimator<deal_II_dimension, deal_II_space_dimension>::
+estimate<VEC,DH<deal_II_dimension,deal_II_space_dimension> > (const Mapping<deal_II_dimension, deal_II_space_dimension>      &,
+          const DH<deal_II_dimension,deal_II_space_dimension>   &,
+          const hp::QCollection<deal_II_dimension-1> &,
+          const FunctionMap<deal_II_space_dimension>::type &,
+          const std::vector<const VEC *>       &,
+          std::vector<Vector<float> *>         &,
+          const ComponentMask &,
+          const Function<deal_II_space_dimension>     *,
+          const unsigned int       ,
+          const unsigned int       ,
+          const types::material_id);
+
+template
+void
+KellyErrorEstimator<deal_II_dimension, deal_II_space_dimension>::
+estimate<VEC,DH<deal_II_dimension,deal_II_space_dimension> > (
+          const DH<deal_II_dimension,deal_II_space_dimension>   &,
+          const hp::QCollection<deal_II_dimension-1> &,
+          const FunctionMap<deal_II_space_dimension>::type &,
+          const std::vector<const VEC *>       &,
+          std::vector<Vector<float> *>         &,
+          const ComponentMask &,
+          const Function<deal_II_space_dimension>     *,
+          const unsigned int       ,
+          const unsigned int       ,
+          const types::material_id);
+
+#endif
+}