]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Merged from trunk
authorheltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 17 Oct 2013 15:29:01 +0000 (15:29 +0000)
committerheltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 17 Oct 2013 15:29:01 +0000 (15:29 +0000)
git-svn-id: https://svn.dealii.org/branches/branch_manifold_id@31284 0785d39b-7218-0410-832d-ea1e28bc413d

17 files changed:
deal.II/cmake/setup_compiler_flags_intel.cmake
deal.II/doc/doxygen/headers/global_dof_index.h
deal.II/doc/news/changes.h
deal.II/examples/step-20/doc/intro.dox
deal.II/examples/step-9/step-9.cc
deal.II/include/deal.II/base/types.h
deal.II/include/deal.II/base/vectorization.h
deal.II/include/deal.II/fe/mapping.h
deal.II/include/deal.II/lac/parallel_block_vector.h
deal.II/include/deal.II/lac/sparsity_pattern.h
deal.II/include/deal.II/numerics/data_out.h
deal.II/include/deal.II/numerics/derivative_approximation.h
deal.II/source/base/data_out_base.cc
deal.II/source/grid/grid_tools.cc
deal.II/source/numerics/data_out.cc
deal.II/source/numerics/derivative_approximation.cc
deal.II/source/numerics/matrix_tools.cc

index 912798af4768f68a4c49ce1faeb4a7138cac68af..98afb1943ba43de15c05da90bb684812a19a0a2e 100644 (file)
@@ -62,6 +62,14 @@ ENABLE_IF_SUPPORTED(CMAKE_CXX_FLAGS "-w2")
 #   -w68   integer conversion resulted in a change of sign
 #          (triggers a lot in functionparser)
 #   -w175  subscript out of range
+#   -w135  class template "dealii::FE_Q_Base<POLY, dim, spacedim>"
+#          has no member "Implementation"
+#          (the compiler is objectively wrong since the warning
+#           triggers also on code of the form
+#           class FE_Q_Base {
+#             struct Implementation; // forward declaration
+#             friend struct Implementation;
+#           };)
 #   -w177  declared but not referenced
 #   -w279  controlling expression is constant
 #   -w327  NULL reference is not allowed
@@ -78,6 +86,7 @@ ENABLE_IF_SUPPORTED(CMAKE_CXX_FLAGS "-w2")
 #   -w2536 type qualifiers are meaningless here
 #
 ENABLE_IF_SUPPORTED(CMAKE_CXX_FLAGS "-wd68")
+ENABLE_IF_SUPPORTED(CMAKE_CXX_FLAGS "-wd135")
 ENABLE_IF_SUPPORTED(CMAKE_CXX_FLAGS "-wd175")
 ENABLE_IF_SUPPORTED(CMAKE_CXX_FLAGS "-wd177")
 ENABLE_IF_SUPPORTED(CMAKE_CXX_FLAGS "-wd279")
