]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Allow to also select which cells to work on by using the material id, not only the...
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 4 Jan 2006 15:09:17 +0000 (15:09 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 4 Jan 2006 15:09:17 +0000 (15:09 +0000)
git-svn-id: https://svn.dealii.org/trunk@11947 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/numerics/error_estimator.h
deal.II/deal.II/source/numerics/error_estimator.cc
deal.II/doc/news/changes.html
tests/bits/error_estimator_02.cc [new file with mode: 0644]

index 5ca9f0be907f433363f644192d63f71f6540a9ba..2c12d93ab21d2057733f4429d22060fea95c46c6 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005 by the deal.II authors
+//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -290,6 +290,13 @@ class KellyErrorEstimator
                                      * 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, it means that
+                                     * indicators will only be computed for
+                                     * cells with this particular material
+                                     * id.
                                      */
     template <typename InputVector>
     static void estimate (const Mapping<dim>      &mapping,
@@ -301,8 +308,9 @@ class KellyErrorEstimator
                          const std::vector<bool> &component_mask = std::vector<bool>(),
                          const Function<dim>     *coefficients   = 0,
                          const unsigned int       n_threads = multithread_info.n_default_threads,
-                          const unsigned int       subdomain_id = deal_II_numbers::invalid_unsigned_int);
-
+                          const unsigned int       subdomain_id = deal_II_numbers::invalid_unsigned_int,
+                          const unsigned int       material_id = deal_II_numbers::invalid_unsigned_int);
+    
                                     /**
                                      * Calls the @p estimate
                                      * function, see above, with
@@ -317,7 +325,8 @@ class KellyErrorEstimator
                          const std::vector<bool> &component_mask = std::vector<bool>(),
                          const Function<dim>     *coefficients   = 0,
                          const unsigned int       n_threads = multithread_info.n_default_threads,
-                          const unsigned int       subdomain_id = deal_II_numbers::invalid_unsigned_int);
+                          const unsigned int       subdomain_id = deal_II_numbers::invalid_unsigned_int,
+                          const unsigned int       material_id = deal_II_numbers::invalid_unsigned_int);
     
                                     /**
                                      * Same function as above, but
@@ -356,7 +365,8 @@ class KellyErrorEstimator
                          const std::vector<bool>     &component_mask = std::vector<bool>(),
                          const Function<dim>         *coefficients   = 0,
                          const unsigned int           n_threads = multithread_info.n_default_threads,
-                          const unsigned int           subdomain_id = deal_II_numbers::invalid_unsigned_int);
+                          const unsigned int           subdomain_id = deal_II_numbers::invalid_unsigned_int,
+                          const unsigned int           material_id = deal_II_numbers::invalid_unsigned_int);
 
                                     /**
                                      * Calls the @p estimate
@@ -372,7 +382,8 @@ class KellyErrorEstimator
                          const std::vector<bool>     &component_mask = std::vector<bool>(),
                          const Function<dim>         *coefficients   = 0,
                          const unsigned int           n_threads = multithread_info.n_default_threads,
-                          const unsigned int           subdomain_id = deal_II_numbers::invalid_unsigned_int);
+                          const unsigned int           subdomain_id = deal_II_numbers::invalid_unsigned_int,
+                          const unsigned int           material_id = deal_II_numbers::invalid_unsigned_int);
 
     
                                     /**
@@ -543,10 +554,15 @@ class KellyErrorEstimator
        std::vector<double>          JxW_values;
 
                                          /**
-                                          * The subdomain ids we are to care
+                                          * The subdomain id we are to care
                                           * for.
                                           */
         const unsigned int subdomain_id;
