]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Implement multithreaded version of the error estimator.
authorrichter <richter@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 1 Feb 2000 10:35:37 +0000 (10:35 +0000)
committerrichter <richter@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 1 Feb 2000 10:35:37 +0000 (10:35 +0000)
git-svn-id: https://svn.dealii.org/trunk@2302 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/numerics/error_estimator.h
deal.II/deal.II/source/numerics/error_estimator.cc

index 69018e41e53cdb43f65207b2552b9cf44b93c97d..8eeb90c9529d53f3ed4ef20fd13646534158b804 100644 (file)
 #include <dofs/dof_handler.h>
 #include <map>
 
-
+                                // if multithreaded set number
+                                // of threads to use.
+                                // The default number of threads can
+                                // be set during the compiling with
+                                // the flag N_THREADS.
+#ifdef DEAL_II_USE_MT
+ #ifndef N_THREADS
+  #define N_THREADS 4
+ #endif
+#else
+ #ifdef N_THREADS
+  #undef N_THREADS
+ #endif
+ #define N_THREADS 1
+#endif
 
 
 /**
 template <int dim>
 class KellyErrorEstimator {
   public:
-                                    /**
+
+                                    /**
                                      * Declare a data type which denotes a
                                      * mapping between a boundary indicator
                                      * and the function denoting the boundary
@@ -180,9 +195,10 @@ class KellyErrorEstimator {
                                      * Only one boundary function may be given
                                      * for each boundary indicator, which is
                                      * guaranteed by the #map# data type.
-                                     */
+                                     */    
     typedef map<unsigned char,const Function<dim>*> FunctionMap;
 
+    
                                     /**
                                      * Implementation of the error
                                      * estimator described above. You
@@ -219,14 +235,28 @@ class KellyErrorEstimator {
                                      * bit-vector with all-set
                                      * entries, or an empty
                                      * bit-vector.
+                                     *
+                                     * The estimator supports
+                                     * multithreading and splits
+                                     * the cells to 4 (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 only
+                                     * implemented in two and three
+                                     * dimensions. 
                                      */
+
     static void estimate (const DoFHandler<dim>   &dof,
                          const Quadrature<dim-1> &quadrature,
                          const FunctionMap       &neumann_bc,
                          const Vector<double>    &solution,
                          Vector<float>           &error,
-                         const vector<bool>      &component_mask = vector<bool>(),
-                         const Function<dim>     *coefficients   = 0);
+                         const vector<bool>      &component_mask_ = vector<bool>(),
+                         const Function<dim>     *coefficients   = 0,
+                         unsigned int            n_threads=N_THREADS);
     
                                     /**
                                      * Exception
@@ -245,8 +275,12 @@ class KellyErrorEstimator {
                                      */
     DeclException0 (ExcInvalidBoundaryFunction);
     
+
+
   private:
-                                    /**
+
+
+                                    /**
                                      * Declare a data type to represent the
                                      * mapping between faces and integrated
                                      * jumps of gradients. See the general
@@ -254,12 +288,107 @@ class KellyErrorEstimator {
                                      */
     typedef map<typename DoFHandler<dim>::face_iterator,double> FaceIntegrals;
 
+
                                     /**
                                      * Redeclare an active cell iterator.
                                      * This is simply for convenience.
                                      */
     typedef DoFHandler<dim>::active_cell_iterator active_cell_iterator;
+
+
+                                    /**
+                                     * All data needed by the several functions
+                                     * of the error estimator is gathered in
+                                     * this struct. It is passed as an reference
+                                     * to the seperate functions in the object.
+                                     * Global data is not possible because no
+                                     * real object is created.
+                                     */
+    struct Data
+    {
+       const DoFHandler<dim>   &dof;
+       const Quadrature<dim-1> &quadrature;
+       const FunctionMap       &neumann_bc;
+       const Vector<double>    &solution;
+       vector<bool>            component_mask;
+       const Function<dim>     *coefficients;
+       unsigned int            n_threads;
+    
+       DoFHandler<dim>::active_cell_iterator endc;
+       unsigned int            n_components;
+       unsigned int            n_q_points;
+       FaceIntegrals           face_integrals;  
+
+                                        /**
+                                         * A vector to store the jump of the
+                                         * normal vectors in the quadrature
+                                         * points.
+                                         * There is one vector for every
+                                         * thread used in the estimator.
+                                         * The allocation of memory has to
+                                         * be global to enable fast use
+                                         * of multithreading     
+                                         */ 
+       vector< vector<vector<double> > >         phi;
+                                        /**
+                                         * A vector for the gradients of
+                                         * the finite element function
+                                         * on one cell
+                                         *
+                                         * 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.
+                                         */
+       vector< vector<vector<Tensor<1,dim> > > > psi;
+                                        /**
+                                         * The same vector for a neighbor cell
+                                         */
+       vector< vector<vector<Tensor<1,dim> > > > neighbor_psi;
+                                        /**
+                                         * The normal vectors of the finite
+                                         * element function on one face
+                                         */
+       vector< vector<Point<dim> > > normal_vectors;
+       
+       vector< vector<double> >                  coefficient_values1;
+       vector< vector<Vector<double> > >         coefficient_values;
+
+       vector< vector<double> >                  JxW_values;
+
+                                        /**
+                                         * A constructor of the
+                                         * class Data. All variables are
+                                         * passed as references.
+                                         */
+       Data(const DoFHandler<dim>   &dof,
+            const Quadrature<dim-1> &quadrature,
+            const FunctionMap       &neumann_bc,
+            const Vector<double>    &solution,
+            vector<bool>            component_mask_,
+            const Function<dim>     *coefficients,
+            unsigned int            n_threads);
+    };
+
     