index 51ebce36d7a40e91181a8ea3cfc2f617f38326c4..2b5d84bddd3c25f9f3f78bd06d007787572c79f6 100644 (file)
 /**
  * @page GlobalDoFIndex When to use types::global_dof_index instead of unsigned int
  *
- * When the 64-bits version of deal.II is used, it becomes necessary to
- * declare as types::global_dof_index (unsigned long long int) instead of
- * unsigned int. Here, we want to clarify when types::global_dof_index must be
- * used.
+ * deal.II can be configured to use 64-bit indices for degrees of freedom,
+ * rather than the usual unsigned integers that default to 32-bit on most
+ * current systems. This is necessary since we want to be able to solve
+ * problems with more than 4 billion unknowns (the limit of what can be
+ * represented with 32-bit unsigned integers). At the same time, we do not
+ * want to indiscriminately replace all integers in deal.II with 64-bit
+ * versions, since this would increase memory use in many places where we
+ * represent quantities that will most definitely not be larger than 4 billion.
+ *
+ * The data type we define for these indices to keep the bulk
+ * of the code base free of <code>#ifdef</code>s is types::global_dof_index.
+ * If deal.II is configured as normal, this type is <code>unsigned int</code>,
+ * but can be switched to <code>unsigned long long int</code> if the right
+ * flag is provided (see the ReadMe file). This page is intended to clarify
+ * when types::global_dof_index must be used and when one can use a regular
+ * unsigned integer:
  *
  * <dl>
  *
  *
  * <dt class="glossary">@anchor GlobalDoFIndexCell <b>Cell</b></dt>
  * <dd>
- * The ID of cell is not unique. Cells with different levels of refinement
+ * The ID of cell is not unique: Cells with different levels of refinement
  * and/or on different processors can have the same ID. Thus, all the data
- * associated to cells can be unsigned int.
+ * associated to cells can be unsigned int because on a single processor,
+ * one one mesh level, there will definitely not be more than 4 billion
+ * cells.
  * </dd>
  *
  * <dt class="glossary">@anchor GlobalDoFIndexDoFHandler
  * <b>DoFHandler</b></dt>
  * <dd>
- * The ID of each degree of freedom is unique. Therefore, degrees of freedom
+ * The ID of each degree of freedom is unique in a parallel computation.
+ * Therefore, degrees of freedom
  * are types::global_dof_index.
  * </dd>
  *
@@ -52,7 +67,8 @@
  * <dd>
  * The numbers of row and column are types::global_dof_index even if it is not
  * expected that someone will create a FullMatrix with so many entries.
- * However, ConstraintMatrix is templated on the matrix type and thus, the
+ * However, some functions of ConstraintMatrix are templated on the matrix
+ * type and thus, the
  * size of a FullMatrix has to be of the same type than the size of
  * SparseMatrix.
  * </dd>
  * <dt class="glossary">@anchor GlobalDoFIndexSparseMatrix
  * <b>SparseMatrix</b></dt>
  * <dd>
- * The size of SparseMatrix can be arbitrary large therefore,
- * types::global_do_index is used. However, even for a large complex problem we
- * can solve now, there is no reason for the number of non-zero entries in a 
- * sparse matrix to go over four billions. Thus, we still use unsigned int
- * for, e.g., row_lengths in the object.
+ * The size of SparseMatrix can be arbitrary large and it is conceivable that
+ * with sufficient memory on a single node, one may generate a matrix with
+ * more than 4 billion rows or columns. Therefore, types::global_dof_index is
+ * used. However, even for a large complex problem we can solve now, it is not
+ * reasonable to expect the number of non-zero entries in a sparse matrix to
+ * go over four billion. Thus, we still use unsigned int for, e.g.,
+ * SparsityPattern::row_lengths and similar functions.
  * </dd>
  *
  * </dl>
index 7cdcdd827eb4a460dd3de67e5709f24bb3ba8491..54d0a2412bf08d629557c48f4ba06edb117c3cbf 100644 (file)
@@ -145,6 +145,12 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+  <li>
+  Fixed: When deriving from DataOut to filter the cells where output is generated, there were two different bugs that result in segmentation faults or wrong cells written (example, step-18).
+  <br>
+  (Timo Heister, 2013/10/16)
+  </li>
+
   <li>
   New: GridIn::read_vtk() reads 2d and 3d meshes in VTK format.
   <br>
index dd4c3829ca99ed8709fa3db0c3c072a25a4b476d..13476bc2015c1e15361cceba3713a80212c96def 100644 (file)
@@ -400,10 +400,10 @@ also a full matrix. On the other hand, the CG algorithm doesn't require
 us to actually have a representation of $S$, it is sufficient to form
 matrix-vector products with it. We can do so in steps: to compute $Sv=B^TM^{-1}Bv=B^T(M^{-1}(Bv))$, we
 <ol>
- <li> form $w = T v$;
+ <li> form $w = B v$;
  <li> solve $My = w$ for $y=M^{-1}w$, using the CG method applied to the
   positive definite and symmetric mass matrix $M$;
- <li> form $z=B^Ty$ to obtain $Sv=z$.
+ <li> form $z=B^Ty$ to obtain $z=Sv$.
 </ol>
 Note how we evaluate the expression $B^TM^{-1}Bv$ right to left to
 avoid matrix-matrix products; this way, all we have to do is evaluate
index 4a9e28152e026260c5ba533e6b56fb9a268d703c..73b56ce91278dda09f40589c70d5b7a050db3db4 100644 (file)
@@ -54,7 +54,7 @@
 // <code>multithread_info</code> of that type) which can be used to query the
 // number of processors in your system, which is often useful when deciding
 // how many threads to start in parallel.
-#include <deal.II/base/thread_management.h>
+#include <deal.II/base/work_stream.h>
 #include <deal.II/base/multithread_info.h>
 
 // The next new include file declares a base class <code>TensorFunction</code>
@@ -74,6 +74,25 @@ namespace Step9
 {
   using namespace dealii;
 
+  namespace Assembler
+  {
+    struct Scratch
+    {
+      Scratch() {};
+    };
+
+    struct CopyData
+    {
+      CopyData() {};
+
+      unsigned int dofs_per_cell;
+      std::vector<types::global_dof_index> local_dof_indices;
+      // We declare cell matrix and cell right hand side...
+      FullMatrix<double> cell_matrix;
+      Vector<double>  cell_rhs;
+    };
+  }
+
   // @sect3{AdvectionProblem class declaration}
 
   // Following we declare the main class of this program. It is very much
@@ -107,8 +126,10 @@ namespace Step9
     // in the module, there are other, and possibly better suited, ways to
     // achieve the same goal.
     void assemble_system ();
-    void assemble_system_interval (const typename DoFHandler<dim>::active_cell_iterator &begin,
-                                   const typename DoFHandler<dim>::active_cell_iterator &end);
+    void build_local_system (typename DoFHandler<dim>::active_cell_iterator const &cell,
+        Assembler::Scratch &scratch,Assembler::CopyData &copy_data);
+    void copy_local_to_global (Assembler::CopyData const &copy_data);
+                                   
 
     // The following functions again are as in previous examples, as are the
     // subsequent variables.
@@ -128,23 +149,6 @@ namespace Step9
 
     Vector<double>       solution;
     Vector<double>       system_rhs;
-
-    // When assembling the matrix in parallel, we have to synchronize when
-    // several threads attempt to write the local contributions of a cell to
-    // the global matrix at the same time. This is done using a
-    // <code>Mutex</code>, which is an object that can be owned by only one
-    // thread at a time. If a thread wants to write to the matrix, it has to
-    // acquire this lock (if it is presently owned by another thread, then it
-    // has to wait), then write to the matrix and finally release the
-    // lock. Note that if the library was not compiled to support
-    // multithreading (which you have to specify at the time you call the
-    // <code>./configure</code> script in the top-level directory), then
-    // the actual data type of the typedef
-    // <code>Threads::Mutex</code> is a dummy class that provides all the
-    // functions needed for a mutex, but does nothing when they are called;
-    // this is reasonable, of course, since if only one thread is running at a
-    // time, there is no need to synchronize with other threads.
-    Threads::Mutex     assembler_lock;
   };
 
 
@@ -526,7 +530,7 @@ namespace Step9
     // is returned by n_threads(). This
     // is also queried by functions inside the library to determine
     // how many threads they shall create.
-    const unsigned int n_threads = multithread_info.n_threads();
+    
     // It is worth noting, however, that this setup determines the load
     // distribution onto processor in a static way: it does not take into
     // account that some other part of our program may also be running
@@ -559,7 +563,6 @@ namespace Step9
     // object. Likewise, the function <code>join</code> that is supposed to
     // wait for all spawned threads to return, returns immediately, as there
     // can't be any threads running.
-    Threads::ThreadGroup<> threads;
 
     // Now we have to split the range of cells into chunks of approximately
     // the same size. Each thread will then assemble the local contributions
@@ -584,22 +587,21 @@ namespace Step9
     // <code>raw_iterator</code>), and in this case the C++ language requires
     // us to specify the template type explicitly. For brevity, we first
     // typedef this data type to an alias.
+
     typedef typename DoFHandler<dim>::active_cell_iterator active_cell_iterator;
-    std::vector<std::pair<active_cell_iterator,active_cell_iterator> >
-    thread_ranges
-      = Threads::split_range<active_cell_iterator> (dof_handler.begin_active (),
-                                                    dof_handler.end (),
-                                                    n_threads);
 
     // Finally, for each of the chunks of iterators we have computed, start
     // one thread (or if not in multithread mode: execute assembly on these
     // chunks sequentially). This is done using the following sequence of
     // function calls:
-    for (unsigned int thread=0; thread<n_threads; ++thread)
-      threads += Threads::new_thread (&AdvectionProblem<dim>::assemble_system_interval,
-                                      *this,
-                                      thread_ranges[thread].first,
-                                      thread_ranges[thread].second);
+
+    Assembler::Scratch scratch;
+    Assembler::CopyData copy_data;
+    WorkStream::run(dof_handler.begin_active(),dof_handler.end(),*this,
+        &AdvectionProblem::build_local_system,&AdvectionProblem::copy_local_to_global,
+        scratch,copy_data);
+
+
     // The reasons and internal workings of these functions can be found in
     // the report on the subject of multithreading, which is available online
     // as well. Suffice it to say that we create a new thread that calls the
@@ -624,7 +626,6 @@ namespace Step9
     //
     // Again, if the library was not configured to use multithreading, then
     // no threads can run in parallel and the function returns immediately.
-    threads.join_all ();
 
 
     // After the matrix has been assembled in parallel, we still have to
@@ -647,8 +648,9 @@ namespace Step9
   template <int dim>
   void
   AdvectionProblem<dim>::
-  assemble_system_interval (const typename DoFHandler<dim>::active_cell_iterator &begin,
-                            const typename DoFHandler<dim>::active_cell_iterator &end)
+  build_local_system (typename DoFHandler<dim>::active_cell_iterator const &cell,
+      Assembler::Scratch &scratch,
+      Assembler::CopyData &copy_data)
   {
     // First of all, we will need some objects that describe boundary values,
     // right hand side function and the advection field. As we will only
@@ -681,17 +683,17 @@ namespace Step9
                                       update_JxW_values | update_normal_vectors);
 
     // Then we define some abbreviations to avoid unnecessarily long lines:
-    const unsigned int   dofs_per_cell   = fe.dofs_per_cell;
+    copy_data.dofs_per_cell = fe.dofs_per_cell;
     const unsigned int   n_q_points      = quadrature_formula.size();
     const unsigned int   n_face_q_points = face_quadrature_formula.size();
 
     // We declare cell matrix and cell right hand side...
-    FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
-    Vector<double>       cell_rhs (dofs_per_cell);
+    copy_data.cell_matrix = FullMatrix<double> (copy_data.dofs_per_cell, copy_data.dofs_per_cell);
+    copy_data.cell_rhs = Vector<double> (copy_data.dofs_per_cell);
 
     // ... an array to hold the global indices of the degrees of freedom of
     // the cell on which we are presently working...
-    std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+    copy_data.local_dof_indices.resize(copy_data.dofs_per_cell);
 
     // ... and array in which the values of right hand side, advection
     // direction, and boundary values will be stored, for cell and face
@@ -701,117 +703,118 @@ namespace Step9
     std::vector<double>         face_boundary_values (n_face_q_points);
     std::vector<Tensor<1,dim> > face_advection_directions (n_face_q_points);
 
-    // Then we start the main loop over the cells:
-    typename DoFHandler<dim>::active_cell_iterator cell;
-    for (cell=begin; cell!=end; ++cell)
-      {
-        // First clear old contents of the cell contributions...
-        cell_matrix = 0;
-        cell_rhs = 0;
-
-        // ... then initialize the <code>FEValues</code> object...
-        fe_values.reinit (cell);
-
-        // ... obtain the values of right hand side and advection directions
-        // at the quadrature points...
-        advection_field.value_list (fe_values.get_quadrature_points(),
-                                    advection_directions);
-        right_hand_side.value_list (fe_values.get_quadrature_points(),
-                                    rhs_values);
-
-        // ... set the value of the streamline diffusion parameter as
-        // described in the introduction...
-        const double delta = 0.1 * cell->diameter ();
-
-        // ... and assemble the local contributions to the system matrix and
-        // right hand side as also discussed above:
-        for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
-          for (unsigned int i=0; i<dofs_per_cell; ++i)
-            {
-              for (unsigned int j=0; j<dofs_per_cell; ++j)
-                cell_matrix(i,j) += ((advection_directions[q_point] *
-                                      fe_values.shape_grad(j,q_point)   *
-                                      (fe_values.shape_value(i,q_point) +
-                                       delta *
-                                       (advection_directions[q_point] *
-                                        fe_values.shape_grad(i,q_point)))) *
-                                     fe_values.JxW(q_point));
-
-              cell_rhs(i) += ((fe_values.shape_value(i,q_point) +
-                               delta *
-                               (advection_directions[q_point] *
-                                fe_values.shape_grad(i,q_point))        ) *
-                              rhs_values[q_point] *
-                              fe_values.JxW (q_point));
-            };
 
-        // Besides the cell terms which we have build up now, the bilinear
-        // form of the present problem also contains terms on the boundary of
-        // the domain. Therefore, we have to check whether any of the faces of
-        // this cell are on the boundary of the domain, and if so assemble the
-        // contributions of this face as well. Of course, the bilinear form
-        // only contains contributions from the <code>inflow</code> part of
-        // the boundary, but to find out whether a certain part of a face of
-        // the present cell is part of the inflow boundary, we have to have
-        // information on the exact location of the quadrature points and on
-        // the direction of flow at this point; we obtain this information
-        // using the FEFaceValues object and only decide within the main loop
-        // whether a quadrature point is on the inflow boundary.
-        for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
-          if (cell->face(face)->at_boundary())
-            {
-              // Ok, this face of the present cell is on the boundary of the
-              // domain. Just as for the usual FEValues object which we have
-              // used in previous examples and also above, we have to
-              // reinitialize the FEFaceValues object for the present face:
-              fe_face_values.reinit (cell, face);
-
-              // For the quadrature points at hand, we ask for the values of
-              // the inflow function and for the direction of flow:
-              boundary_values.value_list (fe_face_values.get_quadrature_points(),
-                                          face_boundary_values);
-              advection_field.value_list (fe_face_values.get_quadrature_points(),
-                                          face_advection_directions);
-
-              // Now loop over all quadrature points and see whether it is on
-              // the inflow or outflow part of the boundary. This is
-              // determined by a test whether the advection direction points
-              // inwards or outwards of the domain (note that the normal
-              // vector points outwards of the cell, and since the cell is at
-              // the boundary, the normal vector points outward of the domain,
-              // so if the advection direction points into the domain, its
-              // scalar product with the normal vector must be negative):
-              for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
-                if (fe_face_values.normal_vector(q_point) *
-                    face_advection_directions[q_point]
-                    < 0)
-                  // If the is part of the inflow boundary, then compute the
-                  // contributions of this face to the global matrix and right
-                  // hand side, using the values obtained from the
-                  // FEFaceValues object and the formulae discussed in the
-                  // introduction:
-                  for (unsigned int i=0; i<dofs_per_cell; ++i)
-                    {
-                      for (unsigned int j=0; j<dofs_per_cell; ++j)
-                        cell_matrix(i,j) -= (face_advection_directions[q_point] *
-                                             fe_face_values.normal_vector(q_point) *
-                                             fe_face_values.shape_value(i,q_point) *
-                                             fe_face_values.shape_value(j,q_point) *
-                                             fe_face_values.JxW(q_point));
-
-                      cell_rhs(i) -= (face_advection_directions[q_point] *
-                                      fe_face_values.normal_vector(q_point) *
-                                      face_boundary_values[q_point]         *
-                                      fe_face_values.shape_value(i,q_point) *
-                                      fe_face_values.JxW(q_point));
-                    };
-            };
+    // ... then initialize the <code>FEValues</code> object...
+    fe_values.reinit (cell);
+
+    // ... obtain the values of right hand side and advection directions
+    // at the quadrature points...
+    advection_field.value_list (fe_values.get_quadrature_points(),
+                                advection_directions);
+    right_hand_side.value_list (fe_values.get_quadrature_points(),
+                                rhs_values);
+
+    // ... set the value of the streamline diffusion parameter as
+    // described in the introduction...
+    const double delta = 0.1 * cell->diameter ();
+
+    // ... and assemble the local contributions to the system matrix and
+    // right hand side as also discussed above:
+    for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+      for (unsigned int i=0; i<copy_data.dofs_per_cell; ++i)
+        {
+          for (unsigned int j=0; j<copy_data.dofs_per_cell; ++j)
+            copy_data.cell_matrix(i,j) += ((advection_directions[q_point] *
+                                            fe_values.shape_grad(j,q_point)   *
+                                            (fe_values.shape_value(i,q_point) +
+                                             delta *
+                                             (advection_directions[q_point] *
+                                              fe_values.shape_grad(i,q_point)))) *
+                                           fe_values.JxW(q_point));
+
+          copy_data.cell_rhs(i) += ((fe_values.shape_value(i,q_point) +
+                                     delta *
+                                     (advection_directions[q_point] *
+                                      fe_values.shape_grad(i,q_point))        ) *
+                                    rhs_values[q_point] *
+                                    fe_values.JxW (q_point));
+        };
+
+    // Besides the cell terms which we have build up now, the bilinear
+    // form of the present problem also contains terms on the boundary of
+    // the domain. Therefore, we have to check whether any of the faces of
+    // this cell are on the boundary of the domain, and if so assemble the
+    // contributions of this face as well. Of course, the bilinear form
+    // only contains contributions from the <code>inflow</code> part of
+    // the boundary, but to find out whether a certain part of a face of
+    // the present cell is part of the inflow boundary, we have to have
+    // information on the exact location of the quadrature points and on
+    // the direction of flow at this point; we obtain this information
+    // using the FEFaceValues object and only decide within the main loop
+    // whether a quadrature point is on the inflow boundary.
+    for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+      if (cell->face(face)->at_boundary())
+        {
+          // Ok, this face of the present cell is on the boundary of the
+          // domain. Just as for the usual FEValues object which we have
+          // used in previous examples and also above, we have to
+          // reinitialize the FEFaceValues object for the present face:
+          fe_face_values.reinit (cell, face);
+
+          // For the quadrature points at hand, we ask for the values of
+          // the inflow function and for the direction of flow:
+          boundary_values.value_list (fe_face_values.get_quadrature_points(),
+                                      face_boundary_values);
+          advection_field.value_list (fe_face_values.get_quadrature_points(),
+                                      face_advection_directions);
+
+          // Now loop over all quadrature points and see whether it is on
+          // the inflow or outflow part of the boundary. This is
+          // determined by a test whether the advection direction points
+          // inwards or outwards of the domain (note that the normal
+          // vector points outwards of the cell, and since the cell is at
+          // the boundary, the normal vector points outward of the domain,
+          // so if the advection direction points into the domain, its
+          // scalar product with the normal vector must be negative):
+          for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
+            if (fe_face_values.normal_vector(q_point) *
+                face_advection_directions[q_point]
+                < 0)
+              // If the is part of the inflow boundary, then compute the
+              // contributions of this face to the global matrix and right
+              // hand side, using the values obtained from the
+              // FEFaceValues object and the formulae discussed in the
+              // introduction:
+              for (unsigned int i=0; i<copy_data.dofs_per_cell; ++i)
+                {
+                  for (unsigned int j=0; j<copy_data.dofs_per_cell; ++j)
+                    copy_data.cell_matrix(i,j) -= (face_advection_directions[q_point] *
+                                                   fe_face_values.normal_vector(q_point) *
+                                                   fe_face_values.shape_value(i,q_point) *
+                                                   fe_face_values.shape_value(j,q_point) *
+                                                   fe_face_values.JxW(q_point));
+
+                  copy_data.cell_rhs(i) -= (face_advection_directions[q_point] *
+                                            fe_face_values.normal_vector(q_point) *
+                                            face_boundary_values[q_point]         *
+                                            fe_face_values.shape_value(i,q_point) *
+                                            fe_face_values.JxW(q_point));
+                };
+        };
+
+
+    // Now go on by transferring the local contributions to the system of
+    // equations into the global objects. The first step was to obtain the
+    // global indices of the degrees of freedom on this cell.
+    cell->get_dof_indices (copy_data.local_dof_indices);
+  }
 
 
-        // Now go on by transferring the local contributions to the system of
-        // equations into the global objects. The first step was to obtain the
-        // global indices of the degrees of freedom on this cell.
-        cell->get_dof_indices (local_dof_indices);
+  
+  template <int dim>
+  void
+  AdvectionProblem<dim>::copy_local_to_global (Assembler::CopyData const &copy_data)
+  {
 
         // Up until now we have not taken care of the fact that this function
         // might run more than once in parallel, as the operations above only
@@ -836,17 +839,15 @@ namespace Step9
         // matrix, and can do so freely. When finished, we release the lock
         // again so as to allow other threads to acquire it and write to the
         // matrix.
-        assembler_lock.acquire ();
-        for (unsigned int i=0; i<dofs_per_cell; ++i)
+        for (unsigned int i=0; i<copy_data.dofs_per_cell; ++i)
           {
-            for (unsigned int j=0; j<dofs_per_cell; ++j)
-              system_matrix.add (local_dof_indices[i],
-                                 local_dof_indices[j],
-                                 cell_matrix(i,j));
+            for (unsigned int j=0; j<copy_data.dofs_per_cell; ++j)
+              system_matrix.add (copy_data.local_dof_indices[i],
+                                 copy_data.local_dof_indices[j],
+                                 copy_data.cell_matrix(i,j));
 
-            system_rhs(local_dof_indices[i]) += cell_rhs(i);
+            system_rhs(copy_data.local_dof_indices[i]) += copy_data.cell_rhs(i);
           };
-        assembler_lock.release ();
         // At this point, the locked operations on the global matrix are done,
         // i.e. other threads can now enter into the protected section by
         // acquiring the lock. Two final notes are in place here, however:
@@ -863,8 +864,8 @@ namespace Step9
         // class; since all threads execute member functions of the same
         // object, they have the same <code>this</code> pointer and therefore
         // also operate on the same <code>lock</code>.
-      };
   }
+  
 
 
 
index 0da2ba52e6f21631bfe88cedde4eb8a69b2a85a1..7fa251bcad7a4e23fa284d9afccdfebccf1b9e7c 100644 (file)
@@ -70,6 +70,9 @@ namespace types
    *
    * The data type always indicates an
    * unsigned integer type.
+   *
+   * See the @ref GlobalDoFIndex page for guidance on when this
+   * type should or should not be used.
    */
   // TODO: we should check that unsigned long long int
   // has the same size as uint64_t
