]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Make the functions in MatrixCreator use multiple threads if supported.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 9 Apr 2001 09:38:34 +0000 (09:38 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 9 Apr 2001 09:38:34 +0000 (09:38 +0000)
git-svn-id: https://svn.dealii.org/trunk@4401 0785d39b-7218-0410-832d-ea1e28bc413d

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

index d0fdc993f72520003573fad03e4a4eb97f37b3b7..9000db266a373b5d6bf1e37e32732b144ce952b7 100644 (file)
@@ -11,6 +11,7 @@
 
 
 #include <base/exceptions.h>
+#include <base/thread_management.h>
 #include <map>
 
 
@@ -205,7 +206,11 @@ class MatrixCreator
                                      * Assemble the mass matrix. If no 
                                      * coefficient is given, it is assumed
                                      * to be unity.
-                                     * 
+                                     *
+                                     * If the library is configured
+                                     * to use multithreading, this
+                                     * function works in parallel.
+                                     *
                                      * See the general doc of this class
                                      * for more information.
                                      */
@@ -213,7 +218,7 @@ class MatrixCreator
                                    const DoFHandler<dim>    &dof,
                                    const Quadrature<dim>    &q,
                                    SparseMatrix<double>     &matrix,
-                                   const Function<dim>      *a = 0);
+                                   const Function<dim> * const a = 0);
 
                                     /**
                                      * Calls the @p{create_mass_matrix}
@@ -223,16 +228,20 @@ class MatrixCreator
     static void create_mass_matrix (const DoFHandler<dim>    &dof,
                                    const Quadrature<dim>    &q,
                                    SparseMatrix<double>     &matrix,
-                                   const Function<dim>      *a = 0);
+                                   const Function<dim> * const a = 0);
 
                                     /**
-                                     * Assemble the mass matrix and a right
-                                     * hand side vector. If no 
-                                     * coefficient is given, it is assumed
-                                     * to be unity.
+                                     * Assemble the mass matrix and a
+                                     * right hand side vector. If no
+                                     * coefficient is given, it is
+                                     * assumed to be unity.
                                      *
-                                     * See the general doc of this class
-                                     * for more information.
+                                     * If the library is configured
+                                     * to use multithreading, this
+                                     * function works in parallel.
+                                     *
+                                     * See the general doc of this
+                                     * class for more information.
                                      */
     static void create_mass_matrix (const Mapping<dim>       &mapping,
                                    const DoFHandler<dim>    &dof,
@@ -240,7 +249,7 @@ class MatrixCreator
                                    SparseMatrix<double>     &matrix,
                                    const Function<dim>      &rhs,
                                    Vector<double>           &rhs_vector,
-                                   const Function<dim>      *a = 0);
+                                   const Function<dim> * const a = 0);
 
                                     /**
                                      * Calls the @p{create_mass_matrix}
@@ -252,52 +261,65 @@ class MatrixCreator
                                    SparseMatrix<double>     &matrix,
                                    const Function<dim>      &rhs,
                                    Vector<double>           &rhs_vector,
-                                   const Function<dim>      *a = 0);
+                                   const Function<dim> * const a = 0);
     
                                     /**
-                                     * Assemble the mass matrix and a right
-                                     * hand side vector along the boundary.
-                                     * If no 
-                                     * coefficient is given, it is assumed
-                                     * to be constant one.
+                                     * Assemble the mass matrix and a
+                                     * right hand side vector along
+                                     * the boundary.  If no
+                                     * coefficient is given, it is
+                                     * assumed to be constant one.
                                      *
-                                     * The matrix is assumed to already be
-                                     * initialized with a suiting sparsity
-                                     * pattern (the @ref{DoFHandler} provides an
+                                     * The matrix is assumed to
+                                     * already be initialized with a
+                                     * suiting sparsity pattern (the
+                                     * @ref{DoFHandler} provides an
                                      * appropriate function).
                                      *
-                                     * See the general doc of this class
-                                     * for more information.
+                                     * If the library is configured
+                                     * to use multithreading, this
+                                     * function works in parallel.
+                                     *
+                                     * See the general doc of this
+                                     * class for more information.
                                      */