+                                    /**
+                                     * Computates the error on all cells
+                                     * of the domain with the number n,
+                                     * satisfying
+                                     * #n=this_thread (mod n_threads)#
+                                     * This enumeration is choosen to
+                                     * generate a random distribution
+                                     * of all cells.
+                                     *
+                                     * This function is only needed in
+                                     * two or three dimensions.
+                                     * The Errorestimator in one dimension
+                                     * is implemented seperatly.
+                                     */
+    static void * estimate_some(Data &data, unsigned int this_thread);
+                               
                                     /**
                                      * Actually do the computation on a face
                                      * which has no hanging nodes (it is
@@ -278,17 +407,16 @@ class KellyErrorEstimator {
                                      * to avoid ending up with a function
                                      * of 500 lines of code.
                                      */
-    static void integrate_over_regular_face (const active_cell_iterator &cell,
+
+
+    static void integrate_over_regular_face (Data &data,
+                                            int this_thread,
+                                            const active_cell_iterator &cell,
                                             const unsigned int   face_no,
-                                            const FunctionMap   &neumann_bc,
-                                            const unsigned int   n_q_points,
                                             FEFaceValues<dim>   &fe_face_values_cell,
-                                            FEFaceValues<dim>   &fe_face_values_neighbor,
-                                            FaceIntegrals       &face_integrals,
-                                            const Vector<double>&solution,
-                                            const vector<bool>  &component_mask,
-                                            const Function<dim> *coefficient);
-
+                                            FEFaceValues<dim>   &fe_face_values_neighbor);
+    
+    
                                     /**
                                      * The same applies as for the function
                                      * above, except that integration is
@@ -297,15 +425,12 @@ class KellyErrorEstimator {
                                      * so that the integration is a bit more
                                      * complex.
                                      */
-    static void integrate_over_irregular_face (const active_cell_iterator &cell,
+    static void integrate_over_irregular_face (Data &data,
+                                              int this_thread,
+                                              const active_cell_iterator &cell,
                                               const unsigned int   face_no,
-                                              const unsigned int    n_q_points,
                                               FEFaceValues<dim>    &fe_face_values,
-                                              FESubfaceValues<dim> &fe_subface_values,
-                                              FaceIntegrals        &face_integrals,
-                                              const Vector<double> &solution,
-                                              const vector<bool>   &component_mask,
-                                              const Function<dim>  *coefficient);
+                                              FESubfaceValues<dim> &fe_subface_values);
 };
 
 