index f2c560f36c502895792e1d4f1cb744bc4063c412..ec48938935dd9d0fba2b98eca8b282a901e3c5f4 100644 (file)
@@ -1952,10 +1952,8 @@ DEAL_II_NAMESPACE_CLOSE
 namespace std
 {
   /**
-   * Computes the sine of a vectorized
-   * data field. The result is return as
-   * vectorized array in the form
-   * <tt>{sin(x[0]), sin(x[1]), ...,
+   * Computes the sine of a vectorized data field. The result is return as
+   * vectorized array in the form <tt>{sin(x[0]), sin(x[1]), ...,
    * sin(x[n_array_elements-1])}</tt>.
    *
    * @relates VectorizedArray
@@ -1974,10 +1972,8 @@ namespace std
 
 
   /**
-   * Computes the tangent of a vectorized
-   * data field. The result is return as
-   * vectorized array in the form
-   * <tt>{tan(x[0]), tan(x[1]), ...,
+   * Computes the tangent of a vectorized data field. The result is return as
+   * vectorized array in the form <tt>{tan(x[0]), tan(x[1]), ...,
    * tan(x[n_array_elements-1])}</tt>.
    *
    * @relates VectorizedArray
@@ -1995,10 +1991,8 @@ namespace std
 
 
   /**
-   * Computes the cosine of a vectorized
-   * data field. The result is return as
-   * vectorized array in the form
-   * <tt>{cos(x[0]), cos(x[1]), ...,
+   * Computes the cosine of a vectorized data field. The result is return as
+   * vectorized array in the form <tt>{cos(x[0]), cos(x[1]), ...,
    * cos(x[n_array_elements-1])}</tt>.
    *
    * @relates VectorizedArray
@@ -2016,10 +2010,8 @@ namespace std
 
 
   /**
-   * Computes the exponential of a vectorized
-   * data field. The result is return as
-   * vectorized array in the form
-   * <tt>{exp(x[0]), exp(x[1]), ...,
+   * Computes the exponential of a vectorized data field. The result is return
+   * as vectorized array in the form <tt>{exp(x[0]), exp(x[1]), ...,
    * exp(x[n_array_elements-1])}</tt>.
    *
    * @relates VectorizedArray
@@ -2037,10 +2029,8 @@ namespace std
 
 
   /**
-   * Computes the natural logarithm of a
-   * vectorized data field. The result is return
-   * as vectorized array in the form
-   * <tt>{log(x[0]), log(x[1]), ...,
+   * Computes the natural logarithm of a vectorized data field. The result is
+   * return as vectorized array in the form <tt>{log(x[0]), log(x[1]), ...,
    * log(x[n_array_elements-1])}</tt>.
    *
    * @relates VectorizedArray
@@ -2059,10 +2049,8 @@ namespace std
 
 
   /**
-   * Computes the square root of a vectorized
-   * data field. The result is return as
-   * vectorized array in the form
-   * <tt>{sqrt(x[0]), sqrt(x[1]), ...,
+   * Computes the square root of a vectorized data field. The result is return
+   * as vectorized array in the form <tt>{sqrt(x[0]), sqrt(x[1]), ...,
    * sqrt(x[n_array_elements-1])}</tt>.
    *
    * @relates VectorizedArray
@@ -2078,11 +2066,9 @@ namespace std
 
 
   /**
-   * Computes the absolute value (modulus) of a
-   * vectorized data field. The result is return
-   * as vectorized array in the form
-   * <tt>{abs(x[0]), abs(x[1]), ...,
-   * abs(x[n_array_elements-1])}</tt>.
+   * Computes the absolute value (modulus) of a vectorized data field. The
+   * result is return as vectorized array in the form <tt>{abs(x[0]),
+   * abs(x[1]), ..., abs(x[n_array_elements-1])}</tt>.
    *
    * @relates VectorizedArray
    */
@@ -2097,11 +2083,9 @@ namespace std
 
 
   /**
-   * Computes the componentwise maximum of two
-   * vectorized data fields. The result is
-   * return as vectorized array in the form
-   * <tt>{max(x[0],y[0]), max(x[1],y[1]),
-   * ...}</tt>.
+   * Computes the componentwise maximum of two vectorized data fields. The
+   * result is return as vectorized array in the form <tt>{max(x[0],y[0]),
+   * max(x[1],y[1]), ...}</tt>.
    *
    * @relates VectorizedArray
    */
@@ -2117,11 +2101,9 @@ namespace std
 
 
   /**
-   * Computes the componentwise minimum of two
-   * vectorized data fields. The result is
-   * return as vectorized array in the form
-   * <tt>{min(x[0],y[0]), min(x[1],y[1]),
-   * ...}</tt>.
+   * Computes the componentwise minimum of two vectorized data fields. The
+   * result is return as vectorized array in the form <tt>{min(x[0],y[0]),
+   * min(x[1],y[1]), ...}</tt>.
    *
    * @relates VectorizedArray
    */
index 647de86ce8adadf7029f82bd85395819efede146..4067d709832dc7c926b0a79ee61bc9bce6ce3d97 100644 (file)
@@ -591,7 +591,7 @@ public:
      mapping_contravariant_gradient() and mapping_piola_gradient()
      are only true as stated for linear mappings.
      If, for example, the mapping is bilinear then there is a missing
-     term associated with the derivative of of J.
+     term associated with the derivative of J.
   */
   virtual
   void
index 747d12a33eb1f00561054fc44c507392e4cfd81b..3be362a6359e73908308ccb266f51928be8ffe3e 100644 (file)
@@ -226,6 +226,41 @@ namespace parallel
       void reinit (const BlockVector<Number2> &V,
                    const bool                 fast=false);
 
+      /**
+       * This function copies the data that has accumulated in the data buffer
+       * for ghost indices to the owning processor. For the meaning of the
+       * argument @p operation, see the entry on @ref GlossCompress
+       * "Compressing distributed vectors and matrices" in the glossary.
+       *
+       * There are two variants for this function. If called with argument @p
+       * VectorOperation::add adds all the data accumulated in ghost elements
+       * to the respective elements on the owning processor and clears the
+       * ghost array afterwards. If called with argument @p
+       * VectorOperation::insert, a set operation is performed. Since setting
+       * elements in a vector with ghost elements is ambiguous (as one can set
+       * both the element on the ghost site as well as the owning site), this
+       * operation makes the assumption that all data is set correctly on the
+       * owning processor. Upon call of compress(VectorOperation::insert), all
+       * ghost entries are therefore simply zeroed out (using
+       * zero_ghost_values()). In debug mode, a check is performed that makes
+       * sure that the data set is actually consistent between processors,
+       * i.e., whenever a non-zero ghost element is found, it is compared to
+       * the value on the owning processor and an exception is thrown if these
+       * elements do not agree.
+       *
+       */
+      void compress (::dealii::VectorOperation::values operation);
+
+      /**
+       * Fills the data field for ghost indices with the values stored in the
+       * respective positions of the owning processor. This function is needed
+       * before reading from ghosts. The function is @p const even though
+       * ghost data is changed. This is needed to allow functions with a @p
+       * const vector to perform the data exchange without creating
+       * temporaries.
+       */
+      void update_ghost_values () const;
+
       /**
        * Return whether the vector contains only elements with value
        * zero. This function is mainly for internal consistency checks and
@@ -543,6 +578,29 @@ namespace parallel
     }
 
 
+
+    template <typename Number>
+    inline
+    void
+    BlockVector<Number>::compress (::dealii::VectorOperation::values operation)
+    {
+      for (unsigned int block=0; block<this->n_blocks(); ++block)
+        this->block(block).compress(operation);
+    }
+
+
+
+    template <typename Number>
+    inline
+    void
+    BlockVector<Number>::update_ghost_values () const
+    {
+      for (unsigned int block=0; block<this->n_blocks(); ++block)
+        this->block(block).update_ghost_values();
+    }
+
+
+
     template <typename Number>
     inline
     bool
index 532162f06b32d41013f9b629e8a89cf98e08b828..414e6ea92b88c2825ed6c9354dc49f862b0de435 100644 (file)
@@ -842,7 +842,9 @@ public:
 
   /**
    * Copy data from an object of type CompressedSparsityPattern,
-   * CompressedSetSparsityPattern or CompressedSimpleSparsityPattern.
+   * CompressedSetSparsityPattern or CompressedSimpleSparsityPattern. Although
+   * not a compressed sparsity pattern, this function is also instantiated
+   * if the argument is of type SparsityPattern (i.e., the current class).
    * Previous content of this object is lost, and the sparsity pattern is in
    * compressed mode afterwards.
    */
index 7dfe9c9bd4fce6f65d093789c87896480f7778ec..961fed7c20382de1a8032cba859ef8c1a4603904 100644 (file)
@@ -253,7 +253,8 @@ private:
   void build_one_patch (const std::pair<cell_iterator, unsigned int> *cell_and_index,
                         internal::DataOut::ParallelData<DH::dimension, DH::space_dimension> &data,
                         ::dealii::DataOutBase::Patch<DH::dimension, DH::space_dimension> &patch,
-                        const CurvedCellRegion curved_cell_region);
+                        const CurvedCellRegion curved_cell_region,
+                        std::vector<dealii::DataOutBase::Patch<DH::dimension, DH::space_dimension> > &patches);
 };
 
 
index 3fd49856c484eddff2fe64a7d55f4f8ca7260821..3b47dd39041b1dd677e1e54af54aab27b6442367 100644 (file)
 
 #include <deal.II/base/config.h>
 #include <deal.II/base/exceptions.h>
+#include <deal.II/base/std_cxx1x/tuple.h>
+#include <deal.II/base/synchronous_iterator.h>
 #include <deal.II/fe/fe_update_flags.h>
 #include <deal.II/fe/mapping.h>
+#include <deal.II/lac/vector.h>
 #include <utility>
 
 DEAL_II_NAMESPACE_OPEN
@@ -673,12 +676,12 @@ private:
   template <class DerivativeDescription, int dim,
            template <int, int> class DH, class InputVector, int spacedim>
   static void
-  approximate (const typename DH<dim,spacedim>::active_cell_iterator &cell,
+  approximate (SynchronousIterators<std_cxx1x::tuple<typename DH<dim,spacedim>
+               ::active_cell_iterator,Vector<float>::iterator> > const &cell,
                const Mapping<dim,spacedim>    &mapping,
                const DH<dim,spacedim>         &dof,
                const InputVector     &solution,
-               const unsigned int     component,
-               Vector<float>         &derivative_norm);
+               const unsigned int     component);
 
   /**
    * Compute the derivative approximation on
index d4d7a03852b55b0752fd85b0f120c843e42e1ae2..e5518afcb83d345b787794f5604b15b601d83f52 100644 (file)
@@ -6523,9 +6523,8 @@ void DataOutInterface<dim,spacedim>::write_vtu_in_parallel (const char *filename
   MPI_File fh;
   err = MPI_File_open(comm, const_cast<char *>(filename),
                 MPI_MODE_CREATE | MPI_MODE_WRONLY, info, &fh);
-  std::cout << err << std::endl;
   AssertThrow(err==0, ExcMessage("Unable to open file with MPI_File_open!"));
-  
+
 
   MPI_File_set_size(fh, 0); // delete the file contents
   // this barrier is necessary, because otherwise others might already
index df146a150b1dbbfb6bac450c64a5604b707bf1a5..669dccb71ba7d2e12d19b20dae822c151ea39211 100644 (file)
@@ -1115,16 +1115,17 @@ next_cell:
                               +
                               1);
 
-    // next we have to build a mapping from the
-    // list of cells to their indices. for
-    // this, use the user_index field
-    std::vector<unsigned int> saved_user_indices;
-    triangulation.save_user_indices (saved_user_indices);
+    // create a map pair<lvl,idx> -> SparsityPattern index
+    // TODO: we are no longer using user_indices for this because we can get
+    // pointer/index clashes when saving/restoring them. The following approach
+    // works, but this map can get quite big. Not sure about more efficient solutions.
+    std::map< std::pair<unsigned int,unsigned int>, unsigned int >
+      indexmap;
     unsigned int index = 0;
     for (typename Triangulation<dim,spacedim>::active_cell_iterator
          cell = triangulation.begin_active();
          cell != triangulation.end(); ++cell, ++index)
-      cell->set_user_index (index);
+      indexmap[std::pair<unsigned int,unsigned int>(cell->level(),cell->index())] = index;
 
     // next loop over all cells and
     // their neighbors to build the
@@ -1153,21 +1154,15 @@ next_cell:
               &&
               (cell->neighbor(f)->has_children() == false))
             {
-              cell_connectivity.add (index,
-                                     cell->neighbor(f)->user_index());
-              cell_connectivity.add (cell->neighbor(f)->user_index(),
-                                     index);
+              unsigned int other_index = indexmap.find(
+                  std::pair<unsigned int,unsigned int>(cell->neighbor(f)->level(),cell->neighbor(f)->index()))->second;
+              cell_connectivity.add (index, other_index);
+              cell_connectivity.add (other_index, index);
             }
       }
 
-    // now compress the so-built connectivity
-    // pattern and restore user indices. the
-    // const-cast is necessary since we treat
-    // the triangulation as constant (we here
-    // return it to its original state)
+    // now compress the so-built connectivity pattern
     cell_connectivity.compress ();
-    const_cast<Triangulation<dim,spacedim>&>(triangulation)
-    .load_user_indices (saved_user_indices);
   }
 
 
index 03e0cce91a1e16a2624a516831e7c2aedc2ced83..be256797bb887e0798432c0ff17a5f0e405608f3 100644 (file)
@@ -59,17 +59,11 @@ namespace internal
 
 
     /**
-     * In a WorkStream context, use this function to append the patch computed
-     * by the parallel stage to the array of patches.
+     * Dummy function used for WorkStream.
      */
     template <int dim, int spacedim>
     void