-    static void create_boundary_mass_matrix (const Mapping<dim>       &mapping,
-                                            const DoFHandler<dim>    &dof,
-                                            const Quadrature<dim-1>  &q,
-                                            SparseMatrix<double>     &matrix,
-                                            const FunctionMap        &boundary_functions,
-                                            Vector<double>           &rhs_vector,
-                                            std::vector<unsigned int>&dof_to_boundary_mapping,
-                                            const Function<dim>      *a = 0);
+    static
+    void create_boundary_mass_matrix (const Mapping<dim>       &mapping,
+                                     const DoFHandler<dim>    &dof,
+                                     const Quadrature<dim-1>  &q,
+                                     SparseMatrix<double>     &matrix,
+                                     const FunctionMap        &boundary_functions,
+                                     Vector<double>           &rhs_vector,
+                                     std::vector<unsigned int>&dof_to_boundary_mapping,
+                                     const Function<dim>      *a = 0);
 
                                     /**
-                                     * Calls the @p{create_boundary_mass_matrix}
+                                     * Calls the
+                                     * @p{create_boundary_mass_matrix}
                                      * function, see above, with
                                      * @p{mapping=MappingQ1<dim>()}.
                                      */
-    static void create_boundary_mass_matrix (const DoFHandler<dim>    &dof,
-                                            const Quadrature<dim-1>  &q,
-                                            SparseMatrix<double>     &matrix,
-                                            const FunctionMap        &boundary_functions,
-                                            Vector<double>           &rhs_vector,
-                                            std::vector<unsigned int>&dof_to_boundary_mapping,
-                                            const Function<dim>      *a = 0);
+    static
+    void create_boundary_mass_matrix (const DoFHandler<dim>    &dof,
+                                     const Quadrature<dim-1>  &q,
+                                     SparseMatrix<double>     &matrix,
+                                     const FunctionMap        &boundary_functions,
+                                     Vector<double>           &rhs_vector,
+                                     std::vector<unsigned int>&dof_to_boundary_mapping,
+                                     const Function<dim>      *a = 0);
 
                                     /**
-                                     * Assemble the Laplace matrix. If no 
-                                     * coefficient is given, it is assumed
-                                     * to be constant one.
+                                     * Assemble the Laplace
+                                     * matrix. If no coefficient is
+                                     * given, it is assumed to be
+                                     * constant one.
                                      * 
-                                     * See the general doc of this class
-                                     * for more information.
+                                     * If the library is configured
+                                     * to use multithreading, this
+                                     * function works in parallel.
+                                     *
+                                     * See the general doc of this
+                                     * class for more information.
                                      */
     static void create_laplace_matrix (const Mapping<dim>       &mapping,
                                       const DoFHandler<dim>    &dof,
@@ -316,10 +338,11 @@ class MatrixCreator
                                       const Function<dim>      *a = 0);
 
                                     /**
-                                     * Generate Laplace matrix for a given level.
+                                     * Generate Laplace matrix for a
+                                     * given level.
                                      * 
-                                     * See the general doc of this class
-                                     * for more information.
+                                     * See the general doc of this
+                                     * class for more information.
                                      */
     static void create_level_laplace_matrix (unsigned int             level,
                                             const MGDoFHandler<dim>& dof,
@@ -328,13 +351,18 @@ class MatrixCreator
                                             const Function<dim>*     a = 0);
 
                                     /**
-                                     * Assemble the Laplace matrix and a right
-                                     * hand side vector. If no 
-                                     * coefficient is given, it is assumed
-                                     * to be constant one.
+                                     * Assemble the Laplace matrix
+                                     * and a right hand side
+                                     * vector. If no coefficient is
+                                     * given, it is assumed to be
+                                     * constant one.
                                      * 
-                                     * See the general doc of this class
-                                     * for more information.
+                                     * If the library is configured
+                                     * to use multithreading, this
+                                     * function works in parallel.
+                                     *
+                                     * See the general doc of this
+                                     * class for more information.
                                      */
     static void create_laplace_matrix (const Mapping<dim>       &mapping,
                                       const DoFHandler<dim>    &dof,
@@ -342,7 +370,7 @@ class MatrixCreator
                                       SparseMatrix<double>     &matrix,
                                       const Function<dim>      &rhs,
                                       Vector<double>           &rhs_vector,
-                                      const Function<dim>      *a = 0);
+                                      const Function<dim> * const a = 0);
 
                                     /**
                                      * Calls the @p{create_laplace_matrix}
@@ -360,6 +388,122 @@ class MatrixCreator
                                      * Exception
                                      */
     DeclException0 (ExcComponentMismatch);