@@ -316,3 +441,5 @@ class KellyErrorEstimator {
 /* end of #ifndef __error_estimator_H */
 #endif
 /*----------------------------   error_estimator.h     ---------------------------*/
+
+
index 0a9be5f7dd0f57f828f91e8b851f0ed7ed75a248..f4ca76d3636b1b57ce4659171e825653942fcd85 100644 (file)
 #include <algorithm>
 #include <cmath>
 #include <vector>
+#include <base/timer.h>
+
+                                // if multithreaded include
+                                // ThreadManager
+#ifdef DEAL_II_USE_MT
+#include <base/thread_manager.h>
+#endif
+
+#if deal_II_dimension == 1
+
+template <>
+KellyErrorEstimator<1>::Data::Data(const DoFHandler<1>     &,
+                                  const Quadrature<0>     &,
+                                  const FunctionMap       &,
+                                  const Vector<double>    &,
+                                  vector<bool>            ,
+                                  const Function<1>       *,
+                                  unsigned int            ):
+               dof(* static_cast <const DoFHandler<1> *> (0)),
+               quadrature(* static_cast <const Quadrature<0> *> (0)),
+               neumann_bc(* static_cast <const FunctionMap *> (0)),
+               solution(* static_cast <const Vector<double> *> (0))
+{
+  Assert (false, ExcInternalError());
+}
+
+#else
+
+template <int dim>
+KellyErrorEstimator<dim>::Data::Data(const DoFHandler<dim>   &dof,
+                                    const Quadrature<dim-1> &quadrature,
+                                    const FunctionMap       &neumann_bc,
+                                    const Vector<double>    &solution,
+                                    vector<bool>            component_mask_,
+                                    const Function<dim>     *coefficients,
+                                    unsigned int            n_threads):
+               dof(dof),
+               quadrature(quadrature),
+               neumann_bc(neumann_bc),
+               solution(solution),
+               coefficients(coefficients),
+               n_threads(n_threads)
+{
+  n_components = dof.get_fe().n_components();
+  
+                                  // if no mask given: treat all components
+  component_mask = ((component_mask_.size() == 0)    ?
+                   vector<bool>(n_components, true) :
+                   component_mask_);
+  
+  Assert (component_mask.size() == n_components, ExcInvalidComponentMask());
+  Assert (count(component_mask.begin(), component_mask.end(), true) > 0,
+         ExcInvalidComponentMask());
+  
+  Assert ((coefficients == 0) ||
+         (coefficients->n_components == n_components) ||
+         (coefficients->n_components == 1),
+         ExcInvalidCoefficient());
+  
+  Assert (neumann_bc.find(255) == neumann_bc.end(),
+         ExcInvalidBoundaryIndicator());
+  
+  for (FunctionMap::const_iterator i=neumann_bc.begin(); i!=neumann_bc.end(); ++i)
+    Assert (i->second->n_components == n_components, ExcInvalidBoundaryFunction());
+  
+                                  // the last cell, often needed
+  endc=dof.end();
+  
+  n_q_points=quadrature.n_quadrature_points;
+  
+                                  // Init the size of a lot of vectors
+                                  // needed in the calculations once
+                                  // per thread.
+  phi.resize(n_threads);
+  psi.resize(n_threads);
+  neighbor_psi.resize(n_threads);
+  normal_vectors.resize(n_threads);
+  coefficient_values1.resize(n_threads);
+  coefficient_values.resize(n_threads);
+  JxW_values.resize(n_threads);
+  
+  for (unsigned int t=0;t<n_threads;++t)
+    {
+      phi[t].resize(n_q_points);
+      psi[t].resize(n_q_points);
+      neighbor_psi[t].resize(n_q_points);
+      normal_vectors[t].resize(n_q_points);
+      coefficient_values1[t].resize(n_q_points);
+      coefficient_values[t].resize(n_q_points);
+      JxW_values[t].resize(n_q_points);
+      
+      for (unsigned int qp=0;qp<n_q_points;++qp)
+       {
+         phi[t][qp].resize(n_components);
+         psi[t][qp].resize(n_components);
+         neighbor_psi[t][qp].resize(n_components);
+         coefficient_values[t][qp].reinit(n_components);
+       }
+    }
+}
+
+#endif
 
 
 inline static double sqr (const double x) {
@@ -24,9 +126,15 @@ inline static double sqr (const double x) {
 };
 
 
-
 #if deal_II_dimension == 1
 
+template <>
+void * KellyErrorEstimator<1>::estimate_some (Data &,const unsigned int )
+{
+  Assert (false, ExcInternalError() );
+  return 0;
+}
+
 template <>
 void KellyErrorEstimator<1>::estimate (const DoFHandler<1>  &dof,
                                       const Quadrature<0>  &,
@@ -34,7 +142,8 @@ void KellyErrorEstimator<1>::estimate (const DoFHandler<1>  &dof,
                                       const Vector<double> &solution,
                                       Vector<float>        &error,
                                       const vector<bool>   &component_mask_,
-                                      const Function<1>    *coefficient)
+                                      const Function<1>    *coefficient,
+                                      unsigned int)
 {
   const unsigned int n_components = dof.get_fe().n_components();
 
@@ -174,56 +283,14 @@ void KellyErrorEstimator<1>::estimate (const DoFHandler<1>  &dof,
     };
 };
        
-#endif
 
+                                // #if deal_II_dimension !=1
+#else
 
 
 template <int dim>
-void KellyErrorEstimator<dim>::estimate (const DoFHandler<dim>    &dof,
-                                        const Quadrature<dim-1>  &quadrature,
-                                        const FunctionMap        &neumann_bc,
-                                        const Vector<double>     &solution,
-                                        Vector<float>            &error,
-                                        const vector<bool>       &component_mask_,
-                                        const Function<dim>      *coefficient)
+void * KellyErrorEstimator<dim>::estimate_some (Data &data, unsigned int this_thread) 
 {
-  const unsigned int n_components = dof.get_fe().n_components();
-
-                                  // if no mask given: treat all components
-  vector<bool> component_mask ((component_mask_.size() == 0)    ?
-                              vector<bool>(n_components, true) :
-                              component_mask_);
-  Assert (component_mask.size() == n_components, ExcInvalidComponentMask());
-  Assert (count(component_mask.begin(), component_mask.end(), true) > 0,
-         ExcInvalidComponentMask());
-  
-  Assert ((coefficient == 0) ||
-         (coefficient->n_components == n_components) ||
-         (coefficient->n_components == 1),
-         ExcInvalidCoefficient());
-
-  Assert (neumann_bc.find(255) == neumann_bc.end(),
-         ExcInvalidBoundaryIndicator());
-
-  for (FunctionMap::const_iterator i=neumann_bc.begin(); i!=neumann_bc.end(); ++i)
-    Assert (i->second->n_components == n_components, ExcInvalidBoundaryFunction());
-
-                                  // create a map of integrals indexed by
-                                  // the corresponding face. In this map
-                                  // we store the integrated jump of the
-                                  // gradient for each face. By doing so,
-                                  // we can check whether we have already
-                                  // integrated along this face by testing
-                                  // whether the respective face is already
-                                  // a key in this map.
-                                  // At the end of the function, we again
-                                  // loop over the cells and collect the
-                                  // conrtibutions of the different faces
-                                  // of the cell.
-  FaceIntegrals face_integrals;
-  
-                                  // number of integration points per face
-  const unsigned int n_q_points = quadrature.n_quadrature_points;
   
                                   // make up a fe face values object for the
                                   // restriction of the finite element function
@@ -235,124 +302,235 @@ void KellyErrorEstimator<dim>::estimate (const DoFHandler<dim>    &dof,
                                   // need not compute all values on the
                                   // neighbor cells, so using two objects
                                   // gives us a performance gain).
-  FEFaceValues<dim> fe_face_values_cell (dof.get_fe(), quadrature,
+  FEFaceValues<dim> fe_face_values_cell (data.dof.get_fe(),
+                                        data.quadrature,
                                         UpdateFlags(update_gradients      |
                                                     update_JxW_values     |
-                                                    ((!neumann_bc.empty() ||
-                                                      (coefficient != 0))  ?
+                                                    ((!data.neumann_bc.empty() ||
+                                                      (data.coefficients != 0))  ?
                                                      update_q_points : 0) |
                                                     update_normal_vectors)); 
-  FEFaceValues<dim> fe_face_values_neighbor (dof.get_fe(), quadrature, update_gradients);
-  FESubfaceValues<dim> fe_subface_values (dof.get_fe(), quadrature, update_gradients);
+  FEFaceValues<dim> fe_face_values_neighbor (data.dof.get_fe(),
+                                            data.quadrature,
+                                            update_gradients);
+  FESubfaceValues<dim> fe_subface_values (data.dof.get_fe(),
+                                         data.quadrature,
+                                         update_gradients);
   
-                                  // loop over all cells
-  const active_cell_iterator endc = dof.end();
-  for (active_cell_iterator cell = dof.begin_active(); cell!=endc; ++cell)
-                                    // loop over all faces of this cell
-    for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
-      {
-                                        // if we already visited this
-                                        // face: do nothing
-       if (face_integrals.find(cell->face(face_no)) !=
-           face_integrals.end())
-         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 (cell->at_boundary(face_no) == false)
-         if (cell->neighbor(face_no)->level() < cell->level())
-           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.
-       const unsigned char boundary_indicator
-         = cell->face(face_no)->boundary_indicator();
-       if ((boundary_indicator != 255) &&
-           neumann_bc.find(boundary_indicator)==neumann_bc.end()) 
-         {
-           face_integrals[cell->face(face_no)] = 0;
+  
+  DoFHandler<dim>::active_cell_iterator cell=data.dof.begin_active();
+
+                                  // calculate the start cell for this
+                                  // thread. the enumeration is choosen
+                                  // in this strange way to generate a
+                                  // "random" distribution of the cells.
+                                  // if the sequence of the iterator would
+                                  // be used, the threads would take widely
+                                  // spread times to calculate their cells.
+  for (unsigned int t=0;t<this_thread;++t,++cell);
+
+                                  // loop over all cells for this thread
+                                  // the iteration of cell is done at the end
+  for (;cell!=data.endc;)
+    {
+      
+                                      // loop over all faces of this cell
+      for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
+       {
+                                          // if we already visited this
+                                          // face: do nothing
+         if (data.face_integrals[cell->face(face_no)] >=0)
            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 (cell->at_boundary(face_no) == false)
+           if (cell->neighbor(face_no)->level() < cell->level())
+             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.
+         const unsigned char boundary_indicator
+           = cell->face(face_no)->boundary_indicator();
+         if ((boundary_indicator != 255) &&
+             data.neumann_bc.find(boundary_indicator)==data.neumann_bc.end()) 
+           {
+             data.face_integrals[cell->face(face_no)] = 0;
+             continue;
+           };
+         
+         
+         
+         if (cell->face(face_no)->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
+           integrate_over_regular_face (data,
+                                        this_thread,
+                                        cell, face_no,
+                                        fe_face_values_cell,
+                                        fe_face_values_neighbor);
+         
+         else
+                                            // otherwise we need to do some
+                                            // special computations which do
+                                            // not fit into the framework of
+                                            // the above function
+           integrate_over_irregular_face (data,
+                                          this_thread,cell, face_no,
+                                          fe_face_values_cell,
+                                          fe_subface_values);
+       };
 
-       
-       if (cell->face(face_no)->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
-         integrate_over_regular_face (cell, face_no, neumann_bc,
-                                      n_q_points,
-                                      fe_face_values_cell,
-                                      fe_face_values_neighbor,
-                                      face_integrals,
+                                      // next cell in this thread
+      for (unsigned int t=0;((t<data.n_threads)&&(cell!=data.endc));++t,++cell) {};
+    };
+  return 0;
+}
+
+template <int dim>
+void KellyErrorEstimator<dim>::estimate (const DoFHandler<dim>   &dof,
+                                        const Quadrature<dim-1> &quadrature,
+                                        const FunctionMap       &neumann_bc,
+                                        const Vector<double>    &solution,
+                                        Vector<float>           &error,
+                                        const vector<bool>      &component_mask,
+                                        const Function<dim>     *coefficients,
+                                        unsigned int            n_threads)
+{
+                                  // if NOT multithreaded, set n_threads to one
+#ifndef DEAL_II_USE_MT
+  n_threads = 1;
+#endif
+  
+                                  // all the data needed in the error-
+                                  // estimator is gathered in this stuct.
+  KellyErrorEstimator<dim>::Data data (dof,
+                                      quadrature,
+                                      neumann_bc,
                                       solution,
                                       component_mask,
-                                      coefficient);
-       else
-                                          // otherwise we need to do some
-                                          // special computations which do
-                                          // not fit into the framework of
-                                          // the above function
-         integrate_over_irregular_face (cell, face_no,
-                                        n_q_points,
-                                        fe_face_values_cell,
-                                        fe_subface_values,
-                                        face_integrals, solution,
-                                        component_mask,
-                                        coefficient);
-      };
-
+                                      coefficients,
+                                      n_threads);
+                                  // map of integrals indexed by
+                                  // the corresponding face. In this map
+                                  // we store the integrated jump of the
+                                  // gradient for each face. By doing so,
+                                  // we can check whether we have already
+                                  // integrated along this face by testing
+                                  // whether the respective face is already
+                                  // a key in this map.
+                                  // At the end of the function, we again
+                                  // loop over the cells and collect the
+                                  // conrtibutions of the different faces
+                                  // of the cell.
+                                  // the values for all faces are set to
+                                  // -10e20. It would cost a lot of time
+                                  // to syncronisise the initialisation
+                                  // of the map in multithreaded mode.
+                                  // negative value indicates that the
+                                  // face is not calculated.
+  
+  for (active_cell_iterator cell=data.dof.begin_active(); cell!=data.endc; ++cell)
+    for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
+      data.face_integrals[cell->face(face_no)]=-10e20;
+  
+  
+                                  // split all cells into threads
+                                  // if multithreading is used
+#ifdef DEAL_II_USE_MT
 
-                                  // finally add up the contributions of the
+  ThreadManager thread_manager;
+  
+                                  // all data needed to start
+                                  // one thread is gathered
+                                  // in this struct.
+  typedef ThreadManager::Fun_Data2
+    <class Data &, unsigned int > FunData;
+  
+                                  // One struct of this type for every thread
+  vector<FunData>
+    fun_data (data.n_threads,
+             FunData (data,0,&KellyErrorEstimator::estimate_some));
+  
+  
+                                   // get start cells for each thread
+  for (unsigned int l=0;l<data.n_threads;++l)
+    {
+      fun_data[l].arg2=l;
+    };
+    
+                                  // now spawn the threads
+  for (unsigned int i=0;i<data.n_threads; ++i)
+    {
+      thread_manager.spawn(&fun_data[i],THR_SCOPE_SYSTEM | THR_DETACHED);
+    };
+                                  // wait for all threads to return
+  thread_manager.wait();
+  
+                                  // ... ifdef DEAL_II_USE_MT
+#else
+                                  // just one thread, calculate
+                                  // error on all cells
+  KellyErrorEstimator::estimate_some(data,0);
+  
+#endif
+  
+  
+                                  // finally add up the contributions of the
                                   // faces for each cell
   
                                   // reserve one slot for each cell and set
                                   // it to zero
-  error.reinit (dof.get_tria().n_active_cells());
+  error.reinit (data.dof.get_tria().n_active_cells());
+  for (unsigned int i=0;i<data.dof.get_tria().n_active_cells();++i)
+    error(i)=0;
 
   unsigned int present_cell=0;
-  for (active_cell_iterator cell=dof.begin_active(); cell!=endc; ++cell, ++present_cell)
+  
+  for (active_cell_iterator cell=data.dof.begin_active();
+       cell!=data.endc;
+       ++cell, ++present_cell)
     {
                                       // loop over all faces of this cell
-      for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no) 
+      for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell;
+          ++face_no)
        {
-         Assert (face_integrals.find(cell->face(face_no)) !=
-                 face_integrals.end(),
-                 ExcInternalError());
-         error(present_cell) += (face_integrals[cell->face(face_no)] *
+         Assert(data.face_integrals.find(cell->face(face_no)) != data.face_integrals.end(),
+                ExcInternalError());
+         error(present_cell) += (data.face_integrals[cell->face(face_no)] *
                                  cell->diameter() / 24);
        };
-          
       error(present_cell) = sqrt(error(present_cell));
     };
 };
 
+#endif
 
 
 #if deal_II_dimension == 1
 
 template <>
-void KellyErrorEstimator<1>::integrate_over_regular_face (const active_cell_iterator &,
-                                                         const unsigned int      ,
-                                                         const FunctionMap      &,
+void KellyErrorEstimator<1>::integrate_over_regular_face (Data &,
+                                                         int ,
+                                                         const active_cell_iterator &,
                                                          const unsigned int      ,
                                                          FEFaceValues<1>        &,
-                                                         FEFaceValues<1>        &,
-                                                         FaceIntegrals          &,
-                                                         const Vector<double>   &,
-                                                         const vector<bool>     &,
-                                                         const Function<1>      *) {
+                                                         FEFaceValues<1>        &)
+{
   Assert (false, ExcInternalError());
 };
 
@@ -360,15 +538,13 @@ void KellyErrorEstimator<1>::integrate_over_regular_face (const active_cell_iter
 
 template <>
 void KellyErrorEstimator<1>::
-integrate_over_irregular_face (const active_cell_iterator &,
-                              const unsigned int          ,
+integrate_over_irregular_face (Data &,
+                              int,
+                              const active_cell_iterator &,
                               const unsigned int          ,
                               FEFaceValues<1>            &,
-                              FESubfaceValues<1>         &,
-                              FaceIntegrals              &,
-                              const Vector<double>       &,
-                              const vector<bool>         &,
-                              const Function<1>          *) {
+                              FESubfaceValues<1>         &)
+{
   Assert (false, ExcInternalError());
 };
 
@@ -378,37 +554,22 @@ integrate_over_irregular_face (const active_cell_iterator &,
 
 template <int dim>
 void KellyErrorEstimator<dim>::
-integrate_over_regular_face (const active_cell_iterator &cell,
-                            const unsigned int          face_no,
-                            const FunctionMap          &neumann_bc,
-                            const unsigned int          n_q_points,
+integrate_over_regular_face (Data &data,
+                            int this_thread,
+                            const active_cell_iterator &cell,
+                            const unsigned int         face_no,
                             FEFaceValues<dim>          &fe_face_values_cell,
-                            FEFaceValues<dim>          &fe_face_values_neighbor,
-                            FaceIntegrals              &face_integrals,
-                            const Vector<double>       &solution,
-                            const vector<bool>         &component_mask,
-                            const Function<dim>        *coefficient)
+                            FEFaceValues<dim>          &fe_face_values_neighbor)
 {
-  const unsigned int n_components = component_mask.size();
   const DoFHandler<dim>::face_iterator face = cell->face(face_no);
   
                                   // initialize data of the restriction
                                   // of this cell to the present face
   fe_face_values_cell.reinit (cell, face_no);
   
-                                  // 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
-  vector<vector<Tensor<1,dim> > > psi(n_q_points, vector<Tensor<1,dim> >(n_components));
-  fe_face_values_cell.get_function_grads (solution, psi);
+                                  // get gradients of the finite element
+                                  // function on this cell
+  fe_face_values_cell.get_function_grads (data.solution, data.psi[this_thread]);
   
                                   // now compute over the other side of
                                   // the face
@@ -435,18 +596,15 @@ integrate_over_regular_face (const active_cell_iterator &cell,
                                       // function of #neighbor# to the
                                       // common face.
       fe_face_values_neighbor.reinit (neighbor, neighbor_neighbor);
-
-                                      // get a list of the gradients of
-                                      // the finite element solution
-                                      // restricted to the neighbor cell
-      vector<vector<Tensor<1,dim> > > neighbor_psi (n_q_points,
-                                                   vector<Tensor<1,dim> >(n_components));
-      fe_face_values_neighbor.get_function_grads (solution, neighbor_psi);
+      
+                                      // get gradients on neighbor cell
+      fe_face_values_neighbor.get_function_grads (data.solution,
+                                                 data.neighbor_psi[this_thread]);
       
                                       // 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)
-         psi[p][component] -= neighbor_psi[p][component];
+      for (unsigned int component=0; component<data.n_components; ++component)
+       for (unsigned int p=0; p<data.n_q_points; ++p)
+         data.psi[this_thread][p][component] -= data.neighbor_psi[this_thread][p][component];
     };
   
   
@@ -467,43 +625,41 @@ integrate_over_regular_face (const active_cell_iterator &cell,
                                   // taking that of the neighbor
                                   // would only change the sign. We take
                                   // the outward normal.
-                                  //
-                                  // let phi be the name of the integrand
-  vector<vector<double> > phi(n_q_points, vector<double>(n_components));
-  const vector<Point<dim> > &normal_vectors(fe_face_values_cell.
-                                           get_normal_vectors());
-
-  for (unsigned int component=0; component<n_components; ++component)
-    for (unsigned int point=0; point<n_q_points; ++point)
-      phi[point][component] = psi[point][component]*normal_vectors[point];
-
+  
+  data.normal_vectors[this_thread]=fe_face_values_cell.get_normal_vectors();
+  
+  for (unsigned int component=0; component<data.n_components; ++component)
+    for (unsigned int point=0; point<data.n_q_points; ++point)
+      data.phi[this_thread][point][component] = data.psi[this_thread][point][component]*
+                                          data.normal_vectors[this_thread][point];
+  
                                   // if a coefficient was given: use that
                                   // to scale the jump in the gradient
-  if (coefficient != 0)
+  if (data.coefficients != 0)
     {
                                       // scalar coefficient
-      if (coefficient->n_components == 1)
+      if (data.coefficients->n_components == 1)
        {
-         vector<double>      coefficient_values (n_q_points);
-         coefficient->value_list (fe_face_values_cell.get_quadrature_points(),
-                                  coefficient_values);
-         for (unsigned int component=0; component<n_components; ++component)
-           for (unsigned int point=0; point<n_q_points; ++point)
-             phi[point][component] *= coefficient_values[point];
+         
+         data.coefficients->value_list (fe_face_values_cell.get_quadrature_points(),
+                                        data.coefficient_values1[this_thread]);
+           for (unsigned int component=0; component<data.n_components; ++component)
+             for (unsigned int point=0; point<data.n_q_points; ++point)
+               data.phi[this_thread][point][component] *=
+                 data.coefficient_values1[this_thread][point];
        }
       else
-                                        // vector-valued coefficient
+                                          // vector-valued coefficient
        {
-         vector<Vector<double> > coefficient_values (n_q_points,
-                                                     Vector<double>(n_components));
-         coefficient->vector_value_list (fe_face_values_cell.get_quadrature_points(),
-                                         coefficient_values);
-         for (unsigned int component=0; component<n_components; ++component)
-           for (unsigned int point=0; point<n_q_points; ++point)
-             phi[point][component] *= coefficient_values[point](component);
+         data.coefficients->vector_value_list (fe_face_values_cell.get_quadrature_points(),
+                                          data.coefficient_values[this_thread]);
+         for (unsigned int component=0; component<data.n_components; ++component)
+           for (unsigned int point=0; point<data.n_q_points; ++point)
+             data.phi[this_thread][point][component] *=
+               data.coefficient_values[this_thread][point](component);
        };
     };
-      
+  
   
   if (face->at_boundary() == true)
                                     // neumann boundary face. compute
@@ -511,20 +667,21 @@ integrate_over_regular_face (const active_cell_iterator &cell,
                                     // derivative and boundary function
     {
       const unsigned char boundary_indicator = face->boundary_indicator();
-
-      Assert (neumann_bc.find(boundary_indicator) != neumann_bc.end(),
+      
+      Assert (data.neumann_bc.find(boundary_indicator) != data.neumann_bc.end(),
              ExcInternalError ());
                                       // get the values of the boundary
                                       // function at the quadrature
                                       // points
-      vector<Vector<double> > g(n_q_points, Vector<double>(n_components));
-      neumann_bc.find(boundary_indicator)->second
+      
+      vector<Vector<double> > g(data.n_q_points, Vector<double>(data.n_components));
+      data.neumann_bc.find(boundary_indicator)->second
        ->vector_value_list (fe_face_values_cell.get_quadrature_points(),
                             g);
       
-      for (unsigned int component=0; component<n_components; ++component)
-       for (unsigned int point=0; point<n_q_points; ++point)
-         phi[point][component] -= g[point](component);
+      for (unsigned int component=0; component<data.n_components; ++component)
+       for (unsigned int point=0; point<data.n_q_points; ++point)
+         data.phi[this_thread][point][component] -= g[point](component);
     };
   
   
@@ -535,18 +692,19 @@ integrate_over_regular_face (const active_cell_iterator &cell,
                                   // each component being the
                                   // mentioned value at one of the
                                   // quadrature points
-  const vector<double> &JxW_values = fe_face_values_cell.get_JxW_values();
+
+  data.JxW_values[this_thread] = fe_face_values_cell.get_JxW_values();
   
                                   // take the square of the phi[i]
                                   // for integration, and sum up
   double face_integral = 0;
-  for (unsigned int component=0; component<n_components; ++component)
-    if (component_mask[component] == true)
-      for (unsigned int p=0; p<n_q_points; ++p)
-       face_integral += sqr(phi[p][component]) *
-                        JxW_values[p];
+  for (unsigned int component=0; component<data.n_components; ++component)
+    if (data.component_mask[component] == true)
+      for (unsigned int p=0; p<data.n_q_points; ++p)
+       face_integral += sqr(data.phi[this_thread][p][component]) *
+                        data.JxW_values[this_thread][p];
   
-  face_integrals[face] = face_integral;
+  data.face_integrals[face] = face_integral;
 };
 
 
@@ -554,18 +712,13 @@ integrate_over_regular_face (const active_cell_iterator &cell,
 
 template <int dim>
 void KellyErrorEstimator<dim>::
-integrate_over_irregular_face (const active_cell_iterator &cell,
+integrate_over_irregular_face (Data &data,
+                              int this_thread,
+                              const active_cell_iterator &cell,
                               const unsigned int          face_no,
-                              const unsigned int          n_q_points,
                               FEFaceValues<dim>          &fe_face_values,
-                              FESubfaceValues<dim>       &fe_subface_values,
-                              FaceIntegrals              &face_integrals,
-                              const Vector<double>       &solution,
-                              const vector<bool>         &component_mask,
-                              const Function<dim>        *coefficient)
+                              FESubfaceValues<dim>       &fe_subface_values)
 {
-  const unsigned int n_components = component_mask.size();
-
   const DoFHandler<dim>::cell_iterator neighbor = cell->neighbor(face_no);
   Assert (neighbor.state() == valid, ExcInternalError());
   Assert (neighbor->has_children(), ExcInternalError());
@@ -580,8 +733,7 @@ integrate_over_irregular_face (const active_cell_iterator &cell,
                                   // finite element, and the first
                                   // index the number of the
                                   // quadrature point
-  vector<vector<Tensor<1,dim> > > psi(n_q_points, vector<Tensor<1,dim> >(n_components));
-
+  
                                   // store which number #cell# has in the
                                   // list of neighbors of #neighbor#
   unsigned int neighbor_neighbor;
@@ -606,26 +758,26 @@ integrate_over_irregular_face (const active_cell_iterator &cell,
       Assert (!neighbor->child(GeometryInfo<dim>::
                               child_cell_on_face(neighbor_neighbor,subface_no))->has_children(),
              ExcInternalError());
-            
+      
                                       // restrict the finite element on the
                                       // present cell to the subface and
                                       // store the gradient of the solution
                                       // in psi
       fe_subface_values.reinit (cell, face_no, subface_no);
-      fe_subface_values.get_function_grads (solution, psi);
+      fe_subface_values.get_function_grads (data.solution, data.psi[this_thread]);
 
                                       // restrict the finite element on the
                                       // neighbor cell to the common #subface#.
                                       // store the gradient in #neighbor_psi#
-      vector<vector<Tensor<1,dim> > > neighbor_psi (n_q_points,
-                                                   vector<Tensor<1,dim> >(n_components));
+      
       fe_face_values.reinit (neighbor_child, neighbor_neighbor);
-      fe_face_values.get_function_grads (solution, neighbor_psi);
+      fe_face_values.get_function_grads (data.solution, data.neighbor_psi[this_thread]);
       
                                       // 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)
-         psi[p][component] -= neighbor_psi[p][component];
+      for (unsigned int component=0; component<data.n_components; ++component)
+       for (unsigned int p=0; p<data.n_q_points; ++p)
+         data.psi[this_thread][p][component] -=
+           data.neighbor_psi[this_thread][p][component];
 
                                       // note that unlike for the
                                       // case of regular faces
@@ -645,53 +797,54 @@ integrate_over_irregular_face (const active_cell_iterator &cell,
                                       // the outward normal.
                                       //
                                       // let phi be the name of the integrand
-      vector<vector<double> > phi(n_q_points, vector<double>(n_components));
-      const vector<Point<dim> > &normal_vectors(fe_face_values.
-                                               get_normal_vectors());
-  
-      for (unsigned int component=0; component<n_components; ++component)
-       for (unsigned int point=0; point<n_q_points; ++point)
-         phi[point][component] = psi[point][component]*normal_vectors[point];
+     
+      data.normal_vectors[this_thread]=fe_face_values.get_normal_vectors();
+      
+      
+      for (unsigned int component=0; component<data.n_components; ++component)
+       for (unsigned int point=0; point<data.n_q_points; ++point)
+         data.phi[this_thread][point][component] =
+           data.psi[this_thread][point][component]*
+           data.normal_vectors[this_thread][point];
       
                                       // if a coefficient was given: use that
                                       // to scale the jump in the gradient
-      if (coefficient != 0)
+      if (data.coefficients != 0)
        {
                                           // scalar coefficient
-         if (coefficient->n_components == 1)
+         if (data.coefficients->n_components == 1)
            {
-             vector<double>      coefficient_values (n_q_points);
-             coefficient->value_list (fe_face_values.get_quadrature_points(),
-                                      coefficient_values);
-             for (unsigned int component=0; component<n_components; ++component)
-               for (unsigned int point=0; point<n_q_points; ++point)
-                 phi[point][component] *= coefficient_values[point];
+             data.coefficients->value_list (fe_face_values.get_quadrature_points(),
+                                            data.coefficient_values1[this_thread]);
+             for (unsigned int component=0; component<data.n_components; ++component)
+               for (unsigned int point=0; point<data.n_q_points; ++point)
+                 data.phi[this_thread][point][component] *=
+                   data.coefficient_values1[this_thread][point];
            }
          else
                                             // vector-valued coefficient
            {
-             vector<Vector<double> > coefficient_values (n_q_points,
-                                                         Vector<double>(n_components));
-             coefficient->vector_value_list (fe_face_values.get_quadrature_points(),
-                                             coefficient_values);
-             for (unsigned int component=0; component<n_components; ++component)
-               for (unsigned int point=0; point<n_q_points; ++point)
-                 phi[point][component] *= coefficient_values[point](component);
+             data.coefficients->vector_value_list (fe_face_values.get_quadrature_points(),
+                                                   data.coefficient_values[this_thread]);
+             for (unsigned int component=0; component<data.n_components; ++component)
+               for (unsigned int point=0; point<data.n_q_points; ++point)
+                 data.phi[this_thread][point][component] *=
+                   data.coefficient_values[this_thread][point](component);
            };
        };
-
-      const vector<double> &JxW_values = fe_face_values.get_JxW_values();
+      
+      data.JxW_values[this_thread] = fe_face_values.get_JxW_values();
       
                                       // take the square of the phi[i]
                                       // for integration, and sum up
       double face_integral = 0;
-      for (unsigned int component=0; component<n_components; ++component)
-       if (component_mask[component] == true)
-         for (unsigned int p=0; p<n_q_points; ++p)
-           face_integral += sqr(phi[p][component]) *
-                            JxW_values[p];
+      for (unsigned int component=0; component<data.n_components; ++component)
+       if (data.component_mask[component] == true)
+         for (unsigned int p=0; p<data.n_q_points; ++p)
+           face_integral += sqr(data.phi[this_thread][p][component]) *
+                            data.JxW_values[this_thread][p];
 
-      face_integrals[neighbor_child->face(neighbor_neighbor)] = face_integral;
+      data.face_integrals[neighbor_child->face(neighbor_neighbor)] = face_integral;
     };
 
   
@@ -704,16 +857,21 @@ integrate_over_irregular_face (const active_cell_iterator &cell,
   for (unsigned int subface_no=0; subface_no<GeometryInfo<dim>::subfaces_per_face;
        ++subface_no) 
     {
-      Assert (face_integrals.find(face->child(subface_no)) !=
-             face_integrals.end(),
+      Assert (data.face_integrals.find(face->child(subface_no)) !=
+             data.face_integrals.end(),
+             ExcInternalError());
+      Assert (data.face_integrals[face->child(subface_no)]>=0,
              ExcInternalError());
-      sum += face_integrals[face->child(subface_no)];
+      sum += data.face_integrals[face->child(subface_no)];
     };
-  face_integrals[face] = sum;
+
+  data.face_integrals[face] = sum;
+
 };
 
 
 
 
 // explicit instantiations
+
 template class KellyErrorEstimator<deal_II_dimension>;

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.