-    append_patch_to_list (const DataOutBase::Patch<dim,spacedim> &patch,
-                          std::vector<DataOutBase::Patch<dim,spacedim> > &patches)
-    {
-      patches.push_back (patch);
-      patches.back().patch_index = patches.size()-1;
-    }
+    copy(const DataOutBase::Patch<dim,spacedim> &patch) {}
   }
 }
 
@@ -81,7 +75,8 @@ DataOut<dim,DH>::
 build_one_patch (const std::pair<cell_iterator, unsigned int> *cell_and_index,
                  internal::DataOut::ParallelData<DH::dimension, DH::space_dimension> &data,
                  DataOutBase::Patch<DH::dimension, DH::space_dimension> &patch,
-                 const CurvedCellRegion curved_cell_region)
+                 const CurvedCellRegion curved_cell_region,
+                 std::vector<DataOutBase::Patch<DH::dimension, DH::space_dimension> > &patches)
 {
   // use ucd_to_deal map as patch vertices are in the old, unnatural
   // ordering. if the mapping does not preserve locations
@@ -306,6 +301,14 @@ build_one_patch (const std::pair<cell_iterator, unsigned int> *cell_and_index,
       patch.neighbors[f]
         = (*data.cell_to_patch_index_map)[neighbor->level()][neighbor->index()];
     }
+
+  const unsigned int patch_idx = cell_and_index->second;
+  // did we mess up the indices?
+  Assert(patch_idx < patches.size(), ExcInternalError());
+  
+  // Put the patch in the patches vector
+  patches[patch_idx] = patch;
+  patches[patch_idx].patch_index = patch_idx;
 }
 
 
@@ -357,19 +360,8 @@ void DataOut<dim,DH>::build_patches (const Mapping<DH::dimension,DH::space_dimen
 
   std::vector<std::pair<cell_iterator, unsigned int> > all_cells;
   {
-    // set the index of the first cell. if first_locally_owned_cell /
-    // next_locally_owned_cell returns non-active cells, then the index is not
-    // usable anyway, but otherwise we should keep track where we are
-    unsigned int index;
     cell_iterator cell = first_locally_owned_cell();
-    if ((cell == this->triangulation->end())
-        ||
-        (cell->has_children()))
-      index = 0;
-    else
-      index = std::distance (this->triangulation->begin_active(),
-                             active_cell_iterator(cell));
-    for ( ; cell != this->triangulation->end(); )
+    for (unsigned int index = 0; cell != this->triangulation->end(); ++index)
       {
         Assert (static_cast<unsigned int>(cell->level()) <
                 cell_to_patch_index_map.size(),
@@ -381,25 +373,12 @@ void DataOut<dim,DH>::build_patches (const Mapping<DH::dimension,DH::space_dimen
         cell_to_patch_index_map[cell->level()][cell->index()] = all_cells.size();
 
         all_cells.push_back (std::make_pair(cell, index));
-
-        // if both this and the next cell are active, then increment the index
-        // that keeps track on which active cell we are sitting correctly. if
-        // one of the cells is not active, then this index doesn't mean
-        // anything anyway, so just ignore it. same if we are at the end of
-        // the range
-        cell_iterator next = next_locally_owned_cell(cell);
-        if (!cell->has_children() &&
-            next != this->triangulation->end() &&
-            !next->has_children())
-          index += std::distance (active_cell_iterator(cell),
-                                  active_cell_iterator(next));
-        cell = next;
+        cell = next_locally_owned_cell(cell);
       }
   }
 
   this->patches.clear ();
-  this->patches.reserve (all_cells.size());
-  Assert (this->patches.size() == 0, ExcInternalError());
+  this->patches.resize(all_cells.size());
 
   // now create a default object for the WorkStream object to work with
   unsigned int n_datasets=this->cell_data.size();
@@ -449,9 +428,9 @@ void DataOut<dim,DH>::build_patches (const Mapping<DH::dimension,DH::space_dimen
                      &all_cells[0]+all_cells.size(),
                      std_cxx1x::bind(&DataOut<dim,DH>::build_one_patch,
                                      *this, std_cxx1x::_1, std_cxx1x::_2, std_cxx1x::_3,
-                                     curved_cell_region),
-                     std_cxx1x::bind(&internal::DataOut::append_patch_to_list<dim,DH::space_dimension>,
-                                     std_cxx1x::_1, std_cxx1x::ref(this->patches)),
+                                     curved_cell_region,std_cxx1x::ref(this->patches)),
+                     std_cxx1x::bind(&internal::DataOut::copy<dim,DH::space_dimension>,
+                                     std_cxx1x::_1),
                      thread_data,
                      sample_patch);
 }
@@ -486,7 +465,7 @@ typename DataOut<dim,DH>::cell_iterator
 DataOut<dim,DH>::first_locally_owned_cell ()
 {
   typename DataOut<dim,DH>::cell_iterator
-  cell = this->triangulation->begin_active ();
+    cell = first_cell();
 
   // skip cells if the current one has no children (is active) and is a ghost
   // or artificial cell
index 397e688c60d98f74201d95f84856167b24e4ee8c..cb8825c807ca698c97508bdd92e82adea7420ee9 100644 (file)
@@ -650,20 +650,25 @@ approximate_derivative (const Mapping<dim,spacedim>    &mapping,
 
   // Only act on the locally owned cells
   typedef FilteredIterator<typename DH<dim,spacedim>::active_cell_iterator> CellFilter;
-
-  // There is no need for a copier because there is no conflict between threads
+  
+  typedef std_cxx1x::tuple<typename DH<dim,spacedim>::active_cell_iterator,Vector<float>::iterator>
+    Iterators;
+  SynchronousIterators<Iterators> begin(Iterators (CellFilter(IteratorFilters::LocallyOwnedCell(),
+            dof_handler.begin_active()),derivative_norm.begin())),
+    end(Iterators (CellFilter(IteratorFilters::LocallyOwnedCell(),dof_handler.end()),
+          derivative_norm.end()));
+  
+  // There is no need for a copier because there is no conflict between threads 
   // to write in derivative_norm. Scratch and CopyData are also useless.
-  WorkStream::run(CellFilter(IteratorFilters::LocallyOwnedCell(),dof_handler.begin_active()),
-      CellFilter(IteratorFilters::LocallyOwnedCell(),dof_handler.end()),
-      static_cast<std_cxx1x::function<void (typename DH<dim,spacedim>::active_cell_iterator const&,
+  WorkStream::run(begin,end,
+      static_cast<std_cxx1x::function<void (SynchronousIterators<Iterators> const&,
         internal::Assembler::Scratch const&,internal::Assembler::CopyData &)> >
       (std_cxx1x::bind(DerivativeApproximation::template approximate<DerivativeDescription,dim,DH,
                        InputVector,spacedim>,
                        std_cxx1x::_1,
                        std_cxx1x::cref(mapping),
                        std_cxx1x::cref(dof_handler),
-                       std_cxx1x::cref(solution),component,
-                       std_cxx1x::ref(derivative_norm))),
+                       std_cxx1x::cref(solution),component)),
     static_cast<std_cxx1x::function<void (internal::Assembler::CopyData const &)> >
       (internal::Assembler::copier),internal::Assembler::Scratch (),
       internal::Assembler::CopyData ());
@@ -674,27 +679,21 @@ approximate_derivative (const Mapping<dim,spacedim>    &mapping,
 template <class DerivativeDescription, int dim,
          template <int, int> class DH, class InputVector, int spacedim>
 void
-DerivativeApproximation::approximate (const typename DH<dim,spacedim>::active_cell_iterator &cell,
+DerivativeApproximation::approximate (SynchronousIterators<std_cxx1x::tuple<typename DH<dim,spacedim>
+                                      ::active_cell_iterator,Vector<float>::iterator> > const &cell,
                                       const Mapping<dim,spacedim>                  &mapping,
                                       const DH<dim,spacedim>                       &dof_handler,
                                       const InputVector                            &solution,
-                                      const unsigned int                            component,
-                                      Vector<float>                                &derivative_norm)
+                                      const unsigned int                            component)
 {
-  // iterators over all cells and the
-  // respective entries in the output
-  // vector:
-  Vector<float>::iterator derivative_norm_on_this_cell = derivative_norm.begin() + 
-    std::distance(dof_handler.begin_active(),cell);
-
   typename DerivativeDescription::Derivative derivative;
   // call the function doing the actual
   // work on this cell
   DerivativeApproximation::template approximate_cell<DerivativeDescription,dim,DH,InputVector>
-    (mapping,dof_handler,solution,component,cell,derivative);
+    (mapping,dof_handler,solution,component,std_cxx1x::get<0>(cell.iterators),derivative);
   // evaluate the norm and fill the vector
-  *derivative_norm_on_this_cell
-    = DerivativeDescription::derivative_norm (derivative);
+  //*derivative_norm_on_this_cell
+  *std_cxx1x::get<1>(cell.iterators) = DerivativeDescription::derivative_norm (derivative);
 }
 
 
index af0e940cc8bbdbc591900945156a28f28d586cb4..75d21fe9e83926acdc8c7e2ea5bb6a634dccaeef 100644 (file)
@@ -16,7 +16,6 @@
 
 #include <deal.II/base/function.h>
 #include <deal.II/base/quadrature.h>
-#include <deal.II/base/thread_management.h>
 #include <deal.II/base/work_stream.h>
 #include <deal.II/base/multithread_info.h>
 #include <deal.II/base/geometry_info.h>
@@ -603,6 +602,41 @@ namespace MatrixCreator
                                                      data.dof_indices,
                                                      *matrix);
     }
+
+
+
+    namespace AssemblerBoundary
+    {
+      struct Scratch
+      {
+        Scratch() {}
+      };
+
+      template <typename DH>
+      struct CopyData
+      {
+        CopyData() {};
+
+        CopyData(CopyData const &data);
+
+        unsigned int dofs_per_cell;
+        std::vector<types::global_dof_index> dofs;
+        std::vector<std::vector<bool> > dof_is_on_face;
+        typename DH::active_cell_iterator cell;
+        std::vector<FullMatrix<double> > cell_matrix;
+        std::vector<Vector<double> > cell_vector;
+      };
+
+      template <typename DH>
+      CopyData<DH>::CopyData(CopyData const &data) :
+       dofs_per_cell(data.dofs_per_cell),
+       dofs(data.dofs),
+       dof_is_on_face(data.dof_is_on_face),
+       cell(data.cell),
+       cell_matrix(data.cell_matrix),
+       cell_vector(data.cell_vector)
+      {}
+    }
   }
 }
 