+
+  private:
+                                    /**
+                                     * Convenience abbreviation for
+                                     * DoF handler cell iterators.
+                                     */
+    typedef typename DoFHandler<dim>::active_cell_iterator active_cell_iterator;
+
+                                    /**
+                                     * Pair of iterators denoting a
+                                     * half-open range.
+                                     */
+    typedef std::pair<active_cell_iterator,active_cell_iterator> IteratorRange;
+    
+
+                                    /**
+                                     * Version of the same function
+                                     * (without suffix @p{_1}) with
+                                     * the same argument list that
+                                     * operates only on an interval
+                                     * of iterators. Used for
+                                     * parallelization. The mutex is
+                                     * used to synchronise access to
+                                     * the matrix.
+                                     */
+    static
+    void create_mass_matrix_1 (const Mapping<dim>       &mapping,
+                              const DoFHandler<dim>    &dof,
+                              const Quadrature<dim>    &q,
+                              SparseMatrix<double>     &matrix,
+                              const Function<dim> * const a,
+                              const IteratorRange      &range,
+                              Threads::ThreadMutex     &mutex);
+
+                                    /**
+                                     * Version of the same function
+                                     * (without suffix @p{_2}) with
+                                     * the same argument list that
+                                     * operates only on an interval
+                                     * of iterators. Used for
+                                     * parallelization. The mutex is
+                                     * used to synchronise access to
+                                     * the matrix.
+                                     */
+    static
+    void create_mass_matrix_2 (const Mapping<dim>       &mapping,
+                              const DoFHandler<dim>    &dof,
+                              const Quadrature<dim>    &q,
+                              SparseMatrix<double>     &matrix,
+                              const Function<dim>      &rhs,
+                              Vector<double>           &rhs_vector,
+                              const Function<dim> * const a,
+                              const IteratorRange      &range,
+                              Threads::ThreadMutex     &mutex);
+
+                                    /**
+                                     * Version of the same function
+                                     * (without suffix @p{_1}) with
+                                     * the same argument list that
+                                     * operates only on an interval
+                                     * of iterators. Used for
+                                     * parallelization. The mutex is
+                                     * used to synchronise access to
+                                     * the matrix.
+                                     */
+    static
+    void create_laplace_matrix_1 (const Mapping<dim>       &mapping,
+                                 const DoFHandler<dim>    &dof,
+                                 const Quadrature<dim>    &q,
+                                 SparseMatrix<double>     &matrix,
+                                 const Function<dim> * const a,
+                                 const IteratorRange      &range,
+                                 Threads::ThreadMutex     &mutex);
+
+                                    /**
+                                     * Version of the same function
+                                     * (without suffix @p{_2}) with
+                                     * the same argument list that
+                                     * operates only on an interval
+                                     * of iterators. Used for
+                                     * parallelization. The mutex is
+                                     * used to synchronise access to
+                                     * the matrix.
+                                     */
+    static
+    void create_laplace_matrix_2 (const Mapping<dim>       &mapping,
+                                 const DoFHandler<dim>    &dof,
+                                 const Quadrature<dim>    &q,
+                                 SparseMatrix<double>     &matrix,
+                                 const Function<dim>      &rhs,
+                                 Vector<double>           &rhs_vector,
+                                 const Function<dim> * const a,
+                                 const IteratorRange      &range,
+                                 Threads::ThreadMutex     &mutex);
+
+                                    /**
+                                     * Version of the same function
+                                     * (without suffix @p{_1}) with
+                                     * the same argument list that
+                                     * operates only on an interval
+                                     * of iterators. Used for
+                                     * parallelization. The mutex is
+                                     * used to synchronise access to
+                                     * the matrix.
+                                     */
+    static
+    void create_boundary_mass_matrix_1 (const Mapping<dim>       &mapping,
+                                       const DoFHandler<dim>    &dof,
+                                       const Quadrature<dim-1>  &q,
+                                       SparseMatrix<double>     &matrix,
+                                       const FunctionMap        &boundary_functions,
+                                       Vector<double>           &rhs_vector,
+                                       std::vector<unsigned int>&dof_to_boundary_mapping,
+                                       const Function<dim> * const a,
+                                       const IteratorRange      &range,
+                                       Threads::ThreadMutex     &mutex);
 };
 
 
