#include <base/exceptions.h>
+#include <base/thread_management.h>
#include <map>
* 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.
*/
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}
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,
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}
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,
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,
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,
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}
* 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);
};
#include <base/function.h>
+#include <base/multithread_info.h>
#include <dofs/dof_handler.h>
#include <dofs/dof_accessor.h>
#include <grid/tria_iterator.h>
-// 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)
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);
// 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 ();
};
};
}
-// 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 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 |
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);
// 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 ();
};
};
#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();
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
//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()) &&
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,
Assert (fabs(cell_vector(j)) <= 1e-10 * max_diag_entry,
ExcInternalError());
};
+ mutex.release ();
};
};
}
-// 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);
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);
// 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 ();
};
};
-// 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 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 |
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);
// 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 ();
};
};
}
+
template <int dim>
template <typename number>
void