@@ -834,37 +868,30 @@ namespace MatrixCreator
   {
     template <int dim, int spacedim>
     void
-    create_boundary_mass_matrix_1 (std_cxx1x::tuple<const Mapping<dim, spacedim> &,
-                                   const DoFHandler<dim,spacedim> &,
-                                   const Quadrature<dim-1> & >  commons,
-                                   SparseMatrix<double>      &matrix,
-                                   const typename FunctionMap<spacedim>::type &boundary_functions,
-                                   Vector<double>            &rhs_vector,
-                                   std::vector<types::global_dof_index> &dof_to_boundary_mapping,
-                                   const Function<spacedim> *const coefficient,
-                                   const std::vector<unsigned int> &component_mapping,
-                                   const MatrixCreator::internal::IteratorRange<DoFHandler<dim,spacedim> >   range,
-                                   Threads::Mutex      &mutex)
+    create_boundary_mass_matrix_1 (typename DoFHandler<dim,spacedim>::active_cell_iterator const &cell,
+                                   MatrixCreator::internal::AssemblerBoundary::Scratch const &scratch,
+                                   MatrixCreator::internal::AssemblerBoundary::CopyData<DoFHandler<dim,
+                                   spacedim> > &copy_data,
+                                   Mapping<dim, spacedim> const &mapping,
+                                   FiniteElement<dim,spacedim> const &fe,
+                                   Quadrature<dim-1> const &q,
+                                   typename FunctionMap<spacedim>::type const &boundary_functions,
+                                   Function<spacedim> const *const coefficient,
+                                   std::vector<unsigned int> const &component_mapping)
+                                   
     {
       // All assertions for this function
       // are in the calling function
       // before creating threads.
-      const Mapping<dim,spacedim> &mapping = std_cxx1x::get<0>(commons);
-      const DoFHandler<dim,spacedim> &dof = std_cxx1x::get<1>(commons);
-      const Quadrature<dim-1>& q = std_cxx1x::get<2>(commons);
-
-      const FiniteElement<dim,spacedim> &fe = dof.get_fe();
-      const unsigned int n_components  = fe.n_components();
+      const unsigned int n_components = fe.n_components();
       const unsigned int n_function_components = boundary_functions.begin()->second->n_components;
-      const bool         fe_is_system  = (n_components != 1);
+      const bool         fe_is_system = (n_components != 1);
       const bool         fe_is_primitive = fe.is_primitive();
 
-      const unsigned int dofs_per_cell = fe.dofs_per_cell,
-                         dofs_per_face = fe.dofs_per_face;
-
-      FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
-      Vector<double>     cell_vector(dofs_per_cell);
-
+      const unsigned int dofs_per_face = fe.dofs_per_face;
+      
+      copy_data.cell = cell;
+      copy_data.dofs_per_cell = fe.dofs_per_cell;
 
       UpdateFlags update_flags = UpdateFlags (update_values     |
                                               update_JxW_values |
@@ -883,236 +910,246 @@ namespace MatrixCreator
 
       std::vector<double>          rhs_values_scalar (fe_values.n_quadrature_points);
       std::vector<Vector<double> > rhs_values_system (fe_values.n_quadrature_points,
-                                                      Vector<double>(n_function_components));
-
-      std::vector<types::global_dof_index> dofs (dofs_per_cell);
-      std::vector<types::global_dof_index> dofs_on_face_vector (dofs_per_face);
+          Vector<double>(n_function_components));
 
-      // for each dof on the cell, have a
-      // flag whether it is on the face
-      std::vector<bool>         dof_is_on_face(dofs_per_cell);
-
-      typename DoFHandler<dim,spacedim>::active_cell_iterator cell = range.first;
-      for (; cell!=range.second; ++cell)
-        for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
-          // check if this face is on that part of
-          // the boundary we are interested in
-          if (boundary_functions.find(cell->face(face)->boundary_indicator()) !=
-              boundary_functions.end())
-            {
-              cell_matrix = 0;
-              cell_vector = 0;
+      copy_data.dofs.resize(copy_data.dofs_per_cell);
+      cell->get_dof_indices (copy_data.dofs);
 
-              fe_values.reinit (cell, face);
+      std::vector<types::global_dof_index> dofs_on_face_vector (dofs_per_face);
 
-              if (fe_is_system)
-                // FE has several components
-                {
-                  boundary_functions.find(cell->face(face)->boundary_indicator())
-                  ->second->vector_value_list (fe_values.get_quadrature_points(),
-                                               rhs_values_system);
-
-                  if (coefficient_is_vector)
-                    // If coefficient is
-                    // vector valued, fill
-                    // all components
-                    coefficient->vector_value_list (fe_values.get_quadrature_points(),
-                                                    coefficient_vector_values);
-                  else
-                    {
-                      // If a scalar
-                      // function is
-                      // geiven, update
-                      // the values, if
-                      // not, use the
-                      // default one set
-                      // in the
-                      // constructor above
-                      if (coefficient != 0)
-                        coefficient->value_list (fe_values.get_quadrature_points(),
-                                                 coefficient_values);
-                      // Copy scalar
-                      // values into vector
-                      for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
-                        coefficient_vector_values[point] = coefficient_values[point];
-                    }
+      // Because CopyData objects are reused and that push_back is used,
+      // dof_is_on_face, cell_matrix, and cell_vector must be cleared before
+      // they are reused
+      copy_data.dof_is_on_face.clear();
+      copy_data.cell_matrix.clear();
+      copy_data.cell_vector.clear();
+
+      for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+        // check if this face is on that part of
+        // the boundary we are interested in
+        if (boundary_functions.find(cell->face(face)->boundary_indicator()) !=
+            boundary_functions.end())
+          {
+            copy_data.cell_matrix.push_back(FullMatrix<double> (copy_data.dofs_per_cell,
+                  copy_data.dofs_per_cell));
+            copy_data.cell_vector.push_back(Vector<double> (copy_data.dofs_per_cell));
+            fe_values.reinit (cell, face);
 
-                  // Special treatment
-                  // for Hdiv and Hcurl
-                  // elements, where only
-                  // the normal or
-                  // tangential component
-                  // should be projected.
-                  std::vector<std::vector<double> > normal_adjustment(fe_values.n_quadrature_points,
-                                                                      std::vector<double>(n_components, 1.));
+            if (fe_is_system)
+              // FE has several components
+              {
+                boundary_functions.find(cell->face(face)->boundary_indicator())
+                ->second->vector_value_list (fe_values.get_quadrature_points(),
+                                             rhs_values_system);
+
+                if (coefficient_is_vector)
+                  // If coefficient is
+                  // vector valued, fill
+                  // all components
+                  coefficient->vector_value_list (fe_values.get_quadrature_points(),
+                                                  coefficient_vector_values);
+                else
+                  {
+                    // If a scalar
+                    // function is
+                    // given, update
+                    // the values, if
+                    // not, use the
+                    // default one set
+                    // in the
+                    // constructor above
+                    if (coefficient != 0)
+                      coefficient->value_list (fe_values.get_quadrature_points(),
+                                               coefficient_values);
+                    // Copy scalar
+                    // values into vector
+                    for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+                      coefficient_vector_values[point] = coefficient_values[point];
+                  }
 
-                  for (unsigned int comp = 0; comp<n_components; ++comp)
-                    {
-                      const FiniteElement<dim,spacedim> &base = fe.base_element(fe.component_to_base_index(comp).first);
-                      const unsigned int bcomp = fe.component_to_base_index(comp).second;
+                // Special treatment
+                // for Hdiv and Hcurl
+                // elements, where only
+                // the normal or
+                // tangential component
+                // should be projected.
+                std::vector<std::vector<double> > normal_adjustment(fe_values.n_quadrature_points,
+                                                                    std::vector<double>(n_components, 1.));
 
-                      if (!base.conforms(FiniteElementData<dim>::H1) &&
-                          base.conforms(FiniteElementData<dim>::Hdiv))
-                        for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
-                          normal_adjustment[point][comp] = fe_values.normal_vector(point)(bcomp)
-                                                           * fe_values.normal_vector(point)(bcomp);
-                    }
+                for (unsigned int comp = 0; comp<n_components; ++comp)
+                  {
+                    const FiniteElement<dim,spacedim> &base = fe.base_element(fe.component_to_base_index(comp).first);
+                    const unsigned int bcomp = fe.component_to_base_index(comp).second;
 
-                  for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
-                    {
-                      const double weight = fe_values.JxW(point);
-                      for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
-                        if (fe_is_primitive)
-                          {
-                            for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
-                              {
-                                if (fe.system_to_component_index(j).first
-                                    == fe.system_to_component_index(i).first)
-                                  {
-                                    cell_matrix(i,j)
-                                    += weight
-                                       * fe_values.shape_value(j,point)
-                                       * fe_values.shape_value(i,point)
-                                       * coefficient_vector_values[point](fe.system_to_component_index(i).first);
-                                  }
-                              }
-                            cell_vector(i) += fe_values.shape_value(i,point)
-                                              * rhs_values_system[point](component_mapping[fe.system_to_component_index(i).first])
-                                              * weight;
-                          }
-                        else
-                          {
-                            for (unsigned int comp=0; comp<n_components; ++comp)
-                              {
-                                for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
-                                  cell_matrix(i,j)
-                                  += fe_values.shape_value_component(j,point,comp)
-                                     * fe_values.shape_value_component(i,point,comp)
-                                     * normal_adjustment[point][comp]
-                                     * weight * coefficient_vector_values[point](comp);
-                                cell_vector(i) += fe_values.shape_value_component(i,point,comp) *
-                                                  rhs_values_system[point](component_mapping[comp])
-                                                  * normal_adjustment[point][comp]
-                                                  * weight;
-                              }
-                          }
-                    }
-                }
-              else
-                // FE is a scalar one
-                {
-                  boundary_functions.find(cell->face(face)->boundary_indicator())
-                  ->second->value_list (fe_values.get_quadrature_points(), rhs_values_scalar);
+                    if (!base.conforms(FiniteElementData<dim>::H1) &&
+                        base.conforms(FiniteElementData<dim>::Hdiv))
+                      for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+                        normal_adjustment[point][comp] = fe_values.normal_vector(point)(bcomp)
+                                                         * fe_values.normal_vector(point)(bcomp);
+                  }
 
-                  if (coefficient != 0)
-                    coefficient->value_list (fe_values.get_quadrature_points(),
-                                             coefficient_values);
-                  for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
-                    {
-                      const double weight = fe_values.JxW(point);
-                      for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
+                for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+                  {
+                    const double weight = fe_values.JxW(point);
+                    for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
+                      if (fe_is_primitive)
                         {
-                          const double v = fe_values.shape_value(i,point);
                           for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
                             {
-                              const double u = fe_values.shape_value(j,point);
-                              cell_matrix(i,j) += (u * v * weight * coefficient_values[point]);
+                              if (fe.system_to_component_index(j).first
+                                  == fe.system_to_component_index(i).first)
+                                {
+                                  copy_data.cell_matrix.back()(i,j)
+                                  += weight
+                                     * fe_values.shape_value(j,point)
+                                     * fe_values.shape_value(i,point)
+                                     * coefficient_vector_values[point](fe.system_to_component_index(i).first);
+                                }
                             }
-                          cell_vector(i) += v * rhs_values_scalar[point] *weight;
+                          copy_data.cell_vector.back()(i) += fe_values.shape_value(i,point)
+                                            * rhs_values_system[point](component_mapping[fe.system_to_component_index(i).first])
+                                            * weight;
                         }
-                    }
-                }
-              // now transfer cell matrix and vector
-              // to the whole boundary matrix
-              //
-              // in the following: dof[i] holds the
-              // global index of the i-th degree of
-              // freedom on the present cell. If it
-              // is also a dof on the boundary, it
-              // must be a nonzero entry in the
-              // dof_to_boundary_mapping and then
-              // the boundary index of this dof is
-              // dof_to_boundary_mapping[dof[i]].
-              //
-              // if dof[i] is not on the boundary,
-              // it should be zero on the boundary
-              // therefore on all quadrature
-              // points and finally all of its
-              // entries in the cell matrix and
-              // vector should be zero. If not, we
-              // throw an error (note: because of
-              // the evaluation of the shape
-              // functions only up to machine
-              // precision, the term "must be zero"
-              // really should mean: "should be
-              // very small". since this is only an
-              // assertion and not part of the
-              // code, we may choose "very small"
-              // quite arbitrarily)
-              //
-              // the main problem here is that the
-              // matrix or vector entry should also
-              // be zero if the degree of freedom
-              // dof[i] is on the boundary, but not
-              // on the present face, i.e. on
-              // another face of the same cell also
-              // on the boundary. We can therefore
-              // not rely on the
-              // dof_to_boundary_mapping[dof[i]]
-              // being !=-1, we really have to
-              // determine whether dof[i] is a
-              // dof on the present face. We do so
-              // by getting the dofs on the
-              // face into @p{dofs_on_face_vector},
-              // a vector as always. Usually,
-              // searching in a vector is
-              // inefficient, so we copy the dofs
-              // into a set, which enables binary
-              // searches.
-              cell->get_dof_indices (dofs);
-              cell->face(face)->get_dof_indices (dofs_on_face_vector);
-              for (unsigned int i=0; i<dofs_per_cell; ++i)
-                dof_is_on_face[i] = (std::find(dofs_on_face_vector.begin(),
-                                               dofs_on_face_vector.end(),
-                                               dofs[i])
-                                     !=
-                                     dofs_on_face_vector.end());
-
-              // lock the matrix
-              Threads::Mutex::ScopedLock lock (mutex);
-              for (unsigned int i=0; i<dofs_per_cell; ++i)
-                {
-                  if (dof_is_on_face[i] && dof_to_boundary_mapping[dofs[i]] != numbers::invalid_dof_index)
-                    {
-                      for (unsigned int j=0; j<dofs_per_cell; ++j)
-                        if (dof_is_on_face[j] && dof_to_boundary_mapping[dofs[j]] != numbers::invalid_dof_index)
+                      else
+                        {
+                          for (unsigned int comp=0; comp<n_components; ++comp)
+                            {
+                              for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
+                                copy_data.cell_matrix.back()(i,j)
+                                += fe_values.shape_value_component(j,point,comp)
+                                   * fe_values.shape_value_component(i,point,comp)
+                                   * normal_adjustment[point][comp]
+                                   * weight * coefficient_vector_values[point](comp);
+                              copy_data.cell_vector.back()(i) += fe_values.shape_value_component(i,point,comp) *
+                                                rhs_values_system[point](component_mapping[comp])
+                                                * normal_adjustment[point][comp]
+                                                * weight;
+                            }
+                        }
+                  }
+              }
+            else
+              // FE is a scalar one
+              {
+                boundary_functions.find(cell->face(face)->boundary_indicator())
+                ->second->value_list (fe_values.get_quadrature_points(), rhs_values_scalar);
+
+                if (coefficient != 0)
+                  coefficient->value_list (fe_values.get_quadrature_points(),
+                                           coefficient_values);
+                for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+                  {
+                    const double weight = fe_values.JxW(point);
+                    for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
+                      {
+                        const double v = fe_values.shape_value(i,point);
+                        for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
                           {
-                            Assert(numbers::is_finite(cell_matrix(i,j)), ExcNumberNotFinite());
-                            matrix.add(dof_to_boundary_mapping[dofs[i]],
-                                       dof_to_boundary_mapping[dofs[j]],
-                                       cell_matrix(i,j));
+                            const double u = fe_values.shape_value(j,point);
+                            copy_data.cell_matrix.back()(i,j) += (u*v*weight*coefficient_values[point]);
                           }
-                      Assert(numbers::is_finite(cell_vector(i)), ExcNumberNotFinite());
-                      rhs_vector(dof_to_boundary_mapping[dofs[i]]) += cell_vector(i);
-                    }
-                }
-            }
+                        copy_data.cell_vector.back()(i) += v * rhs_values_scalar[point] *weight;
+                      }
+                  }
+              }
+
+
+            cell->face(face)->get_dof_indices (dofs_on_face_vector);
+            // for each dof on the cell, have a
+            // flag whether it is on the face
+            copy_data.dof_is_on_face.push_back(std::vector<bool> (copy_data.dofs_per_cell));
+            // check for each of the dofs on this cell
+            // whether it is on the face
+            for (unsigned int i=0; i<copy_data.dofs_per_cell; ++i)
+              copy_data.dof_is_on_face.back()[i] = (std::find(dofs_on_face_vector.begin(),
+                                                              dofs_on_face_vector.end(),
+                                                              copy_data.dofs[i])
+                                                    !=
+                                                    dofs_on_face_vector.end());
+          }   
+    }
+
+    template <int dim,int spacedim>
+    void copy_boundary_mass_matrix_1(MatrixCreator::internal::AssemblerBoundary::CopyData<DoFHandler<dim,
+                                     spacedim> > const &copy_data,
+                                     typename FunctionMap<spacedim>::type const &boundary_functions,
+                                     std::vector<types::global_dof_index> const &dof_to_boundary_mapping,
+                                     SparseMatrix<double> &matrix,
+                                     Vector<double> &rhs_vector)
+    {
+      // now transfer cell matrix and vector to the whole boundary matrix
+      //
+      // in the following: dof[i] holds the global index of the i-th degree of
+      // freedom on the present cell. If it is also a dof on the boundary, it
+      // must be a nonzero entry in the dof_to_boundary_mapping and then
+      // the boundary index of this dof is dof_to_boundary_mapping[dof[i]].
+      //
+      // if dof[i] is not on the boundary, it should be zero on the boundary
+      // therefore on all quadrature points and finally all of its
+      // entries in the cell matrix and vector should be zero. If not, we
+      // throw an error (note: because of the evaluation of the shape
+      // functions only up to machine precision, the term "must be zero"
+      // really should mean: "should be very small". since this is only an
+      // assertion and not part of the code, we may choose "very small"
+      // quite arbitrarily)
+      //
+      // the main problem here is that the matrix or vector entry should also
+      // be zero if the degree of freedom dof[i] is on the boundary, but not
+      // on the present face, i.e. on another face of the same cell also
+      // on the boundary. We can therefore not rely on the
+      // dof_to_boundary_mapping[dof[i]] being !=-1, we really have to
+      // determine whether dof[i] is a dof on the present face. We do so
+      // by getting the dofs on the face into @p{dofs_on_face_vector},
+      // a vector as always. Usually, searching in a vector is
+      // inefficient, so we copy the dofs into a set, which enables binary
+      // searches.
+      unsigned int pos(0);
+      for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+      {
+        // check if this face is on that part of
+        // the boundary we are interested in
+        if (boundary_functions.find(copy_data.cell->face(face)->boundary_indicator()) !=
+            boundary_functions.end())
+          {
+            for (unsigned int i=0; i<copy_data.dofs_per_cell; ++i)
+              {
+                if (copy_data.dof_is_on_face[pos][i] && 
+                    dof_to_boundary_mapping[copy_data.dofs[i]] != numbers::invalid_dof_index)
+                  {
+                    for (unsigned int j=0; j<copy_data.dofs_per_cell; ++j)
+                      if (copy_data.dof_is_on_face[pos][j] && 
+                          dof_to_boundary_mapping[copy_data.dofs[j]] != numbers::invalid_dof_index)
+                        {
+                          Assert(numbers::is_finite(copy_data.cell_matrix[pos](i,j)),
+                              ExcNumberNotFinite());
+                          matrix.add(dof_to_boundary_mapping[copy_data.dofs[i]],
+                                     dof_to_boundary_mapping[copy_data.dofs[j]],
+                                     copy_data.cell_matrix[pos](i,j));
+                        }
+                    Assert(numbers::is_finite(copy_data.cell_vector[pos](i)), ExcNumberNotFinite());
+                    rhs_vector(dof_to_boundary_mapping[copy_data.dofs[i]]) += copy_data.cell_vector[pos](i);
+                  }
+              }
+            ++pos;
+          }
+      }
     }
 
 
+
     template <>
     void
-    create_boundary_mass_matrix_1<2,3> (std_cxx1x::tuple<const Mapping<2,3> &,
-                                        const DoFHandler<2,3> &,
-                                        const Quadrature<1> & > ,
-                                        SparseMatrix<double> &,
-                                        const FunctionMap<3>::type &,
-                                        Vector<double> &,
-                                        std::vector<types::global_dof_index> &,
-                                        const Function<3> *const ,
-                                        const std::vector<unsigned int> &,
-                                        const MatrixCreator::internal::IteratorRange<DoFHandler<2,3> > ,
-                                        Threads::Mutex &)
+    create_boundary_mass_matrix_1<2,3> (DoFHandler<2,3>::active_cell_iterator const &cell,
+                                        MatrixCreator::internal::AssemblerBoundary::Scratch const 
+                                        &scratch,                                                 
+                                        MatrixCreator::internal::AssemblerBoundary::CopyData<DoFHandler<2,
+                                        3> > &copy_data,                                               
+                                        Mapping<2,3> const &mapping,                              
+                                        FiniteElement<2,3> const &fe,                             
+                                        Quadrature<1> const &q,                                   
+                                        FunctionMap<3>::type const &boundary_functions,  
+                                        Function<3> const *const coefficient,                     
+                                        std::vector<unsigned int> const &component_mapping)       
     {
       Assert(false,ExcNotImplemented());
     }
@@ -1121,17 +1158,17 @@ namespace MatrixCreator
 
     template <>
     void
-    create_boundary_mass_matrix_1<1,3> (std_cxx1x::tuple<const Mapping<1,3> &,
-                                        const DoFHandler<1,3> &,
-                                        const Quadrature<0> & > ,
-                                        SparseMatrix<double> &,
-                                        const FunctionMap<3>::type &,
-                                        Vector<double> &,
-                                        std::vector<types::global_dof_index> &,
-                                        const Function<3> *const ,
-                                        const std::vector<unsigned int> &,
-                                        const MatrixCreator::internal::IteratorRange<DoFHandler<1,3> > ,
-                                        Threads::Mutex &)
+    create_boundary_mass_matrix_1<1,3> (DoFHandler<1,3>::active_cell_iterator const &cell,
+                                        MatrixCreator::internal::AssemblerBoundary::Scratch const 
+                                        &scratch,
+                                        MatrixCreator::internal::AssemblerBoundary::CopyData<DoFHandler<1,
+                                        3> > &copy_data,                             
+                                        Mapping<1,3> const &mapping,
+                                        FiniteElement<1,3> const &fe,
+                                        Quadrature<0> const &q,
+                                        FunctionMap<3>::type const &boundary_functions,
+                                        Function<3> const *const coefficient,
+                                        std::vector<unsigned int> const &component_mapping)                
     {
       Assert(false,ExcNotImplemented());
     }
@@ -1184,47 +1221,24 @@ namespace MatrixCreator
     else
       AssertDimension (n_components, component_mapping.size());
 
-    const unsigned int n_threads = multithread_info.n_threads();
-    Threads::ThreadGroup<> threads;
-
-    // define starting and end point
-    // for each thread
-    typedef typename DoFHandler<dim,spacedim>::active_cell_iterator active_cell_iterator;
-    std::vector<std::pair<active_cell_iterator,active_cell_iterator> > thread_ranges
-      = Threads::split_range<active_cell_iterator> (dof.begin_active(),
-                                                    dof.end(), n_threads);
-
-    // mutex to synchronise access to
-    // the matrix
-    Threads::Mutex mutex;
-
-    typedef std_cxx1x::tuple<const Mapping<dim,spacedim> &,
-            const DoFHandler<dim,spacedim> &,
-            const Quadrature<dim-1>&> Commons;
-
-    // then assemble in parallel
-    typedef void (*create_boundary_mass_matrix_1_t)
-    (Commons,
-     SparseMatrix<double>      &matrix,
-     const typename FunctionMap<spacedim>::type &boundary_functions,
-     Vector<double>            &rhs_vector,
-     std::vector<types::global_dof_index> &dof_to_boundary_mapping,
-     const Function<spacedim> *const coefficient,
-     const std::vector<unsigned int> &component_mapping,
-     const MatrixCreator::internal::IteratorRange<DoFHandler<dim,spacedim> >   range,
-     Threads::Mutex      &mutex);
-    create_boundary_mass_matrix_1_t p
-      = &create_boundary_mass_matrix_1<dim,spacedim>;
-
-//TODO: Use WorkStream here
-    for (unsigned int thread=0; thread<n_threads; ++thread)
-      threads += Threads::new_thread (p,
-                                      Commons(mapping, dof, q), matrix,
-                                      boundary_functions, rhs_vector,
-                                      dof_to_boundary_mapping, coefficient,
-                                      component_mapping,
-                                      thread_ranges[thread], mutex);
-    threads.join_all ();
+    MatrixCreator::internal::AssemblerBoundary::Scratch scratch;
+    MatrixCreator::internal::AssemblerBoundary::CopyData<DoFHandler<dim,spacedim> > copy_data;
+
+    WorkStream::run(dof.begin_active(),dof.end(),
+        static_cast<std_cxx1x::function<void (typename DoFHandler<dim,spacedim>::active_cell_iterator
+          const &,MatrixCreator::internal::AssemblerBoundary::Scratch const &,
+          MatrixCreator::internal::AssemblerBoundary::CopyData<DoFHandler<dim,spacedim> > &)> > 
+        (std_cxx1x::bind(create_boundary_mass_matrix_1<dim,spacedim>,std_cxx1x::_1,std_cxx1x::_2,
+                         std_cxx1x::_3,
+                         std_cxx1x::cref(mapping),std_cxx1x::cref(fe),std_cxx1x::cref(q),
+                         std_cxx1x::cref(boundary_functions),coefficient,
+                         std_cxx1x::cref(component_mapping))), 
+        static_cast<std_cxx1x::function<void (MatrixCreator::internal::AssemblerBoundary
+          ::CopyData<DoFHandler<dim,spacedim> > const &)> > (std_cxx1x::bind(
+              copy_boundary_mass_matrix_1<dim,spacedim>,std_cxx1x::_1,
+            std_cxx1x::cref(boundary_functions),std_cxx1x::cref(dof_to_boundary_mapping),
+            std_cxx1x::ref(matrix),std_cxx1x::ref(rhs_vector))),
+        scratch,copy_data);
   }
 
 