index 64bd7d2a6b7ca6bda604773aca0647e05725fb11..88334a0133857b42982d76ef2092eb62a25c090e 100644 (file)
@@ -13,6 +13,7 @@
 
 
 #include <base/function.h>
+#include <base/multithread_info.h>
 #include <dofs/dof_handler.h>
 #include <dofs/dof_accessor.h>
 #include <grid/tria_iterator.h>
@@ -40,14 +41,46 @@ using namespace std;
 
 
 
-// TODO:[RH, WB] extend this function to use vector valued coefficient functions for system elements.
-// TODO:[WB] implement multithreading for this function.
 template <int dim>
 void MatrixCreator<dim>::create_mass_matrix (const Mapping<dim>       &mapping,
                                             const DoFHandler<dim>    &dof,
                                             const Quadrature<dim>    &q,
                                             SparseMatrix<double>     &matrix,
                                             const Function<dim> * const coefficient)
+{
+  const unsigned int n_threads = multithread_info.n_default_threads;
+  Threads::ThreadManager thread_manager;
+
+                                  // define starting and end point
+                                  // for each thread
+  std::vector<IteratorRange> thread_ranges
+    = Threads::split_range<active_cell_iterator> (dof.begin_active(),
+                                                 dof.end(), n_threads);
+
+                                  // mutex to synchronise access to
+                                  // the matrix
+  Threads::ThreadMutex mutex;
+  
+                                  // then assemble in parallel
+  for (unsigned int thread=0; thread<n_threads; ++thread)
+    Threads::spawn (thread_manager,
+                   Threads::encapsulate(&MatrixCreator<dim>::create_mass_matrix_1)
+                   .collect_args (mapping, dof, q, matrix, coefficient,
+                                  thread_ranges[thread], mutex));
+  thread_manager.wait ();  
+};
+
+
+
+// TODO:[RH] extend this function to use vector valued coefficient functions for system elements.
+template <int dim>
+void MatrixCreator<dim>::create_mass_matrix_1 (const Mapping<dim>       &mapping,
+                                              const DoFHandler<dim>    &dof,
+                                              const Quadrature<dim>    &q,
+                                              SparseMatrix<double>     &matrix,
+                                              const Function<dim> * const coefficient,
+                                              const IteratorRange      &range,
+                                              Threads::ThreadMutex     &mutex)
 {
   UpdateFlags update_flags = UpdateFlags(update_values | update_JxW_values);
   if (coefficient != 0)
@@ -65,8 +98,8 @@ void MatrixCreator<dim>::create_mass_matrix (const Mapping<dim>       &mapping,
   
   std::vector<unsigned int> dof_indices (dofs_per_cell);
   
-  typename DoFHandler<dim>::active_cell_iterator cell = dof.begin_active();
-  for (; cell!=dof.end(); ++cell)
+  typename DoFHandler<dim>::active_cell_iterator cell = range.first;
+  for (; cell!=range.second; ++cell)
     {
       fe_values.reinit (cell);
       
@@ -104,10 +137,15 @@ void MatrixCreator<dim>::create_mass_matrix (const Mapping<dim>       &mapping,
 
                                       // transfer everything into the
                                       // global object
+      mutex.acquire ();
       for (unsigned int i=0; i<dofs_per_cell; ++i)
        for (unsigned int j=0; j<dofs_per_cell; ++j)
-         matrix.add (dof_indices[i], dof_indices[j],
-                     cell_matrix(i,j));
+         if ((n_components==1) ||
+             (fe.system_to_component_index(i).first ==
+              fe.system_to_component_index(j).first))
+           matrix.add (dof_indices[i], dof_indices[j],
+                       cell_matrix(i,j));
+      mutex.release ();
     };
 };
 
@@ -123,8 +161,6 @@ void MatrixCreator<dim>::create_mass_matrix (const DoFHandler<dim>    &dof,
 }
 
 
-// TODO:[RH, WB] extend this function to use vector valued coefficient functions for system elements.
-// TODO:[WB] implement multithreading for this function.
 template <int dim>
 void MatrixCreator<dim>::create_mass_matrix (const Mapping<dim>       &mapping,
                                             const DoFHandler<dim>    &dof,
@@ -133,6 +169,45 @@ void MatrixCreator<dim>::create_mass_matrix (const Mapping<dim>       &mapping,
                                             const Function<dim>      &rhs,
                                             Vector<double>           &rhs_vector,
                                             const Function<dim> * const coefficient)
+{
+  const unsigned int n_threads = multithread_info.n_default_threads;
+  Threads::ThreadManager thread_manager;
+
+                                  // define starting and end point
+                                  // for each thread
+  std::vector<IteratorRange> thread_ranges
+    = Threads::split_range<active_cell_iterator> (dof.begin_active(),
+                                                 dof.end(), n_threads);
+
+                                  // mutex to synchronise access to
+                                  // the matrix
+  Threads::ThreadMutex mutex;
+  
+                                  // then assemble in parallel
+  for (unsigned int thread=0; thread<n_threads; ++thread)
+    Threads::spawn (thread_manager,
+                   Threads::encapsulate(&MatrixCreator<dim>::
+                                        create_mass_matrix_2)
+                   .collect_args (mapping, dof, q, matrix, rhs,
+                                  rhs_vector, coefficient,
+                                  thread_ranges[thread], mutex));
+  thread_manager.wait ();  
+};
+
+
+
+// TODO:[RH] extend this function to use vector valued coefficient functions for system elements.
+template <int dim>
+void
+MatrixCreator<dim>::create_mass_matrix_2 (const Mapping<dim>       &mapping,
+                                         const DoFHandler<dim>    &dof,
+                                         const Quadrature<dim>    &q,
+                                         SparseMatrix<double>     &matrix,
+                                         const Function<dim>      &rhs,
+                                         Vector<double>           &rhs_vector,
+                                         const Function<dim> * const coefficient,
+                                         const IteratorRange      &range,
+                                         Threads::ThreadMutex     &mutex)
 {
   UpdateFlags update_flags = UpdateFlags(update_values |
                                         update_q_points |
@@ -154,8 +229,8 @@ void MatrixCreator<dim>::create_mass_matrix (const Mapping<dim>       &mapping,
   
   std::vector<unsigned int> dof_indices (dofs_per_cell);
   
-  typename DoFHandler<dim>::active_cell_iterator cell = dof.begin_active();
-  for (; cell!=dof.end(); ++cell)
+  typename DoFHandler<dim>::active_cell_iterator cell = range.first;
+  for (; cell!=range.second; ++cell)
     {
       fe_values.reinit (cell);
       
@@ -205,12 +280,17 @@ void MatrixCreator<dim>::create_mass_matrix (const Mapping<dim>       &mapping,
 
                                       // transfer everything into the
                                       // global object
+      mutex.acquire ();
       for (unsigned int i=0; i<dofs_per_cell; ++i)
        for (unsigned int j=0; j<dofs_per_cell; ++j)
-         matrix.add (dof_indices[i], dof_indices[j],
-                     cell_matrix(i,j));
+         if ((n_components==1) ||
+             (fe.system_to_component_index(i).first ==
+              fe.system_to_component_index(j).first))
+           matrix.add (dof_indices[i], dof_indices[j],
+                       cell_matrix(i,j));
       for (unsigned int i=0; i<dofs_per_cell; ++i)
        rhs_vector(dof_indices[i]) += local_rhs(i);
+      mutex.release ();
     };
 };
 
@@ -260,17 +340,79 @@ void MatrixCreator<1>::create_boundary_mass_matrix (const Mapping<1>          &,
 #endif
 
 
-// TODO:[RH, WB] extend this function to use vector valued coefficient functions for system elements.
-// TODO:[WB] implement multithreading for this function.
+// TODO:[RH] extend this function to use vector valued coefficient functions for system elements.
 template <int dim>
-void MatrixCreator<dim>::create_boundary_mass_matrix (const Mapping<dim>        &mapping,
-                                                     const DoFHandler<dim>     &dof,
-                                                     const Quadrature<dim-1>   &q,
-                                                     SparseMatrix<double>      &matrix,
-                                                     const FunctionMap         &boundary_functions,
-                                                     Vector<double>            &rhs_vector,
-                                                     std::vector<unsigned int> &dof_to_boundary_mapping,
-                                                     const Function<dim>       *a)
+void
+MatrixCreator<dim>::create_boundary_mass_matrix (const Mapping<dim>        &mapping,
+                                                const DoFHandler<dim>     &dof,
+                                                const Quadrature<dim-1>   &q,
+                                                SparseMatrix<double>      &matrix,
+                                                const FunctionMap         &boundary_functions,
+                                                Vector<double>            &rhs_vector,
+                                                std::vector<unsigned int> &dof_to_boundary_mapping,
+                                                const Function<dim>       *a)
+{
+  const unsigned int n_threads = multithread_info.n_default_threads;
+  Threads::ThreadManager thread_manager;
+
+                                  // define starting and end point
+                                  // for each thread
+  std::vector<IteratorRange> thread_ranges
+    = Threads::split_range<active_cell_iterator> (dof.begin_active(),
+                                                 dof.end(), n_threads);
+
+                                  // mutex to synchronise access to
+                                  // the matrix
+  Threads::ThreadMutex mutex;
+  
+                                  // then assemble in parallel
+  for (unsigned int thread=0; thread<n_threads; ++thread)
+    Threads::spawn (thread_manager,
+                   Threads::encapsulate(&MatrixCreator<dim>::
+                                        create_boundary_mass_matrix_1)
+                   .collect_args (mapping, dof, q, matrix,
+                                  boundary_functions, rhs_vector,
+                                  dof_to_boundary_mapping, a,
+                                  thread_ranges[thread], mutex));
+  thread_manager.wait ();  
+};
+
+
+#if deal_II_dimension == 1
+
+template <>
+void
+MatrixCreator<1>::
+create_boundary_mass_matrix_1 (const Mapping<1>          &,
+                              const DoFHandler<1>       &,
+                              const Quadrature<0>       &,
+                              SparseMatrix<double>      &,
+                              const FunctionMap         &,
+                              Vector<double>            &,
+                              std::vector<unsigned int> &,
+                              const Function<1> * const ,
+                              const IteratorRange       &,
+                              Threads::ThreadMutex      &)
+{
+  Assert (false, ExcNotImplemented());
+};
+
+#endif
+
+
+template <int dim>
+void
+MatrixCreator<dim>::
+create_boundary_mass_matrix_1 (const Mapping<dim>        &mapping,
+                              const DoFHandler<dim>     &dof,
+                              const Quadrature<dim-1>   &q,
+                              SparseMatrix<double>      &matrix,
+                              const FunctionMap         &boundary_functions,
+                              Vector<double>            &rhs_vector,
+                              std::vector<unsigned int> &dof_to_boundary_mapping,
+                              const Function<dim> * const a,
+                              const IteratorRange       &range,
+                              Threads::ThreadMutex      &mutex)
 {
   const FiniteElement<dim> &fe = dof.get_fe();
   const unsigned int n_components  = fe.n_components();
@@ -325,9 +467,8 @@ void MatrixCreator<dim>::create_boundary_mass_matrix (const Mapping<dim>
   std::vector<unsigned int> dofs_on_face_vector (dofs_per_face);
   std::set<int> dofs_on_face;
 
-  DoFHandler<dim>::active_cell_iterator cell = dof.begin_active (),
-                                       endc = dof.end ();
-  for (; cell!=endc; ++cell)
+  DoFHandler<dim>::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
@@ -516,7 +657,8 @@ void MatrixCreator<dim>::create_boundary_mass_matrix (const Mapping<dim>
 //in the innermost loop, we traverse a set twice in "find", but
 //this could be made much faster, e.g. by first building a vector<bool>
 //that already stores the result of find beforehand
-// (then remove dofs_on_face altogether)         
+// (then remove dofs_on_face altogether)
+         mutex.acquire ();
          for (unsigned int i=0; i<dofs_per_cell; ++i)
            for (unsigned int j=0; j<dofs_per_cell; ++j)
              if ((dofs_on_face.find(dofs[i]) != dofs_on_face.end()) &&
@@ -526,6 +668,7 @@ void MatrixCreator<dim>::create_boundary_mass_matrix (const Mapping<dim>
                           cell_matrix(i,j));
              else
                {
+//TODO:[?] We assume here that shape functions that have support also on the boundary also have their support point on the boundary (by having their indices in dofs_on_face). This is not true, however, e.g. for discontinuous elements.
                                                   // compare here for relative
                                                   // smallness
                  Assert (fabs(cell_matrix(i,j)) <= 1e-10 * max_diag_entry,
@@ -542,6 +685,7 @@ void MatrixCreator<dim>::create_boundary_mass_matrix (const Mapping<dim>
                Assert (fabs(cell_vector(j)) <= 1e-10 * max_diag_entry,
                        ExcInternalError());
              };
+         mutex.release ();
        };
 };
 
@@ -561,14 +705,47 @@ void MatrixCreator<dim>::create_boundary_mass_matrix (const DoFHandler<dim>
 }
 
 
-// TODO:[RH, WB] extend this function to use vector valued coefficient functions for system elements.
-// TODO:[WB] implement multithreading for this function.
+
 template <int dim>
 void MatrixCreator<dim>::create_laplace_matrix (const Mapping<dim>       &mapping,
                                                const DoFHandler<dim>    &dof,
                                                const Quadrature<dim>    &q,
                                                SparseMatrix<double>     &matrix,
                                                const Function<dim> * const coefficient)
+{
+  const unsigned int n_threads = multithread_info.n_default_threads;
+  Threads::ThreadManager thread_manager;
+
+                                  // define starting and end point
+                                  // for each thread
+  std::vector<IteratorRange> thread_ranges
+    = Threads::split_range<active_cell_iterator> (dof.begin_active(),
+                                                 dof.end(), n_threads);
+
+                                  // mutex to synchronise access to
+                                  // the matrix
+  Threads::ThreadMutex mutex;
+  
+                                  // then assemble in parallel
+  for (unsigned int thread=0; thread<n_threads; ++thread)
+    Threads::spawn (thread_manager,
+                   Threads::encapsulate(&MatrixCreator<dim>::create_laplace_matrix_1)
+                   .collect_args (mapping, dof, q, matrix, coefficient,
+                                  thread_ranges[thread], mutex));
+  thread_manager.wait ();  
+};
+
+
+
+// TODO:[RH] extend this function to use vector valued coefficient functions for system elements.
+template <int dim>
+void MatrixCreator<dim>::create_laplace_matrix_1 (const Mapping<dim>       &mapping,
+                                                 const DoFHandler<dim>    &dof,
+                                                 const Quadrature<dim>    &q,
+                                                 SparseMatrix<double>     &matrix,
+                                                 const Function<dim> * const coefficient,
+                                                 const IteratorRange      &range,
+                                                 Threads::ThreadMutex     &mutex)
 {
   UpdateFlags update_flags = UpdateFlags(update_JxW_values |
                                         update_gradients);
@@ -587,8 +764,8 @@ void MatrixCreator<dim>::create_laplace_matrix (const Mapping<dim>       &mappin
   
   std::vector<unsigned int> dof_indices (dofs_per_cell);
   
-  typename DoFHandler<dim>::active_cell_iterator cell = dof.begin_active();
-  for (; cell!=dof.end(); ++cell)
+  typename DoFHandler<dim>::active_cell_iterator cell = range.first;
+  for (; cell!=range.second; ++cell)
     {
       fe_values.reinit (cell);
       
@@ -627,10 +804,15 @@ void MatrixCreator<dim>::create_laplace_matrix (const Mapping<dim>       &mappin
 
                                       // transfer everything into the
                                       // global object
+      mutex.acquire ();
       for (unsigned int i=0; i<dofs_per_cell; ++i)
        for (unsigned int j=0; j<dofs_per_cell; ++j)
-         matrix.add (dof_indices[i], dof_indices[j],
-                     cell_matrix(i,j));
+         if ((n_components==1) ||
+             (fe.system_to_component_index(i).first ==
+              fe.system_to_component_index(j).first))
+           matrix.add (dof_indices[i], dof_indices[j],
+                       cell_matrix(i,j));
+      mutex.release ();
     };
 };
 
@@ -684,8 +866,6 @@ void MatrixCreator<dim>::create_level_laplace_matrix (unsigned int level,
 
 
 
-// TODO:[RH, WB] extend this function to use vector valued coefficient functions for system elements.
-// TODO:[WB] implement multithreading for this function.
 template <int dim>
 void MatrixCreator<dim>::create_laplace_matrix (const Mapping<dim>       &mapping,
                                                const DoFHandler<dim>    &dof,
@@ -694,6 +874,45 @@ void MatrixCreator<dim>::create_laplace_matrix (const Mapping<dim>       &mappin
                                                const Function<dim>      &rhs,
                                                Vector<double>           &rhs_vector,
                                                const Function<dim> * const coefficient)
+{
+  const unsigned int n_threads = multithread_info.n_default_threads;
+  Threads::ThreadManager thread_manager;
+
+                                  // define starting and end point
+                                  // for each thread
+  std::vector<IteratorRange> thread_ranges
+    = Threads::split_range<active_cell_iterator> (dof.begin_active(),
+                                                 dof.end(), n_threads);
+
+                                  // mutex to synchronise access to
+                                  // the matrix
+  Threads::ThreadMutex mutex;
+  
+                                  // then assemble in parallel
+  for (unsigned int thread=0; thread<n_threads; ++thread)
+    Threads::spawn (thread_manager,
+                   Threads::encapsulate(&MatrixCreator<dim>::
+                                        create_laplace_matrix_2)
+                   .collect_args (mapping, dof, q, matrix, rhs,
+                                  rhs_vector, coefficient,
+                                  thread_ranges[thread], mutex));
+  thread_manager.wait ();  
+};
+
+
+
+// TODO:[RH] extend this function to use vector valued coefficient functions for system elements.
+template <int dim>
+void
+MatrixCreator<dim>::create_laplace_matrix_2 (const Mapping<dim>       &mapping,
+                                            const DoFHandler<dim>    &dof,
+                                            const Quadrature<dim>    &q,
+                                            SparseMatrix<double>     &matrix,
+                                            const Function<dim>      &rhs,
+                                            Vector<double>           &rhs_vector,
+                                            const Function<dim> * const coefficient,
+                                            const IteratorRange      &range,
+                                            Threads::ThreadMutex     &mutex)
 {
   UpdateFlags update_flags = UpdateFlags(update_values    |
                                         update_gradients |
@@ -716,8 +935,8 @@ void MatrixCreator<dim>::create_laplace_matrix (const Mapping<dim>       &mappin
   
   std::vector<unsigned int> dof_indices (dofs_per_cell);
   
-  typename DoFHandler<dim>::active_cell_iterator cell = dof.begin_active();
-  for (; cell!=dof.end(); ++cell)
+  typename DoFHandler<dim>::active_cell_iterator cell = range.first;
+  for (; cell!=range.second; ++cell)
     {
       fe_values.reinit (cell);
       
@@ -769,12 +988,17 @@ void MatrixCreator<dim>::create_laplace_matrix (const Mapping<dim>       &mappin
 
                                       // transfer everything into the
                                       // global object
+      mutex.acquire ();
       for (unsigned int i=0; i<dofs_per_cell; ++i)
        for (unsigned int j=0; j<dofs_per_cell; ++j)
-         matrix.add (dof_indices[i], dof_indices[j],
-                     cell_matrix(i,j));
+         if ((n_components==1) ||
+             (fe.system_to_component_index(i).first ==
+              fe.system_to_component_index(j).first))
+           matrix.add (dof_indices[i], dof_indices[j],
+                       cell_matrix(i,j));
       for (unsigned int i=0; i<dofs_per_cell; ++i)
        rhs_vector(dof_indices[i]) += local_rhs(i);
+      mutex.release ();
     };
 };
 
@@ -793,6 +1017,7 @@ void MatrixCreator<dim>::create_laplace_matrix (const DoFHandler<dim>    &dof,
 }
 
 
+
 template <int dim>
 template <typename number>
 void

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.