+                                         /**
+                                          * The material id we are to care
+                                          * for.
+                                          */
+        const unsigned int material_id;
         
                                         /**
                                          * Constructor.
@@ -554,7 +570,8 @@ class KellyErrorEstimator
        PerThreadData (const unsigned int n_solution_vectors,
                       const unsigned int n_components,
                       const unsigned int n_q_points,
-                       const unsigned int subdomain_id);
+                       const unsigned int subdomain_id,
+                       const unsigned int material_id);
     };
 
 
@@ -699,23 +716,20 @@ class KellyErrorEstimator<1>
                                      * 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.
+                                     * 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
@@ -740,7 +754,8 @@ class KellyErrorEstimator<1>
                          const std::vector<bool> &component_mask = std::vector<bool>(),
                          const Function<1>     *coefficients   = 0,
                          const unsigned int       n_threads = multithread_info.n_default_threads,
-                          const unsigned int       subdomain_id = deal_II_numbers::invalid_unsigned_int);
+                          const unsigned int       subdomain_id = deal_II_numbers::invalid_unsigned_int,
+                          const unsigned int       material_id = deal_II_numbers::invalid_unsigned_int);
 
                                     /**
                                      * Calls the @p estimate
@@ -756,7 +771,8 @@ class KellyErrorEstimator<1>
                          const std::vector<bool> &component_mask = std::vector<bool>(),
                          const Function<1>     *coefficients   = 0,
                          const unsigned int       n_threads = multithread_info.n_default_threads,
-                          const unsigned int       subdomain_id = deal_II_numbers::invalid_unsigned_int);
+                          const unsigned int       subdomain_id = deal_II_numbers::invalid_unsigned_int,
+                          const unsigned int       material_id = deal_II_numbers::invalid_unsigned_int);
     
                                     /**
                                      * Same function as above, but
@@ -795,7 +811,8 @@ class KellyErrorEstimator<1>
                          const std::vector<bool>     &component_mask = std::vector<bool>(),
                          const Function<1>         *coefficients   = 0,
                          const unsigned int           n_threads = multithread_info.n_default_threads,
-                          const unsigned int           subdomain_id = deal_II_numbers::invalid_unsigned_int);
+                          const unsigned int           subdomain_id = deal_II_numbers::invalid_unsigned_int,
+                          const unsigned int           material_id = deal_II_numbers::invalid_unsigned_int);
 
                                     /**
                                      * Calls the @p estimate
@@ -811,7 +828,8 @@ class KellyErrorEstimator<1>
                          const std::vector<bool>     &component_mask = std::vector<bool>(),
                          const Function<1>         *coefficients   = 0,
                          const unsigned int           n_threads = multithread_info.n_default_threads,
-                          const unsigned int           subdomain_id = deal_II_numbers::invalid_unsigned_int);
+                          const unsigned int           subdomain_id = deal_II_numbers::invalid_unsigned_int,
+                          const unsigned int           material_id = deal_II_numbers::invalid_unsigned_int);
 
     
                                     /**
index 05dbb1f59ca5be213930d9fb894e96e897d625df..1b56a8111ca4c52fad389e06204bc449fcd0e0ae 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005 by the deal.II authors
+//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -76,13 +76,14 @@ estimate (const Mapping<1>      &mapping,
           const std::vector<bool> &component_mask,
           const Function<1>     *coefficients,
           const unsigned int       n_threads,
-          const unsigned int       subdomain_id)
+          const unsigned int       subdomain_id,
+          const unsigned int       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);
+           component_mask, coefficients, n_threads, subdomain_id, material_id);
 }
 
 
@@ -97,12 +98,13 @@ estimate (const DoFHandler<1>   &dof_handler,
           const std::vector<bool> &component_mask,
           const Function<1>     *coefficients,
           const unsigned int       n_threads,
-          const unsigned int       subdomain_id)
+          const unsigned int       subdomain_id,
+          const unsigned int       material_id)
 {
   Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
   static const MappingQ1<1> mapping;
   estimate(mapping, dof_handler, quadrature, neumann_bc, solution,
-          error, component_mask, coefficients, n_threads, subdomain_id);
+          error, component_mask, coefficients, n_threads, subdomain_id, material_id);
 }
 
 
@@ -118,12 +120,13 @@ estimate (const DoFHandler<1>   &dof_handler,
           const std::vector<bool> &component_mask,
           const Function<1>     *coefficients,
           const unsigned int       n_threads,
-          const unsigned int       subdomain_id)
+          const unsigned int       subdomain_id,
+          const unsigned int       material_id)
 {
   Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
   static const MappingQ1<1> mapping;
   estimate(mapping, dof_handler, quadrature, neumann_bc, solutions,
-          errors, component_mask, coefficients, n_threads, subdomain_id);
+          errors, component_mask, coefficients, n_threads, subdomain_id, material_id);
 }
 
 
@@ -139,7 +142,8 @@ estimate (const Mapping<1>                    &mapping,
           const std::vector<bool>                  &component_mask_,
           const Function<1>                   *coefficient,
           const unsigned int,
-          const unsigned int                  subdomain_id)
+          const unsigned int                  subdomain_id,
+          const unsigned int                  material_id)
 {
   const unsigned int n_components       = dof_handler.get_fe().n_components();
   const unsigned int n_solution_vectors = solutions.size();
@@ -232,9 +236,13 @@ estimate (const Mapping<1>                    &mapping,
   DoFHandler<1>::active_cell_iterator cell = dof_handler.begin_active();
   for (unsigned int cell_index=0; cell != dof_handler.end();
        ++cell, ++cell_index)
-    if ((subdomain_id == deal_II_numbers::invalid_unsigned_int)
-        ||
-        (cell->subdomain_id() == subdomain_id))
+    if (((subdomain_id == deal_II_numbers::invalid_unsigned_int)
+         ||
+         (cell->subdomain_id() == subdomain_id))
+        &&
+        ((material_id == deal_II_numbers::invalid_unsigned_int)
+         ||
+         (cell->material_id() == material_id)))
       {
         for (unsigned int n=0; n<n_solution_vectors; ++n)
           (*errors[n])(cell_index) = 0;
@@ -357,9 +365,11 @@ KellyErrorEstimator<dim>::PerThreadData::
 PerThreadData (const unsigned int n_solution_vectors,
               const unsigned int n_components,
               const unsigned int n_q_points,
-               const unsigned int subdomain_id)
+               const unsigned int subdomain_id,
+               const unsigned int material_id)
                 :
-                subdomain_id (subdomain_id)
+                subdomain_id (subdomain_id),
+                material_id (material_id)
 {
                                   // Init the size of a lot of vectors
                                   // needed in the calculations once
@@ -409,13 +419,14 @@ estimate (const Mapping<dim>      &mapping,
           const std::vector<bool> &component_mask,
           const Function<dim>     *coefficients,
           const unsigned int       n_threads,
-          const unsigned int       subdomain_id)
+          const unsigned int       subdomain_id,
+          const unsigned int       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);
+           component_mask, coefficients, n_threads, subdomain_id, material_id);
 }
 
 
@@ -431,12 +442,14 @@ estimate (const DoFHandler<dim>   &dof_handler,
           const std::vector<bool> &component_mask,
           const Function<dim>     *coefficients,
           const unsigned int       n_threads,
-          const unsigned int       subdomain_id)
+          const unsigned int       subdomain_id,
+          const unsigned int       material_id)
 {
   Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
   static const MappingQ1<dim> mapping;
   estimate(mapping, dof_handler, quadrature, neumann_bc, solution,
-          error, component_mask, coefficients, n_threads, subdomain_id);
+          error, component_mask, coefficients, n_threads,
+           subdomain_id, material_id);
 }
 
 
@@ -459,6 +472,7 @@ estimate_some (const Mapping<dim>                  &mapping,
   const unsigned int n_solution_vectors = solutions.size();
 
   const unsigned int subdomain_id = per_thread_data.subdomain_id;
+  const unsigned int material_id  = per_thread_data.material_id;
   
                                   // make up a fe face values object
                                   // for the restriction of the
@@ -589,15 +603,20 @@ estimate_some (const Mapping<dim>                  &mapping,
              continue;
            }
 
-                                           // finally: note that we only
-                                           // have to do something if either
-                                           // the present cell is on the
-                                           // subdomain we care for, or if one
-                                           // of the neighbors behind the face
-                                           // is on the subdomain we care for
-          if ( ! ((subdomain_id == deal_II_numbers::invalid_unsigned_int)
-                  ||
-                  (cell->subdomain_id() == subdomain_id)))
+                                           // 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 == deal_II_numbers::invalid_unsigned_int)
+                    ||
+                    (cell->subdomain_id() == subdomain_id))
+                   &&
+                   ((material_id == deal_II_numbers::invalid_unsigned_int)
+                    ||
+                    (cell->material_id() == material_id))) )
             {
                                                // ok, cell is unwanted, but
                                                // maybe its neighbor behind
@@ -609,14 +628,28 @@ estimate_some (const Mapping<dim>                  &mapping,
 
               bool care_for_cell = false;
               if (cell->neighbor(face_no)->has_children() == false)
-                care_for_cell |= (cell->neighbor(face_no)->subdomain_id()
-                                  == subdomain_id);
+                care_for_cell |= ((cell->neighbor(face_no)->subdomain_id()
+                                   == subdomain_id) ||
+                                  (subdomain_id == deal_II_numbers::invalid_unsigned_int))
+                                 &&
+                                 ((cell->neighbor(face_no)->material_id()
+                                   == material_id) ||
+                                  (material_id == deal_II_numbers::invalid_unsigned_int));
               else
                 {
                   for (unsigned int sf=0;
                        sf<GeometryInfo<dim>::subfaces_per_face; ++sf)
-                    if (cell->neighbor_child_on_subface(face_no,sf)
-                        ->subdomain_id() == subdomain_id)
+                    if (((cell->neighbor_child_on_subface(face_no,sf)
+                          ->subdomain_id() == subdomain_id)
+                         &&
+                         (material_id ==
+                          deal_II_numbers::invalid_unsigned_int))
+                        ||
+                        ((cell->neighbor_child_on_subface(face_no,sf)
+                          ->material_id() == material_id)
+                         &&
+                         (subdomain_id ==
+                          deal_II_numbers::invalid_unsigned_int)))
                       {
                         care_for_cell = true;
                         break;
@@ -624,8 +657,9 @@ estimate_some (const Mapping<dim>                  &mapping,
                 }
 
                                                // so if none of the neighbors
-                                               // cares for this subdomain
-                                               // either, then try next face
+                                               // cares for this subdomain or
+                                               // material either, then try
+                                               // next face
               if (care_for_cell == false)
                 continue;
             }
@@ -684,7 +718,8 @@ estimate (const Mapping<dim>                  &mapping,
           const std::vector<bool>                  &component_mask_,
           const Function<dim>                 *coefficients,
           const unsigned int                   n_threads_,
-          const unsigned int                   subdomain_id)
+          const unsigned int                   subdomain_id,
+          const unsigned int                   material_id)
 {
   const unsigned int n_components = dof_handler.get_fe().n_components();
 
@@ -774,7 +809,7 @@ estimate (const Mapping<dim>                  &mapping,
       = new PerThreadData (solutions.size(),
                            dof_handler.get_fe().n_components(),
                            quadrature.n_quadrature_points,
-                           subdomain_id);
+                           subdomain_id, material_id);
   
                                   // split all cells into threads if
                                   // multithreading is used and run
@@ -837,9 +872,13 @@ estimate (const Mapping<dim>                  &mapping,
   for (active_cell_iterator cell=dof_handler.begin_active();
        cell!=dof_handler.end();
        ++cell, ++present_cell)
-    if ((subdomain_id == deal_II_numbers::invalid_unsigned_int)
-        ||
-        (cell->subdomain_id() == subdomain_id))
+    if ( ((subdomain_id == deal_II_numbers::invalid_unsigned_int)
+          ||
+          (cell->subdomain_id() == subdomain_id))
+         &&
+         ((material_id == deal_II_numbers::invalid_unsigned_int)
+          ||
+          (cell->material_id() == material_id)))
       {
                                          // loop over all faces of this cell
         for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell;
@@ -879,12 +918,14 @@ void KellyErrorEstimator<dim>::estimate (const DoFHandler<dim>               &do
                                         const std::vector<bool>                  &component_mask,
                                         const Function<dim>                 *coefficients,
                                         const unsigned int                   n_threads,
-                                         const unsigned int       subdomain_id)
+                                         const unsigned int       subdomain_id,
+                                         const unsigned int       material_id)
 {
   Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));  
   static const MappingQ1<dim> mapping;
   estimate(mapping, dof_handler, quadrature, neumann_bc, solutions,
-          errors, component_mask, coefficients, n_threads, subdomain_id);
+          errors, component_mask, coefficients, n_threads,
+           subdomain_id, material_id);
 }
 
 
@@ -1295,6 +1336,7 @@ estimate<InputVector > (const Mapping<deal_II_dimension>      &,    \
           const std::vector<bool> &,    \
           const Function<deal_II_dimension>     *,    \
           const unsigned int       , \
+          const unsigned int       , \
           const unsigned int);    \
     \
 template    \
@@ -1308,6 +1350,7 @@ estimate<InputVector > (const DoFHandler<deal_II_dimension>   &,    \
           const std::vector<bool> &,    \
           const Function<deal_II_dimension>     *,    \
           const unsigned int       , \
+          const unsigned int       , \
           const unsigned int);    \
         \
 template    \
@@ -1322,6 +1365,7 @@ estimate<InputVector > (const Mapping<deal_II_dimension>          &,    \
           const std::vector<bool>     &,    \
           const Function<deal_II_dimension>         *,    \
           const unsigned int           , \
+          const unsigned int           , \
           const unsigned int);    \
     \
 template    \
@@ -1335,6 +1379,7 @@ estimate<InputVector > (const DoFHandler<deal_II_dimension>       &,    \
           const std::vector<bool>     &,    \
           const Function<deal_II_dimension>         *,    \
           const unsigned int           , \
+          const unsigned int           , \
           const unsigned int)
 
 INSTANTIATE(Vector<double>);
index fc91f2879bb928c6deda7e4ebcdc88b7d163b072..779ed6ff1ab19d91349c70d6d6deddeb6d94db7f 100644 (file)
@@ -258,6 +258,14 @@ inconvenience this causes.
 <h3>deal.II</h3>
 
 <ol>
+  <li> <p>
+       Improved: The <code>KellyErrorEstimator</code> class now also allows to
+       select the cells on which it is supposed to work by giving a material
+       ID, instead of only a subdomain ID.
+       <br> 
+       (WB 2006/1/4)
+       </p>
+
   <li> <p>
        Improved: A new <code>TriaAccessor::ExcCellHasNoChildren</code>
        exception will be raised if the <code
diff --git a/tests/bits/error_estimator_02.cc b/tests/bits/error_estimator_02.cc
new file mode 100644 (file)
index 0000000..6e938f4
--- /dev/null
@@ -0,0 +1,215 @@
+//----------------------------  error_estimator_02.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2004, 2006 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  error_estimator_02.cc  ---------------------------
+
+// same as error_estimator_02, but check for material_id instead of
+// subdomain_id
+
+
+
+#include "../tests.h"
+#include <base/logstream.h>
+#include <base/function_lib.h>
+#include <base/quadrature_lib.h>
+#include <lac/vector.h>
+#include <grid/tria.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_accessor.h>
+#include <grid/grid_generator.h>
+#include <dofs/dof_handler.h>
+#include <fe/fe_q.h>
+#include <fe/mapping_q.h>
+#include <numerics/vectors.h>
+#include <numerics/error_estimator.h>
+
+#include <fstream>
+
+
+template<int dim>
+class MySquareFunction : public Function<dim>
+{
+  public:
+    MySquareFunction () : Function<dim>(2) {};
+    
+    virtual double value (const Point<dim>   &p,
+                         const unsigned int  component) const
+      {        return (component+1)*p.square(); };
+    
+    virtual void   vector_value (const Point<dim>   &p,
+                                Vector<double>     &values) const
+      { values(0) = value(p,0);
+       values(1) = value(p,1); };
+};
+
+
+
+template <int dim>
+Quadrature<dim-1> &
+get_q_face (Function<dim>&)
+{
+  static QGauss4<dim-1> q;
+  return q;
+}
+
+
+Quadrature<0> &
+get_q_face (Function<1>&)
+{
+  Quadrature<0> *q = 0;
+  return *q;
+}
+
+
+
+template <int dim>
+void make_mesh (Triangulation<dim> &tria)
+{
+  
+  GridGenerator::hyper_cube(tria, -1, 1);
+
+                                   // refine the mesh in a random way so as to
+                                   // generate as many cells with
+                                   // hanging nodes as possible
+  tria.refine_global (4-dim);
+  const double steps[4] = { /*d=0*/ 0, 7, 3, 3 };
+  for (unsigned int i=0; i<steps[dim]; ++i)
+    {
+      typename Triangulation<dim>::active_cell_iterator
+        cell = tria.begin_active();
+      for (unsigned int index=0; cell != tria.end(); ++cell, ++index)
+        if (index % (3*dim) == 0)
+          cell->set_refine_flag();
+      tria.execute_coarsening_and_refinement ();
+    }
+
+                                  // we now have a number of cells,
+                                  // flag them with some material
+                                  // ids based on their position, in
+                                  // particular we take the quadrant
+                                  // (octant)
+  typename Triangulation<dim>::active_cell_iterator
+    cell = tria.begin_active (),
+    endc = tria.end ();
+  for (; cell!=endc; ++cell)
+    {
+      unsigned int material = 0;
+      for (unsigned int d=0; d<dim; ++d)
+       if (cell->center()(d) > 0)
+         material |= (1<<d);
+      Assert (material < (1<<dim), ExcInternalError());
+
+      cell->set_material_id (material);
+    }
+}
+
+
+
+
+template <int dim>
+void
+check ()
+{
+  Functions::CosineFunction<dim> function;
+  
+  Triangulation<dim> tria;
+  make_mesh (tria);
+  
+  FE_Q<dim> element(3);
+  DoFHandler<dim> dof(tria);
+  dof.distribute_dofs(element);
+
+  MappingQ<dim> mapping(3);
+  Quadrature<dim-1> &q_face = get_q_face(function);
+
+  std::map<unsigned char,const Function<dim>*> neumann_bc;
+  neumann_bc[0] = &function;
+  
+  Vector<double> v (dof.n_dofs());
+  VectorTools::interpolate (mapping, dof, function, v);
+
+  Vector<float> error1 (tria.n_active_cells());
+  Vector<float> error2 (tria.n_active_cells());
+
+                                   // compute error by looking at all cells at
+                                   // once and output this as the base line
+                                   // results. scale results so that they show
+                                   // up in the output file as a reasonable
+                                   // number
+  KellyErrorEstimator<dim>::estimate (mapping, dof, q_face, neumann_bc,
+                                      v, error1);
+  const double scaling_factor = 10000.*error1.size()/error1.l1_norm();
+  error1 *= scaling_factor;
+    
+  deallog << "Estimated error indicators:" << std::endl;
+  for (unsigned int i=0; i<error1.size(); ++i)
+    deallog << error1(i) << std::endl;
+
+                                   // then do the same with different
+                                   // material ids and add up the result
+  for (unsigned int material=0; material<(1<<dim); ++material)
+    {
+      deallog << "Material id=" << material << std::endl;
+        
+      Vector<float> this_error (tria.n_active_cells());
+      KellyErrorEstimator<dim>::estimate (mapping, dof, q_face, neumann_bc,
+                                          v, this_error,
+                                          std::vector<bool>(), 0,
+                                          multithread_info.n_default_threads,
+                                          deal_II_numbers::invalid_unsigned_int,
+                                          material);
+      this_error *= scaling_factor;
+
+                                       // copy the result into error2. since
+                                       // every invokation of the kelly
+                                       // estimator should only operate on
+                                       // the cells of one material,
+                                       // whenever there is something in
+                                       // this_error, the corresponding
+                                       // entry in error2 should still be
+                                       // empty
+      for (unsigned int i=0; i<this_error.size(); ++i)
+        {
+          deallog << i << ' ' << this_error(i) << std::endl;
+            
+          Assert ((this_error(i)==0) || (error2(i)==0),
+                  ExcInternalError());
+          if (this_error(i) != 0)
+            error2(i) = this_error(i);
+        }
+    }
+
+                                   // now compare the results of the two
+                                   // computations
+  Assert (error1 == error2, ExcInternalError());
+
+  deallog << "OK" << std::endl;
+}
+
+
+int main ()
+{
+  std::ofstream logfile ("error_estimator_02/output");
+  logfile.precision (2);
+  logfile.setf(std::ios::fixed);  
+  deallog.attach(logfile);
+  deallog.depth_console (0);
+
+  deallog.push ("1d");
+  check<1> ();
+  deallog.pop ();
+  deallog.push ("2d");
+  check<2> ();
+  deallog.pop ();
+  deallog.push ("3d");
+  check<3> ();
+  deallog.pop ();
+}

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

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.