@@ -1234,37 +1248,28 @@ namespace MatrixCreator
 
     template <int dim, int spacedim>
     void
-    create_boundary_mass_matrix_1 (std_cxx1x::tuple<const hp::MappingCollection<dim,spacedim> &,
-                                   const hp::DoFHandler<dim,spacedim> &,
-                                   const hp::QCollection<dim-1> &> commons,
-                                   SparseMatrix<double>      &matrix,
-                                   const typename FunctionMap<spacedim>::type &boundary_functions,
-                                   Vector<double>            &rhs_vector,
-                                   std::vector<types::global_dof_index> &dof_to_boundary_mapping,
-                                   const Function<spacedim> *const coefficient,
-                                   const std::vector<unsigned int> &component_mapping,
-                                   const MatrixCreator::internal::IteratorRange<hp::DoFHandler<dim,spacedim> >   range,
-                                   Threads::Mutex      &mutex)
+    create_hp_boundary_mass_matrix_1 (typename hp::DoFHandler<dim,spacedim>::active_cell_iterator const 
+                                      &cell,
+                                      MatrixCreator::internal::AssemblerBoundary::Scratch const &scratch,
+                                      MatrixCreator::internal::AssemblerBoundary
+                                      ::CopyData<hp::DoFHandler<dim,spacedim> > &copy_data,
+                                      hp::MappingCollection<dim,spacedim> const &mapping,
+                                      hp::FECollection<dim,spacedim> const &fe_collection,
+                                      hp::QCollection<dim-1> const &q,
+                                      const typename FunctionMap<spacedim>::type &boundary_functions,
+                                      Function<spacedim> const *const coefficient,
+                                      std::vector<unsigned int> const &component_mapping)
     {
-      const hp::MappingCollection<dim,spacedim> &mapping = std_cxx1x::get<0>(commons);
-      const hp::DoFHandler<dim,spacedim> &dof = std_cxx1x::get<1>(commons);
-      const hp::QCollection<dim-1>& q = std_cxx1x::get<2>(commons);
-      const hp::FECollection<dim,spacedim> &fe_collection = dof.get_fe();
       const unsigned int n_components  = fe_collection.n_components();
       const unsigned int n_function_components = boundary_functions.begin()->second->n_components;
       const bool         fe_is_system  = (n_components != 1);
-#ifdef DEBUG
-      if (true)
-        {
-          types::global_dof_index max_element = static_cast<types::global_dof_index>(0);
-          for (std::vector<types::global_dof_index>::const_iterator i=dof_to_boundary_mapping.begin();
-               i!=dof_to_boundary_mapping.end(); ++i)
-            if ((*i != hp::DoFHandler<dim,spacedim>::invalid_dof_index) &&
-                (*i > max_element))
-              max_element = *i;
-          Assert (max_element  == matrix.n()-1, ExcInternalError());
-        };
-#endif
+      const FiniteElement<dim,spacedim> &fe = cell->get_fe();
+      const unsigned int dofs_per_face = fe.dofs_per_face;
+
+      copy_data.cell = cell;
+      copy_data.dofs_per_cell = fe.dofs_per_cell;
+      copy_data.dofs.resize(copy_data.dofs_per_cell);
+      cell->get_dof_indices (copy_data.dofs);
 
       const unsigned int max_dofs_per_cell = fe_collection.max_dofs_per_cell(),
                          max_dofs_per_face = fe_collection.max_dofs_per_face();
@@ -1287,148 +1292,131 @@ namespace MatrixCreator
       std::vector<double>          rhs_values_scalar;
       std::vector<Vector<double> > rhs_values_system;
 
-      std::vector<types::global_dof_index> dofs (max_dofs_per_cell);
-      std::vector<types::global_dof_index> dofs_on_face_vector (max_dofs_per_face);
+      std::vector<types::global_dof_index> dofs_on_face_vector (dofs_per_face);
 
-      // for each dof on the cell, have a
-      // flag whether it is on the face
-      std::vector<bool>         dof_is_on_face(max_dofs_per_cell);
+      copy_data.dofs.resize(copy_data.dofs_per_cell);
+      cell->get_dof_indices (copy_data.dofs);
 
+      // Because CopyData objects are reused and that push_back is used,
+      // dof_is_on_face, cell_matrix, and cell_vector must be cleared before
+      // they are reused
+      copy_data.dof_is_on_face.clear();
+      copy_data.cell_matrix.clear();
+      copy_data.cell_vector.clear();
 
-      typename hp::DoFHandler<dim,spacedim>::active_cell_iterator cell = range.first;
-      for (; cell!=range.second; ++cell)
-        for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
-          // check if this face is on that part of
-          // the boundary we are interested in
-          if (boundary_functions.find(cell->face(face)->boundary_indicator()) !=
-              boundary_functions.end())
-            {
-              x_fe_values.reinit (cell, face);
 
-              const FEFaceValues<dim,spacedim> &fe_values = x_fe_values.get_present_fe_values ();
+      for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+        // check if this face is on that part of
+        // the boundary we are interested in
+        if (boundary_functions.find(cell->face(face)->boundary_indicator()) !=
+            boundary_functions.end())
+          {
+            x_fe_values.reinit (cell, face);
 
-              const FiniteElement<dim,spacedim> &fe = cell->get_fe();
-              const unsigned int dofs_per_cell = fe.dofs_per_cell;
-              const unsigned int dofs_per_face = fe.dofs_per_face;
+            const FEFaceValues<dim,spacedim> &fe_values = x_fe_values.get_present_fe_values ();
 
-              cell_matrix.reinit (dofs_per_cell, dofs_per_cell);
-              cell_vector.reinit (dofs_per_cell);
-              cell_matrix = 0;
-              cell_vector = 0;
+            copy_data.cell_matrix.push_back(FullMatrix<double> (copy_data.dofs_per_cell,
+                  copy_data.dofs_per_cell));
+            copy_data.cell_vector.push_back(Vector<double> (copy_data.dofs_per_cell));
 
-              if (fe_is_system)
-                // FE has several components
-                {
-                  rhs_values_system.resize (fe_values.n_quadrature_points,
-                                            Vector<double>(n_function_components));
-                  boundary_functions.find(cell->face(face)->boundary_indicator())
-                  ->second->vector_value_list (fe_values.get_quadrature_points(),
-                                               rhs_values_system);
+            if (fe_is_system)
+              // FE has several components
+              {
+                rhs_values_system.resize (fe_values.n_quadrature_points,
+                                          Vector<double>(n_function_components));
+                boundary_functions.find(cell->face(face)->boundary_indicator())
+                ->second->vector_value_list (fe_values.get_quadrature_points(),
+                                             rhs_values_system);
 
-                  if (coefficient != 0)
-                    {
-                      if (coefficient->n_components==1)
-                        {
-                          coefficient_values.resize (fe_values.n_quadrature_points);
-                          coefficient->value_list (fe_values.get_quadrature_points(),
-                                                   coefficient_values);
-                          for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
-                            {
-                              const double weight = fe_values.JxW(point);
-                              for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
-                                {
-                                  const double v = fe_values.shape_value(i,point);
-                                  for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
-                                    if (fe.system_to_component_index(i).first ==
-                                        fe.system_to_component_index(j).first)
-                                      {
-                                        const double u = fe_values.shape_value(j,point);
-                                        cell_matrix(i,j)
-                                        += (u * v * weight * coefficient_values[point]);
-                                      }
-
-                                  cell_vector(i) += v *
-                                                    rhs_values_system[point](
-                                                      component_mapping[fe.system_to_component_index(i).first]) * weight;
-                                }
-                            }
-                        }
-                      else
-                        {
-                          coefficient_vector_values.resize (fe_values.n_quadrature_points,
-                                                            Vector<double>(n_components));
-                          coefficient->vector_value_list (fe_values.get_quadrature_points(),
-                                                          coefficient_vector_values);
-                          for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
-                            {
-                              const double weight = fe_values.JxW(point);
-                              for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
-                                {
-                                  const double v = fe_values.shape_value(i,point);
-                                  const unsigned int component_i=
-                                    fe.system_to_component_index(i).first;
-                                  for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
-                                    if (fe.system_to_component_index(j).first ==
-                                        component_i)
-                                      {
-                                        const double u = fe_values.shape_value(j,point);
-                                        cell_matrix(i,j) +=
-                                          (u * v * weight * coefficient_vector_values[point](component_i));
-                                      }
-                                  cell_vector(i) += v * rhs_values_system[point](component_mapping[component_i]) * weight;
-                                }
-                            }
-                        }
-                    }
-                  else  //      if (coefficient == 0)
-                    for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+                if (coefficient != 0)
+                  {
+                    if (coefficient->n_components==1)
                       {
-                        const double weight = fe_values.JxW(point);
-                        for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
+                        coefficient_values.resize (fe_values.n_quadrature_points);
+                        coefficient->value_list (fe_values.get_quadrature_points(),
+                                                 coefficient_values);
+                        for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
                           {
-                            const double v = fe_values.shape_value(i,point);
-                            for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
-                              if (fe.system_to_component_index(i).first ==
-                                  fe.system_to_component_index(j).first)
-                                {
-                                  const double u = fe_values.shape_value(j,point);
-                                  cell_matrix(i,j) += (u * v * weight);
-                                }
-                            cell_vector(i) += v *
-                                              rhs_values_system[point](
-                                                fe.system_to_component_index(i).first) *
-                                              weight;
+                            const double weight = fe_values.JxW(point);
+                            for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
+                              {
+                                const double v = fe_values.shape_value(i,point);
+                                for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
+                                  if (fe.system_to_component_index(i).first ==
+                                      fe.system_to_component_index(j).first)
+                                    {
+                                      const double u = fe_values.shape_value(j,point);
+                                      copy_data.cell_matrix.back()(i,j)
+                                      += (u * v * weight * coefficient_values[point]);
+                                    }
+
+                                copy_data.cell_vector.back()(i) += v *
+                                                  rhs_values_system[point](
+                                                    component_mapping[fe.system_to_component_index(i).first]) * weight;
+                              }
                           }
                       }
-                }
-              else
-                // FE is a scalar one
-                {
-                  rhs_values_scalar.resize (fe_values.n_quadrature_points);
-                  boundary_functions.find(cell->face(face)->boundary_indicator())
-                  ->second->value_list (fe_values.get_quadrature_points(), rhs_values_scalar);
-
-                  if (coefficient != 0)
+                    else
+                      {
+                        coefficient_vector_values.resize (fe_values.n_quadrature_points,
+                                                          Vector<double>(n_components));
+                        coefficient->vector_value_list (fe_values.get_quadrature_points(),
+                                                        coefficient_vector_values);
+                        for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+                          {
+                            const double weight = fe_values.JxW(point);
+                            for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
+                              {
+                                const double v = fe_values.shape_value(i,point);
+                                const unsigned int component_i=
+                                  fe.system_to_component_index(i).first;
+                                for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
+                                  if (fe.system_to_component_index(j).first ==
+                                      component_i)
+                                    {
+                                      const double u = fe_values.shape_value(j,point);
+                                      copy_data.cell_matrix.back()(i,j) +=
+                                        (u * v * weight * coefficient_vector_values[point](component_i));
+                                    }
+                                copy_data.cell_vector.back()(i) += v * 
+                                  rhs_values_system[point](component_mapping[component_i]) * weight;
+                              }
+                          }
+                      }
+                  }
+                else  //      if (coefficient == 0)
+                  for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
                     {
-                      coefficient_values.resize (fe_values.n_quadrature_points);
-                      coefficient->value_list (fe_values.get_quadrature_points(),
-                                               coefficient_values);
-                      for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+                      const double weight = fe_values.JxW(point);
+                      for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
                         {
-                          const double weight = fe_values.JxW(point);
-                          for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
-                            {
-                              const double v = fe_values.shape_value(i,point);
-                              for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
-                                {
-                                  const double u = fe_values.shape_value(j,point);
-                                  cell_matrix(i,j) += (u * v * weight * coefficient_values[point]);
-                                }
-                              cell_vector(i) += v * rhs_values_scalar[point] *weight;
-                            }
+                          const double v = fe_values.shape_value(i,point);
+                          for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
+                            if (fe.system_to_component_index(i).first ==
+                                fe.system_to_component_index(j).first)
+                              {
+                                const double u = fe_values.shape_value(j,point);
+                                copy_data.cell_matrix.back()(i,j) += (u * v * weight);
+                              }
+                          copy_data.cell_vector.back()(i) += v *
+                                            rhs_values_system[point](
+                                              fe.system_to_component_index(i).first) *
+                                            weight;
                         }
                     }
-                  else
+              }
+            else
+              // FE is a scalar one
+              {
+                rhs_values_scalar.resize (fe_values.n_quadrature_points);
+                boundary_functions.find(cell->face(face)->boundary_indicator())
+                ->second->value_list (fe_values.get_quadrature_points(), rhs_values_scalar);
+
+                if (coefficient != 0)
+                  {
+                    coefficient_values.resize (fe_values.n_quadrature_points);
+                    coefficient->value_list (fe_values.get_quadrature_points(),
+                                             coefficient_values);
                     for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
                       {
                         const double weight = fe_values.JxW(point);
@@ -1438,164 +1426,154 @@ namespace MatrixCreator
                             for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
                               {
                                 const double u = fe_values.shape_value(j,point);
-                                cell_matrix(i,j) += (u * v * weight);
+                                copy_data.cell_matrix.back()(i,j) += (u * v * weight * 
+                                    coefficient_values[point]);
                               }
-                            cell_vector(i) += v * rhs_values_scalar[point] * weight;
+                            copy_data.cell_vector.back()(i) += v * rhs_values_scalar[point] *weight;
                           }
                       }
-                }
+                  }
+                else
+                  for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+                    {
+                      const double weight = fe_values.JxW(point);
+                      for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
+                        {
+                          const double v = fe_values.shape_value(i,point);
+                          for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
+                            {
+                              const double u = fe_values.shape_value(j,point);
+                              copy_data.cell_matrix.back()(i,j) += (u * v * weight);
+                            }
+                          copy_data.cell_vector.back()(i) += v * rhs_values_scalar[point] * weight;
+                        }
+                    }
+              }
+
+            cell->face(face)->get_dof_indices (dofs_on_face_vector,
+                                               cell->active_fe_index());
+            // for each dof on the cell, have a
+            // flag whether it is on the face
+            copy_data.dof_is_on_face.push_back(std::vector<bool> (copy_data.dofs_per_cell));
+            // check for each of the dofs on this cell
+            // whether it is on the face
+            for (unsigned int i=0; i<copy_data.dofs_per_cell; ++i)
+              copy_data.dof_is_on_face.back()[i] = (std::find(dofs_on_face_vector.begin(),
+                                                              dofs_on_face_vector.end(),
+                                                              copy_data.dofs[i])
+                                                    !=
+                                                    dofs_on_face_vector.end());
+          }
+    }
+
 
-              // now transfer cell matrix and vector
-              // to the whole boundary matrix
-              //
-              // in the following: dof[i] holds the
-              // global index of the i-th degree of
-              // freedom on the present cell. If it
-              // is also a dof on the boundary, it
-              // must be a nonzero entry in the
-              // dof_to_boundary_mapping and then
-              // the boundary index of this dof is
-              // dof_to_boundary_mapping[dof[i]].
-              //
-              // if dof[i] is not on the boundary,
-              // it should be zero on the boundary
-              // therefore on all quadrature
-              // points and finally all of its
-              // entries in the cell matrix and
-              // vector should be zero. If not, we
-              // throw an error (note: because of
-              // the evaluation of the shape
-              // functions only up to machine
-              // precision, the term "must be zero"
-              // really should mean: "should be
-              // very small". since this is only an
-              // assertion and not part of the
-              // code, we may choose "very small"
-              // quite arbitrarily)
-              //
-              // the main problem here is that the
-              // matrix or vector entry should also
-              // be zero if the degree of freedom
-              // dof[i] is on the boundary, but not
-              // on the present face, i.e. on
-              // another face of the same cell also
-              // on the boundary. We can therefore
-              // not rely on the
-              // dof_to_boundary_mapping[dof[i]]
-              // being !=-1, we really have to
-              // determine whether dof[i] is a
-              // dof on the present face. We do so
-              // by getting the dofs on the
-              // face into @p{dofs_on_face_vector},
-              // a vector as always. Usually,
-              // searching in a vector is
-              // inefficient, so we copy the dofs
-              // into a set, which enables binary
-              // searches.
-              dofs.resize (dofs_per_cell);
-              dofs_on_face_vector.resize (dofs_per_face);
-              dof_is_on_face.resize (dofs_per_cell);
-
-              cell->get_dof_indices (dofs);
-              cell->face(face)->get_dof_indices (dofs_on_face_vector,
-                                                 cell->active_fe_index());
-
-              // check for each of the
-              // dofs on this cell
-              // whether it is on the
-              // face
-              for (unsigned int i=0; i<dofs_per_cell; ++i)
-                dof_is_on_face[i] = (std::find(dofs_on_face_vector.begin(),
-                                               dofs_on_face_vector.end(),
-                                               dofs[i])
-                                     !=
-                                     dofs_on_face_vector.end());
-
-              // in debug mode: compute an element
-              // in the matrix which is
-              // guaranteed to belong to a boundary
-              // dof. We do this to check that the
-              // entries in the cell matrix are
-              // guaranteed to be zero if the
-              // respective dof is not on the
-              // boundary. Since because of
-              // round-off, the actual
-              // value of the matrix entry may be
-              // only close to zero, we assert that
-              // it is small relative to an element
-              // which is guaranteed to be nonzero.
-              // (absolute smallness does not
-              // suffice since the size of the
-              // domain scales in here)
-              //
-              // for this purpose we seek the
-              // diagonal of the matrix, where there
-              // must be an element belonging to
-              // the boundary. we take the maximum
-              // diagonal entry.
+
+    template <int dim,int spacedim>
+    void copy_hp_boundary_mass_matrix_1(MatrixCreator::internal::AssemblerBoundary
+                                        ::CopyData<hp::DoFHandler<dim,spacedim> > const &copy_data,
+                                        typename FunctionMap<spacedim>::type const &boundary_functions,
+                                        std::vector<types::global_dof_index> const &dof_to_boundary_mapping,
+                                        SparseMatrix<double> &matrix,
+                                        Vector<double> &rhs_vector)
+    {
+      // now transfer cell matrix and vector to the whole boundary matrix
+      //
+      // in the following: dof[i] holds the  global index of the i-th degree of
+      // freedom on the present cell. If it is also a dof on the boundary, it
+      // must be a nonzero entry in the dof_to_boundary_mapping and then
+      // the boundary index of this dof is dof_to_boundary_mapping[dof[i]].
+      //
+      // if dof[i] is not on the boundary, it should be zero on the boundary
+      // therefore on all quadrature points and finally all of its
+      // entries in the cell matrix and vector should be zero. If not, we
+      // throw an error (note: because of the evaluation of the shape
+      // functions only up to machine precision, the term "must be zero"
+      // really should mean: "should be very small". since this is only an
+      // assertion and not part of the code, we may choose "very small"
+      // quite arbitrarily)
+      //
+      // the main problem here is that the matrix or vector entry should also
+      // be zero if the degree of freedom dof[i] is on the boundary, but not
+      // on the present face, i.e. on another face of the same cell also
+      // on the boundary. We can therefore not rely on the
+      // dof_to_boundary_mapping[dof[i]] being !=-1, we really have to
+      // determine whether dof[i] is a dof on the present face. We do so
+      // by getting the dofs on the face into @p{dofs_on_face_vector},
+      // a vector as always. Usually, searching in a vector is
+      // inefficient, so we copy the dofs into a set, which enables binary
+      // searches.
+      unsigned int pos(0);
+      for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+      {
+        // check if this face is on that part of
+        // the boundary we are interested in
+        if (boundary_functions.find(copy_data.cell->face(face)->boundary_indicator()) !=
+            boundary_functions.end())
+          {
 #ifdef DEBUG
-              double max_diag_entry = 0;
-              for (unsigned int i=0; i<dofs_per_cell; ++i)
-                if (std::fabs(cell_matrix(i,i)) > max_diag_entry)
-                  max_diag_entry = std::fabs(cell_matrix(i,i));
+            // in debug mode: compute an element in the matrix which is
+            // guaranteed to belong to a boundary dof. We do this to check that the
+            // entries in the cell matrix are guaranteed to be zero if the
+            // respective dof is not on the boundary. Since because of
+            // round-off, the actual value of the matrix entry may be
+            // only close to zero, we assert that it is small relative to an element
+            // which is guaranteed to be nonzero. (absolute smallness does not
+            // suffice since the size of the domain scales in here)
+            //
+            // for this purpose we seek the diagonal of the matrix, where there
+            // must be an element belonging to the boundary. we take the maximum
+            // diagonal entry.
+            types::global_dof_index max_element = static_cast<types::global_dof_index>(0);
+            for (std::vector<types::global_dof_index>::const_iterator i=dof_to_boundary_mapping.begin();
+                 i!=dof_to_boundary_mapping.end(); ++i)
+              if ((*i != hp::DoFHandler<dim,spacedim>::invalid_dof_index) &&
+                  (*i > max_element))
+                max_element = *i;
+            Assert (max_element  == matrix.n()-1, ExcInternalError());
+
+            double max_diag_entry = 0;
+            for (unsigned int i=0; i<copy_data.dofs_per_cell; ++i)
+              if (std::fabs(copy_data.cell_matrix[pos](i,i)) > max_diag_entry)
+                max_diag_entry = std::fabs(copy_data.cell_matrix[pos](i,i));
 #endif
 
-              // lock the matrix
-              Threads::Mutex::ScopedLock lock (mutex);
-              for (unsigned int i=0; i<dofs_per_cell; ++i)
-                for (unsigned int j=0; j<dofs_per_cell; ++j)
-                  if (dof_is_on_face[i] && dof_is_on_face[j])
-                    matrix.add(dof_to_boundary_mapping[dofs[i]],
-                               dof_to_boundary_mapping[dofs[j]],
-                               cell_matrix(i,j));
-                  else
-                    {
-                      // assume that all
-                      // shape functions
-                      // that are nonzero
-                      // on the boundary
-                      // are also listed
-                      // in the
-                      // @p{dof_to_boundary}
-                      // mapping. if that
-                      // is not the case,
-                      // then the
-                      // boundary mass
-                      // matrix does not
-                      // make that much
-                      // sense anyway, as
-                      // it only contains
-                      // entries for
-                      // parts of the
-                      // functions living
-                      // on the boundary
-                      //
-                      // these, we may
-                      // compare here for
-                      // relative
-                      // smallness of all
-                      // entries in the
-                      // local matrix
-                      // which are not
-                      // taken over to
-                      // the global one
-                      Assert (std::fabs(cell_matrix(i,j)) <= 1e-10 * max_diag_entry,
-                              ExcInternalError ());
-                    };
-
-              for (unsigned int j=0; j<dofs_per_cell; ++j)
-                if (dof_is_on_face[j])
-                  rhs_vector(dof_to_boundary_mapping[dofs[j]]) += cell_vector(j);
+            for (unsigned int i=0; i<copy_data.dofs_per_cell; ++i)
+              for (unsigned int j=0; j<copy_data.dofs_per_cell; ++j)
+              {
+                if (copy_data.dof_is_on_face[pos][i] && copy_data.dof_is_on_face[pos][j])
+                  matrix.add(dof_to_boundary_mapping[copy_data.dofs[i]],
+                             dof_to_boundary_mapping[copy_data.dofs[j]],
+                             copy_data.cell_matrix[pos](i,j));
                 else
                   {
-                    // compare here for relative
-                    // smallness
-                    Assert (std::fabs(cell_vector(j)) <= 1e-10 * max_diag_entry,
-                            ExcInternalError());
+                    // assume that all shape functions that are nonzero on the boundary
+                    // are also listed in the @p{dof_to_boundary} mapping. if that
+                    // is not the case, then the boundary mass matrix does not
+                    // make that much sense anyway, as it only contains entries for
+                    // parts of the functions living on the boundary
+                    //
+                    // these, we may compare here for relative smallness of all
+                    // entries in the local matrix which are not taken over to
+                    // the global one
+                    Assert (std::fabs(copy_data.cell_matrix[pos](i,j)) <= 1e-10 * max_diag_entry,
+                            ExcInternalError ());
                   }
-            }
-    }
+              }
 
+            for (unsigned int j=0; j<copy_data.dofs_per_cell; ++j)
+              if (copy_data.dof_is_on_face[pos][j])
+                rhs_vector(dof_to_boundary_mapping[copy_data.dofs[j]]) += copy_data.cell_vector[pos](j);
+              else
+                {
+                  // compare here for relative
+                  // smallness
+                  Assert (std::fabs(copy_data.cell_vector[pos](j)) <= 1e-10 * max_diag_entry,
+                          ExcInternalError());
+                }
+            ++pos;
+          }
+      }
+    }
   }
 
 
@@ -1660,46 +1638,24 @@ namespace MatrixCreator
     else
       AssertDimension (n_components, component_mapping.size());
 
-    const unsigned int n_threads = multithread_info.n_threads();
-    Threads::ThreadGroup<> threads;
-
-    // define starting and end point
-    // for each thread
-    typedef typename hp::DoFHandler<dim,spacedim>::active_cell_iterator active_cell_iterator;
-    std::vector<std::pair<active_cell_iterator,active_cell_iterator> > thread_ranges
-      = Threads::split_range<active_cell_iterator> (dof.begin_active(),
-                                                    dof.end(), n_threads);
-
-    typedef std_cxx1x::tuple<const hp::MappingCollection<dim,spacedim> &,
-            const hp::DoFHandler<dim,spacedim> &,
-            const hp::QCollection<dim-1>&> Commons;
-
-    // mutex to synchronise access to
-    // the matrix
-    Threads::Mutex mutex;
-
-    // then assemble in parallel
-    typedef void (*create_boundary_mass_matrix_1_t)
-    (Commons,
-     SparseMatrix<double>      &matrix,
-     const typename FunctionMap<spacedim>::type &boundary_functions,
-     Vector<double>            &rhs_vector,
-     std::vector<types::global_dof_index> &dof_to_boundary_mapping,
-     const Function<spacedim> *const coefficient,
-     const std::vector<unsigned int> &component_mapping,
-     const MatrixCreator::internal::IteratorRange<hp::DoFHandler<dim,spacedim> >   range,
-     Threads::Mutex      &mutex);
-    create_boundary_mass_matrix_1_t p = &create_boundary_mass_matrix_1<dim,spacedim>;
-
-//TODO: Use WorkStream here
-    for (unsigned int thread=0; thread<n_threads; ++thread)
-      threads += Threads::new_thread (p,
-                                      Commons(mapping, dof, q), matrix,
-                                      boundary_functions, rhs_vector,
-                                      dof_to_boundary_mapping, coefficient,
-                                      component_mapping,
-                                      thread_ranges[thread], mutex);
-    threads.join_all ();
+    MatrixCreator::internal::AssemblerBoundary::Scratch scratch;
+    MatrixCreator::internal::AssemblerBoundary::CopyData<hp::DoFHandler<dim,spacedim> > copy_data;
+
+    WorkStream::run(dof.begin_active(),dof.end(),
+        static_cast<std_cxx1x::function<void (typename hp::DoFHandler<dim,spacedim>::active_cell_iterator
+          const &,MatrixCreator::internal::AssemblerBoundary::Scratch const &,
+          MatrixCreator::internal::AssemblerBoundary::CopyData<hp::DoFHandler<dim,spacedim> > &)> > 
+        (std_cxx1x::bind(create_hp_boundary_mass_matrix_1<dim,spacedim>,std_cxx1x::_1,std_cxx1x::_2,
+                         std_cxx1x::_3,
+                         std_cxx1x::cref(mapping),std_cxx1x::cref(fe_collection),std_cxx1x::cref(q),
+                         std_cxx1x::cref(boundary_functions),coefficient,
+                         std_cxx1x::cref(component_mapping))), 
+        static_cast<std_cxx1x::function<void (MatrixCreator::internal::AssemblerBoundary
+          ::CopyData<hp::DoFHandler<dim,spacedim> > const &)> > (std_cxx1x::bind(
+              copy_hp_boundary_mass_matrix_1<dim,spacedim>,std_cxx1x::_1,
+            std_cxx1x::cref(boundary_functions),std_cxx1x::cref(dof_to_boundary_mapping),
+            std_cxx1x::ref(matrix),std_cxx1x::ref(rhs_vector))),
+        scratch,copy_data);
   }
 
 

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.