From: bangerth Date: Thu, 15 Sep 2011 05:28:30 +0000 (+0000) Subject: Make MatrixCreator a namespace. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a0f0aa22d0d0d3ccf1e734b66f0914f05b479c5a;p=dealii-svn.git Make MatrixCreator a namespace. git-svn-id: https://svn.dealii.org/trunk@24325 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/include/deal.II/numerics/matrices.h b/deal.II/include/deal.II/numerics/matrices.h index 9348157e24..bb3ae80fbd 100644 --- a/deal.II/include/deal.II/numerics/matrices.h +++ b/deal.II/include/deal.II/numerics/matrices.h @@ -1,7 +1,7 @@ //--------------------------------------------------------------------------- // $Id$ // -// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2010 by the deal.II authors +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2010, 2011 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -78,10 +78,9 @@ namespace TrilinosWrappers /** - * This class provides functions that assemble certain standard matrices for a + * This namespace provides functions that assemble certain standard matrices for a * given triangulation, using a given finite element, a given mapping and a - * quadrature formula. All functions are static, so it is not necessary to - * create an object of this type. + * quadrature formula. * * *

Conventions for all functions

@@ -208,9 +207,8 @@ namespace TrilinosWrappers * @ingroup numerics * @author Wolfgang Bangerth, 1998, Ralf Hartmann, 2001 */ -class MatrixCreator +namespace MatrixCreator { - public: /** * Assemble the mass matrix. If no * coefficient is given, it is assumed @@ -224,7 +222,7 @@ class MatrixCreator * for more information. */ template - static void create_mass_matrix (const Mapping &mapping, + void create_mass_matrix (const Mapping &mapping, const DoFHandler &dof, const Quadrature &q, SparseMatrix &matrix, @@ -236,7 +234,7 @@ class MatrixCreator * mapping=MappingQ1@(). */ template - static void create_mass_matrix (const DoFHandler &dof, + void create_mass_matrix (const DoFHandler &dof, const Quadrature &q, SparseMatrix &matrix, const Function * const a = 0); @@ -255,7 +253,7 @@ class MatrixCreator * class for more information. */ template - static void create_mass_matrix (const Mapping &mapping, + void create_mass_matrix (const Mapping &mapping, const DoFHandler &dof, const Quadrature &q, SparseMatrix &matrix, @@ -269,7 +267,7 @@ class MatrixCreator * mapping=MappingQ1@(). */ template - static void create_mass_matrix (const DoFHandler &dof, + void create_mass_matrix (const DoFHandler &dof, const Quadrature &q, SparseMatrix &matrix, const Function &rhs, @@ -280,7 +278,7 @@ class MatrixCreator * objects. */ template - static void create_mass_matrix (const hp::MappingCollection &mapping, + void create_mass_matrix (const hp::MappingCollection &mapping, const hp::DoFHandler &dof, const hp::QCollection &q, SparseMatrix &matrix, @@ -291,7 +289,7 @@ class MatrixCreator * objects. */ template - static void create_mass_matrix (const hp::DoFHandler &dof, + void create_mass_matrix (const hp::DoFHandler &dof, const hp::QCollection &q, SparseMatrix &matrix, const Function * const a = 0); @@ -301,7 +299,7 @@ class MatrixCreator * objects. */ template - static void create_mass_matrix (const hp::MappingCollection &mapping, + void create_mass_matrix (const hp::MappingCollection &mapping, const hp::DoFHandler &dof, const hp::QCollection &q, SparseMatrix &matrix, @@ -314,7 +312,7 @@ class MatrixCreator * objects. */ template - static void create_mass_matrix (const hp::DoFHandler &dof, + void create_mass_matrix (const hp::DoFHandler &dof, const hp::QCollection &q, SparseMatrix &matrix, const Function &rhs, @@ -362,7 +360,7 @@ class MatrixCreator * cell-dependent shape functions. */ template - static + void create_boundary_mass_matrix (const Mapping &mapping, const DoFHandler &dof, const Quadrature &q, @@ -375,7 +373,7 @@ class MatrixCreator // * Same function, but for 1d. // */ -// static +// // void create_boundary_mass_matrix (const Mapping<1,1> &mapping, // const DoFHandler<1,1> &dof, // const Quadrature<0> &q, @@ -385,7 +383,7 @@ class MatrixCreator // std::vector&dof_to_boundary_mapping, // const Function<1> * const a = 0); // //codimension 1 -// static +// // void create_boundary_mass_matrix (const Mapping<1,2> &mapping, // const DoFHandler<1,2> &dof, // const Quadrature<0> &q, @@ -404,7 +402,7 @@ class MatrixCreator * mapping=MappingQ1@(). */ template - static + void create_boundary_mass_matrix (const DoFHandler &dof, const Quadrature &q, SparseMatrix &matrix, @@ -419,7 +417,7 @@ class MatrixCreator * objects. */ template - static + void create_boundary_mass_matrix (const hp::MappingCollection &mapping, const hp::DoFHandler &dof, const hp::QCollection &q, @@ -434,7 +432,7 @@ class MatrixCreator * Same function as above, but for hp * objects. */ -// static +// // void create_boundary_mass_matrix (const hp::MappingCollection<1,1> &mapping, // const hp::DoFHandler<1,1> &dof, // const hp::QCollection<0> &q, @@ -449,7 +447,7 @@ class MatrixCreator * objects. */ template - static + void create_boundary_mass_matrix (const hp::DoFHandler &dof, const hp::QCollection &q, SparseMatrix &matrix, @@ -473,7 +471,7 @@ class MatrixCreator * class for more information. */ template - static void create_laplace_matrix (const Mapping &mapping, + void create_laplace_matrix (const Mapping &mapping, const DoFHandler &dof, const Quadrature &q, SparseMatrix &matrix, @@ -486,7 +484,7 @@ class MatrixCreator * mapping=MappingQ1@(). */ template - static void create_laplace_matrix (const DoFHandler &dof, + void create_laplace_matrix (const DoFHandler &dof, const Quadrature &q, SparseMatrix &matrix, const Function * const a = 0); @@ -506,7 +504,7 @@ class MatrixCreator * class for more information. */ template - static void create_laplace_matrix (const Mapping &mapping, + void create_laplace_matrix (const Mapping &mapping, const DoFHandler &dof, const Quadrature &q, SparseMatrix &matrix, @@ -521,7 +519,7 @@ class MatrixCreator * mapping=MappingQ1@(). */ template - static void create_laplace_matrix (const DoFHandler &dof, + void create_laplace_matrix (const DoFHandler &dof, const Quadrature &q, SparseMatrix &matrix, const Function &rhs, @@ -534,7 +532,7 @@ class MatrixCreator * collections. */ template - static void create_laplace_matrix (const hp::MappingCollection &mapping, + void create_laplace_matrix (const hp::MappingCollection &mapping, const hp::DoFHandler &dof, const hp::QCollection &q, SparseMatrix &matrix, @@ -546,7 +544,7 @@ class MatrixCreator * collections. */ template - static void create_laplace_matrix (const hp::DoFHandler &dof, + void create_laplace_matrix (const hp::DoFHandler &dof, const hp::QCollection &q, SparseMatrix &matrix, const Function * const a = 0); @@ -557,7 +555,7 @@ class MatrixCreator * collections. */ template - static void create_laplace_matrix (const hp::MappingCollection &mapping, + void create_laplace_matrix (const hp::MappingCollection &mapping, const hp::DoFHandler &dof, const hp::QCollection &q, SparseMatrix &matrix, @@ -571,7 +569,7 @@ class MatrixCreator * collections. */ template - static void create_laplace_matrix (const hp::DoFHandler &dof, + void create_laplace_matrix (const hp::DoFHandler &dof, const hp::QCollection &q, SparseMatrix &matrix, const Function &rhs, @@ -582,93 +580,7 @@ class MatrixCreator * Exception */ DeclException0 (ExcComponentMismatch); - private: - /** - * Convenience abbreviation for - * pairs of DoF handler cell - * iterators. This type works - * just like a - * std::pair - * but is templatized on the - * dof handler that shouls be used. - */ - template - struct IteratorRange - { - /** - * Typedef for the iterator type. - */ - typedef typename DH::active_cell_iterator active_cell_iterator; - - /** - * Abbreviation for a pair of - * iterators. - */ - typedef std::pair iterator_pair; - - /** - * Constructor. Initialize - * the two values by the - * given values. - */ - IteratorRange (const active_cell_iterator &first, - const active_cell_iterator &second); - - /** - * Constructor taking a pair - * of values for - * initialization. - */ - IteratorRange (const iterator_pair &ip); - - /** - * Pair of iterators denoting - * a half-open range. - */ - active_cell_iterator first, second; - }; - - /** - * 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. - */ - template - static - void create_boundary_mass_matrix_1 (std_cxx1x::tuple &, - const DoFHandler &, - const Quadrature & > commons, - SparseMatrix &matrix, - const typename FunctionMap::type &boundary_functions, - Vector &rhs_vector, - std::vector&dof_to_boundary_mapping, - const Function * const coefficient, - const std::vector& component_mapping, - const IteratorRange > range, - Threads::ThreadMutex &mutex); - - /** - * Same as above, but for hp. - */ - template - static - void create_boundary_mass_matrix_1 (std_cxx1x::tuple &, - const hp::DoFHandler&, - const hp::QCollection & > commons, - SparseMatrix &matrix, - const typename FunctionMap::type &boundary_functions, - Vector &rhs_vector, - std::vector&dof_to_boundary_mapping, - const Function * const coefficient, - const std::vector& component_mapping, - const IteratorRange > range, - Threads::ThreadMutex &mutex); -}; +} @@ -873,8 +785,9 @@ class MatrixCreator * @ingroup numerics * @author Wolfgang Bangerth, 1998, 2000, 2004, 2005 */ -class MatrixTools : public MatrixCreator +class MatrixTools { +// using namespace MatrixCreator public: /** * Apply dirichlet boundary conditions diff --git a/deal.II/source/numerics/matrices.cc b/deal.II/source/numerics/matrices.cc index 3d7a1e2dca..5e5ae0962f 100644 --- a/deal.II/source/numerics/matrices.cc +++ b/deal.II/source/numerics/matrices.cc @@ -62,10 +62,80 @@ DEAL_II_NAMESPACE_OPEN -namespace internal +namespace MatrixCreator { - namespace MatrixCreator + namespace internal { + /** + * Convenience abbreviation for + * pairs of DoF handler cell + * iterators. This type works + * just like a + * std::pair + * but is templatized on the + * dof handler that shouls be used. + */ + template + struct IteratorRange + { + /** + * Typedef for the iterator type. + */ + typedef typename DH::active_cell_iterator active_cell_iterator; + + /** + * Abbreviation for a pair of + * iterators. + */ + typedef std::pair iterator_pair; + + /** + * Constructor. Initialize + * the two values by the + * given values. + */ + IteratorRange (const active_cell_iterator &first, + const active_cell_iterator &second); + + /** + * Constructor taking a pair + * of values for + * initialization. + */ + IteratorRange (const iterator_pair &ip); + + /** + * Pair of iterators denoting + * a half-open range. + */ + active_cell_iterator first, second; + }; + + + + + template + inline + IteratorRange:: + IteratorRange (const active_cell_iterator &first, + const active_cell_iterator &second) + : + first (first), + second (second) + {} + + + + template + inline + IteratorRange::IteratorRange (const iterator_pair &ip) + : + first (ip.first), + second (ip.second) + {} + + + namespace AssemblerData { template void mass_assembler (const CellIterator &cell, - internal::MatrixCreator::AssemblerData::Scratch &data, - internal::MatrixCreator::AssemblerData::CopyData ©_data) + MatrixCreator::internal::AssemblerData::Scratch &data, + MatrixCreator::internal::AssemblerData::CopyData ©_data) { data.x_fe_values.reinit (cell); const FEValues &fe_values = data.x_fe_values.get_present_fe_values (); @@ -329,8 +399,8 @@ namespace internal int spacedim, typename CellIterator> void laplace_assembler (const CellIterator &cell, - internal::MatrixCreator::AssemblerData::Scratch &data, - internal::MatrixCreator::AssemblerData::CopyData ©_data) + MatrixCreator::internal::AssemblerData::Scratch &data, + MatrixCreator::internal::AssemblerData::CopyData ©_data) { data.x_fe_values.reinit (cell); const FEValues &fe_values = data.x_fe_values.get_present_fe_values (); @@ -532,699 +602,573 @@ namespace internal } - -template -inline -MatrixCreator::IteratorRange:: -IteratorRange (const active_cell_iterator &first, - const active_cell_iterator &second) - : - first (first), - second (second) -{} - - - -template -inline -MatrixCreator::IteratorRange::IteratorRange (const iterator_pair &ip) - : - first (ip.first), - second (ip.second) -{} - - - - -template -void MatrixCreator::create_mass_matrix (const Mapping &mapping, - const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - const Function * const coefficient) +namespace MatrixCreator { - Assert (matrix.m() == dof.n_dofs(), - ExcDimensionMismatch (matrix.m(), dof.n_dofs())); - Assert (matrix.n() == dof.n_dofs(), - ExcDimensionMismatch (matrix.n(), dof.n_dofs())); - - hp::FECollection fe_collection (dof.get_fe()); - hp::QCollection q_collection (q); - hp::MappingCollection mapping_collection (mapping); - internal::MatrixCreator::AssemblerData::Scratch - assembler_data (fe_collection, - update_values | update_JxW_values | - (coefficient != 0 ? update_quadrature_points : UpdateFlags(0)), - coefficient, /*rhs_function=*/0, - q_collection, mapping_collection); - - internal::MatrixCreator::AssemblerData::CopyData copy_data; - copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(), - assembler_data.fe_collection.max_dofs_per_cell()); - copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell()); - copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell()); - - WorkStream::run (dof.begin_active(), - static_cast::active_cell_iterator>(dof.end()), - &internal::MatrixCreator::mass_assembler::active_cell_iterator>, - std_cxx1x::bind (&internal::MatrixCreator:: - copy_local_to_global, Vector >, - std_cxx1x::_1, &matrix, (Vector*)0), - assembler_data, - copy_data); -} + template + void create_mass_matrix (const Mapping &mapping, + const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const Function * const coefficient) + { + Assert (matrix.m() == dof.n_dofs(), + ExcDimensionMismatch (matrix.m(), dof.n_dofs())); + Assert (matrix.n() == dof.n_dofs(), + ExcDimensionMismatch (matrix.n(), dof.n_dofs())); + + hp::FECollection fe_collection (dof.get_fe()); + hp::QCollection q_collection (q); + hp::MappingCollection mapping_collection (mapping); + MatrixCreator::internal::AssemblerData::Scratch + assembler_data (fe_collection, + update_values | update_JxW_values | + (coefficient != 0 ? update_quadrature_points : UpdateFlags(0)), + coefficient, /*rhs_function=*/0, + q_collection, mapping_collection); + + MatrixCreator::internal::AssemblerData::CopyData copy_data; + copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(), + assembler_data.fe_collection.max_dofs_per_cell()); + copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell()); + copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell()); + + WorkStream::run (dof.begin_active(), + static_cast::active_cell_iterator>(dof.end()), + &MatrixCreator::internal::mass_assembler::active_cell_iterator>, + std_cxx1x::bind (&MatrixCreator::internal:: + copy_local_to_global, Vector >, + std_cxx1x::_1, &matrix, (Vector*)0), + assembler_data, + copy_data); + } -template -void MatrixCreator::create_mass_matrix (const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - const Function * const coefficient) -{ - Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); - create_mass_matrix(StaticMappingQ1::mapping, dof, - q, matrix, coefficient); -} + template + void create_mass_matrix (const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const Function * const coefficient) + { + Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); + create_mass_matrix(StaticMappingQ1::mapping, dof, + q, matrix, coefficient); + } -template -void MatrixCreator::create_mass_matrix (const Mapping &mapping, - const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - const Function &rhs, - Vector &rhs_vector, - const Function * const coefficient) -{ - Assert (matrix.m() == dof.n_dofs(), - ExcDimensionMismatch (matrix.m(), dof.n_dofs())); - Assert (matrix.n() == dof.n_dofs(), - ExcDimensionMismatch (matrix.n(), dof.n_dofs())); - - hp::FECollection fe_collection (dof.get_fe()); - hp::QCollection q_collection (q); - hp::MappingCollection mapping_collection (mapping); - internal::MatrixCreator::AssemblerData::Scratch - assembler_data (fe_collection, - update_values | - update_JxW_values | update_quadrature_points, - coefficient, &rhs, - q_collection, mapping_collection); - internal::MatrixCreator::AssemblerData::CopyData copy_data; - copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(), - assembler_data.fe_collection.max_dofs_per_cell()); - copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell()); - copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell()); - - WorkStream::run (dof.begin_active(), - static_cast::active_cell_iterator>(dof.end()), - &internal::MatrixCreator::mass_assembler::active_cell_iterator>, - std_cxx1x::bind(&internal::MatrixCreator:: - copy_local_to_global, Vector >, - std_cxx1x::_1, &matrix, &rhs_vector), - assembler_data, - copy_data); -} + template + void create_mass_matrix (const Mapping &mapping, + const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const Function &rhs, + Vector &rhs_vector, + const Function * const coefficient) + { + Assert (matrix.m() == dof.n_dofs(), + ExcDimensionMismatch (matrix.m(), dof.n_dofs())); + Assert (matrix.n() == dof.n_dofs(), + ExcDimensionMismatch (matrix.n(), dof.n_dofs())); + + hp::FECollection fe_collection (dof.get_fe()); + hp::QCollection q_collection (q); + hp::MappingCollection mapping_collection (mapping); + MatrixCreator::internal::AssemblerData::Scratch + assembler_data (fe_collection, + update_values | + update_JxW_values | update_quadrature_points, + coefficient, &rhs, + q_collection, mapping_collection); + MatrixCreator::internal::AssemblerData::CopyData copy_data; + copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(), + assembler_data.fe_collection.max_dofs_per_cell()); + copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell()); + copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell()); + + WorkStream::run (dof.begin_active(), + static_cast::active_cell_iterator>(dof.end()), + &MatrixCreator::internal::mass_assembler::active_cell_iterator>, + std_cxx1x::bind(&MatrixCreator::internal:: + copy_local_to_global, Vector >, + std_cxx1x::_1, &matrix, &rhs_vector), + assembler_data, + copy_data); + } -template -void MatrixCreator::create_mass_matrix (const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - const Function &rhs, - Vector &rhs_vector, - const Function * const coefficient) -{ - Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); - create_mass_matrix(StaticMappingQ1::mapping, - dof, q, matrix, rhs, rhs_vector, coefficient); -} + template + void create_mass_matrix (const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const Function &rhs, + Vector &rhs_vector, + const Function * const coefficient) + { + Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); + create_mass_matrix(StaticMappingQ1::mapping, + dof, q, matrix, rhs, rhs_vector, coefficient); + } -template -void MatrixCreator::create_mass_matrix (const hp::MappingCollection &mapping, - const hp::DoFHandler &dof, - const hp::QCollection &q, - SparseMatrix &matrix, - const Function * const coefficient) -{ - Assert (matrix.m() == dof.n_dofs(), - ExcDimensionMismatch (matrix.m(), dof.n_dofs())); - Assert (matrix.n() == dof.n_dofs(), - ExcDimensionMismatch (matrix.n(), dof.n_dofs())); - - internal::MatrixCreator::AssemblerData::Scratch - assembler_data (dof.get_fe(), - update_values | update_JxW_values | - (coefficient != 0 ? update_quadrature_points : UpdateFlags(0)), - coefficient, /*rhs_function=*/0, - q, mapping); - internal::MatrixCreator::AssemblerData::CopyData copy_data; - copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(), - assembler_data.fe_collection.max_dofs_per_cell()); - copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell()); - copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell()); - - WorkStream::run (dof.begin_active(), - static_cast::active_cell_iterator>(dof.end()), - &internal::MatrixCreator::mass_assembler::active_cell_iterator>, - std_cxx1x::bind (&internal::MatrixCreator:: - copy_local_to_global, Vector >, - std_cxx1x::_1, &matrix, (Vector*)0), - assembler_data, - copy_data); -} + template + void create_mass_matrix (const hp::MappingCollection &mapping, + const hp::DoFHandler &dof, + const hp::QCollection &q, + SparseMatrix &matrix, + const Function * const coefficient) + { + Assert (matrix.m() == dof.n_dofs(), + ExcDimensionMismatch (matrix.m(), dof.n_dofs())); + Assert (matrix.n() == dof.n_dofs(), + ExcDimensionMismatch (matrix.n(), dof.n_dofs())); + + MatrixCreator::internal::AssemblerData::Scratch + assembler_data (dof.get_fe(), + update_values | update_JxW_values | + (coefficient != 0 ? update_quadrature_points : UpdateFlags(0)), + coefficient, /*rhs_function=*/0, + q, mapping); + MatrixCreator::internal::AssemblerData::CopyData copy_data; + copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(), + assembler_data.fe_collection.max_dofs_per_cell()); + copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell()); + copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell()); + + WorkStream::run (dof.begin_active(), + static_cast::active_cell_iterator>(dof.end()), + &MatrixCreator::internal::mass_assembler::active_cell_iterator>, + std_cxx1x::bind (&MatrixCreator::internal:: + copy_local_to_global, Vector >, + std_cxx1x::_1, &matrix, (Vector*)0), + assembler_data, + copy_data); + } -template -void MatrixCreator::create_mass_matrix (const hp::DoFHandler &dof, - const hp::QCollection &q, - SparseMatrix &matrix, - const Function * const coefficient) -{ - Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); - create_mass_matrix(hp::StaticMappingQ1::mapping_collection, dof, q, matrix, coefficient); -} + template + void create_mass_matrix (const hp::DoFHandler &dof, + const hp::QCollection &q, + SparseMatrix &matrix, + const Function * const coefficient) + { + Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); + create_mass_matrix(hp::StaticMappingQ1::mapping_collection, dof, q, matrix, coefficient); + } -template -void MatrixCreator::create_mass_matrix (const hp::MappingCollection &mapping, - const hp::DoFHandler &dof, - const hp::QCollection &q, - SparseMatrix &matrix, - const Function &rhs, - Vector &rhs_vector, - const Function * const coefficient) -{ - Assert (matrix.m() == dof.n_dofs(), - ExcDimensionMismatch (matrix.m(), dof.n_dofs())); - Assert (matrix.n() == dof.n_dofs(), - ExcDimensionMismatch (matrix.n(), dof.n_dofs())); - - internal::MatrixCreator::AssemblerData::Scratch - assembler_data (dof.get_fe(), - update_values | - update_JxW_values | update_quadrature_points, - coefficient, &rhs, - q, mapping); - internal::MatrixCreator::AssemblerData::CopyData copy_data; - copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(), - assembler_data.fe_collection.max_dofs_per_cell()); - copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell()); - copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell()); - - WorkStream::run (dof.begin_active(), - static_cast::active_cell_iterator>(dof.end()), - &internal::MatrixCreator::mass_assembler::active_cell_iterator>, - std_cxx1x::bind (&internal::MatrixCreator:: - copy_local_to_global, Vector >, - std_cxx1x::_1, &matrix, &rhs_vector), - assembler_data, - copy_data); -} + template + void create_mass_matrix (const hp::MappingCollection &mapping, + const hp::DoFHandler &dof, + const hp::QCollection &q, + SparseMatrix &matrix, + const Function &rhs, + Vector &rhs_vector, + const Function * const coefficient) + { + Assert (matrix.m() == dof.n_dofs(), + ExcDimensionMismatch (matrix.m(), dof.n_dofs())); + Assert (matrix.n() == dof.n_dofs(), + ExcDimensionMismatch (matrix.n(), dof.n_dofs())); + + MatrixCreator::internal::AssemblerData::Scratch + assembler_data (dof.get_fe(), + update_values | + update_JxW_values | update_quadrature_points, + coefficient, &rhs, + q, mapping); + MatrixCreator::internal::AssemblerData::CopyData copy_data; + copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(), + assembler_data.fe_collection.max_dofs_per_cell()); + copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell()); + copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell()); + + WorkStream::run (dof.begin_active(), + static_cast::active_cell_iterator>(dof.end()), + &MatrixCreator::internal::mass_assembler::active_cell_iterator>, + std_cxx1x::bind (&MatrixCreator::internal:: + copy_local_to_global, Vector >, + std_cxx1x::_1, &matrix, &rhs_vector), + assembler_data, + copy_data); + } -template -void MatrixCreator::create_mass_matrix (const hp::DoFHandler &dof, - const hp::QCollection &q, - SparseMatrix &matrix, - const Function &rhs, - Vector &rhs_vector, - const Function * const coefficient) -{ - Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); - create_mass_matrix(hp::StaticMappingQ1::mapping_collection, dof, q, - matrix, rhs, rhs_vector, coefficient); -} + template + void create_mass_matrix (const hp::DoFHandler &dof, + const hp::QCollection &q, + SparseMatrix &matrix, + const Function &rhs, + Vector &rhs_vector, + const Function * const coefficient) + { + Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); + create_mass_matrix(hp::StaticMappingQ1::mapping_collection, dof, q, + matrix, rhs, rhs_vector, coefficient); + } -template -void -MatrixCreator::create_boundary_mass_matrix (const Mapping &mapping, - const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - const typename FunctionMap::type &boundary_functions, - Vector &rhs_vector, - std::vector &dof_to_boundary_mapping, - const Function * const coefficient, - std::vector component_mapping) -{ - // what would that be in 1d? the - // identity matrix on the boundary - // dofs? - if (dim == 1) - { - Assert (false, ExcNotImplemented()); - return; - } - const FiniteElement &fe = dof.get_fe(); - const unsigned int n_components = fe.n_components(); - - Assert (matrix.n() == dof.n_boundary_dofs(boundary_functions), - ExcInternalError()); - Assert (matrix.n() == matrix.m(), ExcInternalError()); - Assert (matrix.n() == rhs_vector.size(), ExcInternalError()); - Assert (boundary_functions.size() != 0, ExcInternalError()); - Assert (dof_to_boundary_mapping.size() == dof.n_dofs(), - ExcInternalError()); - Assert (coefficient ==0 || - coefficient->n_components==1 || - coefficient->n_components==n_components, ExcComponentMismatch()); - - if (component_mapping.size() == 0) + namespace + { + template + void + create_boundary_mass_matrix_1 (std_cxx1x::tuple &, + const DoFHandler &, + const Quadrature & > commons, + SparseMatrix &matrix, + const typename FunctionMap::type &boundary_functions, + Vector &rhs_vector, + std::vector &dof_to_boundary_mapping, + const Function * const coefficient, + const std::vector& component_mapping, + const MatrixCreator::internal::IteratorRange > range, + Threads::ThreadMutex &mutex) { - AssertDimension (n_components, boundary_functions.begin()->second->n_components); - for (unsigned int i=0;i threads; - - // define starting and end point - // for each thread - typedef typename DoFHandler::active_cell_iterator active_cell_iterator; - std::vector > thread_ranges - = Threads::split_range (dof.begin_active(), - dof.end(), n_threads); - - // mutex to synchronise access to - // the matrix - Threads::ThreadMutex mutex; - - typedef std_cxx1x::tuple&, - const DoFHandler&, - const Quadrature&> Commons; - - // then assemble in parallel - typedef void (*create_boundary_mass_matrix_1_t) - (Commons, - SparseMatrix &matrix, - const typename FunctionMap::type &boundary_functions, - Vector &rhs_vector, - std::vector &dof_to_boundary_mapping, - const Function * const coefficient, - const std::vector& component_mapping, - const IteratorRange > range, - Threads::ThreadMutex &mutex); - create_boundary_mass_matrix_1_t p = &MatrixCreator::template create_boundary_mass_matrix_1; - -//TODO: Use WorkStream here - for (unsigned int thread=0; thread -void -MatrixCreator:: -create_boundary_mass_matrix_1<2,3> (std_cxx1x::tuple &, - const DoFHandler<2,3> &, - const Quadrature<1> & > , - SparseMatrix &, - const FunctionMap<3>::type &, - Vector &, - std::vector &, - const Function<3> * const , - const std::vector &, - const IteratorRange > , - Threads::ThreadMutex &) -{ - Assert(false,ExcNotImplemented()); -} - - - -template -void -MatrixCreator:: -create_boundary_mass_matrix_1 (std_cxx1x::tuple &, - const DoFHandler &, - const Quadrature & > commons, - SparseMatrix &matrix, - const typename FunctionMap::type &boundary_functions, - Vector &rhs_vector, - std::vector &dof_to_boundary_mapping, - const Function * const coefficient, - const std::vector& component_mapping, - const IteratorRange > range, - Threads::ThreadMutex &mutex) -{ - // All assertions for this function - // are in the calling function - // before creating threads. - const Mapping& mapping = std_cxx1x::get<0>(commons); - const DoFHandler& dof = std_cxx1x::get<1>(commons); - const Quadrature& q = std_cxx1x::get<2>(commons); - - const FiniteElement &fe = dof.get_fe(); - 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_primitive = fe.is_primitive(); - - const unsigned int dofs_per_cell = fe.dofs_per_cell, - dofs_per_face = fe.dofs_per_face; - - FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell); - Vector cell_vector(dofs_per_cell); - - - UpdateFlags update_flags = UpdateFlags (update_values | - update_JxW_values | - update_normal_vectors | - update_quadrature_points); - FEFaceValues fe_values (mapping, fe, q, update_flags); - - // two variables for the coefficient, - // one for the two cases indicated in - // the name - std::vector coefficient_values (fe_values.n_quadrature_points, 1.); - std::vector > coefficient_vector_values (fe_values.n_quadrature_points, - Vector(n_components)); - const bool coefficient_is_vector = (coefficient != 0 && coefficient->n_components != 1) - ? true : false; - - std::vector rhs_values_scalar (fe_values.n_quadrature_points); - std::vector > rhs_values_system (fe_values.n_quadrature_points, - Vector(n_function_components)); - - std::vector dofs (dofs_per_cell); - std::vector dofs_on_face_vector (dofs_per_face); - - // for each dof on the cell, have a - // flag whether it is on the face - std::vector dof_is_on_face(dofs_per_cell); - - typename DoFHandler::active_cell_iterator cell = range.first; - for (; cell!=range.second; ++cell) - for (unsigned int face=0; face::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; - - fe_values.reinit (cell, face); - - if (fe_is_system) - // FE has several components + // All assertions for this function + // are in the calling function + // before creating threads. + const Mapping& mapping = std_cxx1x::get<0>(commons); + const DoFHandler& dof = std_cxx1x::get<1>(commons); + const Quadrature& q = std_cxx1x::get<2>(commons); + + const FiniteElement &fe = dof.get_fe(); + 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_primitive = fe.is_primitive(); + + const unsigned int dofs_per_cell = fe.dofs_per_cell, + dofs_per_face = fe.dofs_per_face; + + FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell); + Vector cell_vector(dofs_per_cell); + + + UpdateFlags update_flags = UpdateFlags (update_values | + update_JxW_values | + update_normal_vectors | + update_quadrature_points); + FEFaceValues fe_values (mapping, fe, q, update_flags); + + // two variables for the coefficient, + // one for the two cases indicated in + // the name + std::vector coefficient_values (fe_values.n_quadrature_points, 1.); + std::vector > coefficient_vector_values (fe_values.n_quadrature_points, + Vector(n_components)); + const bool coefficient_is_vector = (coefficient != 0 && coefficient->n_components != 1) + ? true : false; + + std::vector rhs_values_scalar (fe_values.n_quadrature_points); + std::vector > rhs_values_system (fe_values.n_quadrature_points, + Vector(n_function_components)); + + std::vector dofs (dofs_per_cell); + std::vector dofs_on_face_vector (dofs_per_face); + + // for each dof on the cell, have a + // flag whether it is on the face + std::vector dof_is_on_face(dofs_per_cell); + + typename DoFHandler::active_cell_iterator cell = range.first; + for (; cell!=range.second; ++cell) + for (unsigned int face=0; face::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()) { - 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 > normal_adjustment(fe_values.n_quadrature_points, - std::vector(n_components, 1.)); + fe_values.reinit (cell, face); - for (unsigned int comp = 0;comp& base = fe.base_element(fe.component_to_base_index(comp).first); - const unsigned int bcomp = fe.component_to_base_index(comp).second; + 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::H1) && - base.conforms(FiniteElementData::Hdiv)) - for (unsigned int point=0; point > normal_adjustment(fe_values.n_quadrature_points, + std::vector(n_components, 1.)); - for (unsigned int point=0; point& base = fe.base_element(fe.component_to_base_index(comp).first); + const unsigned int bcomp = fe.component_to_base_index(comp).second; + + if (!base.conforms(FiniteElementData::H1) && + base.conforms(FiniteElementData::Hdiv)) + for (unsigned int point=0; pointface(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; pointface(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; pointget_dof_indices (dofs); - cell->face(face)->get_dof_indices (dofs_on_face_vector); - for (unsigned int i=0; iget_dof_indices (dofs); + cell->face(face)->get_dof_indices (dofs_on_face_vector); + for (unsigned int i=0; i -void MatrixCreator::create_boundary_mass_matrix (const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - const typename FunctionMap::type &rhs, - Vector &rhs_vector, - std::vector &dof_to_boundary_mapping, - const Function * const a, - std::vector component_mapping) -{ - Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); - create_boundary_mass_matrix(StaticMappingQ1::mapping, dof, q, - matrix,rhs, rhs_vector, dof_to_boundary_mapping, a, component_mapping); -} + template <> + void + create_boundary_mass_matrix_1<2,3> (std_cxx1x::tuple &, + const DoFHandler<2,3> &, + const Quadrature<1> & > , + SparseMatrix &, + const FunctionMap<3>::type &, + Vector &, + std::vector &, + const Function<3> * const , + const std::vector &, + const MatrixCreator::internal::IteratorRange > , + Threads::ThreadMutex &) + { + Assert(false,ExcNotImplemented()); + } + } -template -void -MatrixCreator::create_boundary_mass_matrix (const hp::MappingCollection &mapping, - const hp::DoFHandler &dof, - const hp::QCollection &q, - SparseMatrix &matrix, - const typename FunctionMap::type &boundary_functions, - Vector &rhs_vector, - std::vector &dof_to_boundary_mapping, - const Function * const coefficient, - std::vector component_mapping) -{ - // what would that be in 1d? the - // identity matrix on the boundary - // dofs? - if (dim == 1) - { - Assert (false, ExcNotImplemented()); - return; - } + template + void + create_boundary_mass_matrix (const Mapping &mapping, + const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const typename FunctionMap::type &boundary_functions, + Vector &rhs_vector, + std::vector &dof_to_boundary_mapping, + const Function * const coefficient, + std::vector component_mapping) + { + // what would that be in 1d? the + // identity matrix on the boundary + // dofs? + if (dim == 1) + { + Assert (false, ExcNotImplemented()); + return; + } - const hp::FECollection &fe_collection = dof.get_fe(); - const unsigned int n_components = fe_collection.n_components(); - - Assert (matrix.n() == dof.n_boundary_dofs(boundary_functions), - ExcInternalError()); - Assert (matrix.n() == matrix.m(), ExcInternalError()); - Assert (matrix.n() == rhs_vector.size(), ExcInternalError()); - Assert (boundary_functions.size() != 0, ExcInternalError()); - Assert (dof_to_boundary_mapping.size() == dof.n_dofs(), - ExcInternalError()); - Assert (coefficient ==0 || - coefficient->n_components==1 || - coefficient->n_components==n_components, ExcComponentMismatch()); - - if (component_mapping.size() == 0) - { - AssertDimension (n_components, boundary_functions.begin()->second->n_components); - for (unsigned int i=0;i &fe = dof.get_fe(); + const unsigned int n_components = fe.n_components(); + + Assert (matrix.n() == dof.n_boundary_dofs(boundary_functions), + ExcInternalError()); + Assert (matrix.n() == matrix.m(), ExcInternalError()); + Assert (matrix.n() == rhs_vector.size(), ExcInternalError()); + Assert (boundary_functions.size() != 0, ExcInternalError()); + Assert (dof_to_boundary_mapping.size() == dof.n_dofs(), + ExcInternalError()); + Assert (coefficient ==0 || + coefficient->n_components==1 || + coefficient->n_components==n_components, ExcComponentMismatch()); + + if (component_mapping.size() == 0) + { + AssertDimension (n_components, boundary_functions.begin()->second->n_components); + for (unsigned int i=0;i threads; + const unsigned int n_threads = multithread_info.n_default_threads; + Threads::ThreadGroup<> threads; - // define starting and end point - // for each thread - typedef typename hp::DoFHandler::active_cell_iterator active_cell_iterator; - std::vector > thread_ranges - = Threads::split_range (dof.begin_active(), - dof.end(), n_threads); + // define starting and end point + // for each thread + typedef typename DoFHandler::active_cell_iterator active_cell_iterator; + std::vector > thread_ranges + = Threads::split_range (dof.begin_active(), + dof.end(), n_threads); - typedef std_cxx1x::tuple&, - const hp::DoFHandler&, - const hp::QCollection&> Commons; + // mutex to synchronise access to + // the matrix + Threads::ThreadMutex mutex; - // mutex to synchronise access to - // the matrix - Threads::ThreadMutex mutex; + typedef std_cxx1x::tuple&, + const DoFHandler&, + const Quadrature&> Commons; - // then assemble in parallel - typedef void (*create_boundary_mass_matrix_1_t) + // then assemble in parallel + typedef void (*create_boundary_mass_matrix_1_t) (Commons, SparseMatrix &matrix, const typename FunctionMap::type &boundary_functions, @@ -1232,635 +1176,747 @@ MatrixCreator::create_boundary_mass_matrix (const hp::MappingCollection &dof_to_boundary_mapping, const Function * const coefficient, const std::vector& component_mapping, - const IteratorRange > range, + const MatrixCreator::internal::IteratorRange > range, Threads::ThreadMutex &mutex); - create_boundary_mass_matrix_1_t p = &MatrixCreator::template create_boundary_mass_matrix_1; + create_boundary_mass_matrix_1_t p + = &create_boundary_mass_matrix_1; //TODO: Use WorkStream here - for (unsigned int thread=0; thread -void -MatrixCreator:: -create_boundary_mass_matrix_1 (std_cxx1x::tuple &, - const hp::DoFHandler &, - const hp::QCollection &> commons, - SparseMatrix &matrix, - const typename FunctionMap::type &boundary_functions, - Vector &rhs_vector, - std::vector &dof_to_boundary_mapping, - const Function * const coefficient, - const std::vector& component_mapping, - const IteratorRange > range, - Threads::ThreadMutex &mutex) -{ - const hp::MappingCollection& mapping = std_cxx1x::get<0>(commons); - const hp::DoFHandler& dof = std_cxx1x::get<1>(commons); - const hp::QCollection& q = std_cxx1x::get<2>(commons); - const hp::FECollection &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) + namespace + { + + template + void + create_boundary_mass_matrix_1 (std_cxx1x::tuple &, + const hp::DoFHandler &, + const hp::QCollection &> commons, + SparseMatrix &matrix, + const typename FunctionMap::type &boundary_functions, + Vector &rhs_vector, + std::vector &dof_to_boundary_mapping, + const Function * const coefficient, + const std::vector& component_mapping, + const MatrixCreator::internal::IteratorRange > range, + Threads::ThreadMutex &mutex) { - unsigned int max_element = 0; - for (std::vector::const_iterator i=dof_to_boundary_mapping.begin(); - i!=dof_to_boundary_mapping.end(); ++i) - if ((*i != hp::DoFHandler::invalid_dof_index) && - (*i > max_element)) - max_element = *i; - Assert (max_element == matrix.n()-1, ExcInternalError()); - }; + const hp::MappingCollection& mapping = std_cxx1x::get<0>(commons); + const hp::DoFHandler& dof = std_cxx1x::get<1>(commons); + const hp::QCollection& q = std_cxx1x::get<2>(commons); + const hp::FECollection &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) + { + unsigned int max_element = 0; + for (std::vector::const_iterator i=dof_to_boundary_mapping.begin(); + i!=dof_to_boundary_mapping.end(); ++i) + if ((*i != hp::DoFHandler::invalid_dof_index) && + (*i > max_element)) + max_element = *i; + Assert (max_element == matrix.n()-1, ExcInternalError()); + }; #endif - const unsigned int max_dofs_per_cell = fe_collection.max_dofs_per_cell(), - max_dofs_per_face = fe_collection.max_dofs_per_face(); + const unsigned int max_dofs_per_cell = fe_collection.max_dofs_per_cell(), + max_dofs_per_face = fe_collection.max_dofs_per_face(); - FullMatrix cell_matrix(max_dofs_per_cell, max_dofs_per_cell); - Vector cell_vector(max_dofs_per_cell); + FullMatrix cell_matrix(max_dofs_per_cell, max_dofs_per_cell); + Vector cell_vector(max_dofs_per_cell); - UpdateFlags update_flags = UpdateFlags (update_values | - update_JxW_values | - update_quadrature_points); - hp::FEFaceValues x_fe_values (mapping, fe_collection, q, update_flags); + UpdateFlags update_flags = UpdateFlags (update_values | + update_JxW_values | + update_quadrature_points); + hp::FEFaceValues x_fe_values (mapping, fe_collection, q, update_flags); - // two variables for the coefficient, - // one for the two cases indicated in - // the name - std::vector coefficient_values; - std::vector > coefficient_vector_values; + // two variables for the coefficient, + // one for the two cases indicated in + // the name + std::vector coefficient_values; + std::vector > coefficient_vector_values; - std::vector rhs_values_scalar; - std::vector > rhs_values_system; + std::vector rhs_values_scalar; + std::vector > rhs_values_system; - std::vector dofs (max_dofs_per_cell); - std::vector dofs_on_face_vector (max_dofs_per_face); + std::vector dofs (max_dofs_per_cell); + std::vector dofs_on_face_vector (max_dofs_per_face); - // for each dof on the cell, have a - // flag whether it is on the face - std::vector dof_is_on_face(max_dofs_per_cell); + // for each dof on the cell, have a + // flag whether it is on the face + std::vector dof_is_on_face(max_dofs_per_cell); - typename hp::DoFHandler::active_cell_iterator cell = range.first; - for (; cell!=range.second; ++cell) - for (unsigned int face=0; face::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 &fe_values = x_fe_values.get_present_fe_values (); + typename hp::DoFHandler::active_cell_iterator cell = range.first; + for (; cell!=range.second; ++cell) + for (unsigned int face=0; face::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 &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 &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; + const FiniteElement &fe = cell->get_fe(); + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int dofs_per_face = fe.dofs_per_face; - if (fe_is_system) - // FE has several components - { - rhs_values_system.resize (fe_values.n_quadrature_points, - Vector(n_function_components)); - boundary_functions.find(cell->face(face)->boundary_indicator()) - ->second->vector_value_list (fe_values.get_quadrature_points(), - rhs_values_system); + cell_matrix.reinit (dofs_per_cell, dofs_per_cell); + cell_vector.reinit (dofs_per_cell); + cell_matrix = 0; + cell_vector = 0; - if (coefficient != 0) + if (fe_is_system) + // FE has several components { - if (coefficient->n_components==1) + rhs_values_system.resize (fe_values.n_quadrature_points, + Vector(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) { - coefficient_values.resize (fe_values.n_quadrature_points); - coefficient->value_list (fe_values.get_quadrature_points(), - coefficient_values); - for (unsigned int point=0; pointn_components==1) { - const double weight = fe_values.JxW(point); - for (unsigned int i=0; ivalue_list (fe_values.get_quadrature_points(), + coefficient_values); + for (unsigned int point=0; point(n_components)); - coefficient->vector_value_list (fe_values.get_quadrature_points(), - coefficient_vector_values); - for (unsigned int point=0; point(n_components)); + coefficient->vector_value_list (fe_values.get_quadrature_points(), + coefficient_vector_values); + for (unsigned int point=0; pointface(face)->boundary_indicator()) - ->second->value_list (fe_values.get_quadrature_points(), rhs_values_scalar); - - if (coefficient != 0) + } + else + // FE is a scalar one { - coefficient_values.resize (fe_values.n_quadrature_points); - coefficient->value_list (fe_values.get_quadrature_points(), - coefficient_values); - for (unsigned int point=0; pointface(face)->boundary_indicator()) + ->second->value_list (fe_values.get_quadrature_points(), rhs_values_scalar); + + if (coefficient != 0) { - const double weight = fe_values.JxW(point); - for (unsigned int i=0; ivalue_list (fe_values.get_quadrature_points(), + coefficient_values); + for (unsigned int point=0; pointget_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; iget_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 max_diag_entry) - max_diag_entry = std::fabs(cell_matrix(i,i)); + double max_diag_entry = 0; + for (unsigned int i=0; i max_diag_entry) + max_diag_entry = std::fabs(cell_matrix(i,i)); #endif - // lock the matrix - Threads::ThreadMutex::ScopedLock lock (mutex); - for (unsigned int i=0; i -void MatrixCreator::create_boundary_mass_matrix (const hp::DoFHandler &dof, - const hp::QCollection &q, - SparseMatrix &matrix, - const typename FunctionMap::type &rhs, - Vector &rhs_vector, - std::vector &dof_to_boundary_mapping, - const Function * const a, - std::vector component_mapping) -{ - Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); - create_boundary_mass_matrix(hp::StaticMappingQ1::mapping_collection, dof, q, - matrix,rhs, rhs_vector, dof_to_boundary_mapping, a, component_mapping); -} + template + void create_boundary_mass_matrix (const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const typename FunctionMap::type &rhs, + Vector &rhs_vector, + std::vector &dof_to_boundary_mapping, + const Function * const a, + std::vector component_mapping) + { + Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); + create_boundary_mass_matrix(StaticMappingQ1::mapping, dof, q, + matrix,rhs, rhs_vector, dof_to_boundary_mapping, a, component_mapping); + } -template -void MatrixCreator::create_laplace_matrix (const Mapping &mapping, - const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - const Function * const coefficient) -{ - Assert (matrix.m() == dof.n_dofs(), - ExcDimensionMismatch (matrix.m(), dof.n_dofs())); - Assert (matrix.n() == dof.n_dofs(), - ExcDimensionMismatch (matrix.n(), dof.n_dofs())); - - hp::FECollection fe_collection (dof.get_fe()); - hp::QCollection q_collection (q); - hp::MappingCollection mapping_collection (mapping); - internal::MatrixCreator::AssemblerData::Scratch - assembler_data (fe_collection, - update_gradients | update_JxW_values | - (coefficient != 0 ? update_quadrature_points : UpdateFlags(0)), - coefficient, /*rhs_function=*/0, - q_collection, mapping_collection); - internal::MatrixCreator::AssemblerData::CopyData copy_data; - copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(), - assembler_data.fe_collection.max_dofs_per_cell()); - copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell()); - copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell()); - - void (*copy_local_to_global) (const internal::MatrixCreator::AssemblerData::CopyData &, - SparseMatrix *, - Vector *) - = &internal::MatrixCreator:: - copy_local_to_global, Vector >; - - WorkStream::run (dof.begin_active(), - static_cast::active_cell_iterator>(dof.end()), - &internal::MatrixCreator::laplace_assembler::active_cell_iterator>, - std_cxx1x::bind (copy_local_to_global, - std_cxx1x::_1, &matrix, (Vector*)0), - assembler_data, - copy_data); -} + template + void + create_boundary_mass_matrix (const hp::MappingCollection &mapping, + const hp::DoFHandler &dof, + const hp::QCollection &q, + SparseMatrix &matrix, + const typename FunctionMap::type &boundary_functions, + Vector &rhs_vector, + std::vector &dof_to_boundary_mapping, + const Function * const coefficient, + std::vector component_mapping) + { + // what would that be in 1d? the + // identity matrix on the boundary + // dofs? + if (dim == 1) + { + Assert (false, ExcNotImplemented()); + return; + } + const hp::FECollection &fe_collection = dof.get_fe(); + const unsigned int n_components = fe_collection.n_components(); + + Assert (matrix.n() == dof.n_boundary_dofs(boundary_functions), + ExcInternalError()); + Assert (matrix.n() == matrix.m(), ExcInternalError()); + Assert (matrix.n() == rhs_vector.size(), ExcInternalError()); + Assert (boundary_functions.size() != 0, ExcInternalError()); + Assert (dof_to_boundary_mapping.size() == dof.n_dofs(), + ExcInternalError()); + Assert (coefficient ==0 || + coefficient->n_components==1 || + coefficient->n_components==n_components, ExcComponentMismatch()); + + if (component_mapping.size() == 0) + { + AssertDimension (n_components, boundary_functions.begin()->second->n_components); + for (unsigned int i=0;i -void MatrixCreator::create_laplace_matrix (const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - const Function * const coefficient) -{ - Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); - create_laplace_matrix(StaticMappingQ1::mapping, dof, q, matrix, coefficient); -} + const unsigned int n_threads = multithread_info.n_default_threads; + Threads::ThreadGroup<> threads; + // define starting and end point + // for each thread + typedef typename hp::DoFHandler::active_cell_iterator active_cell_iterator; + std::vector > thread_ranges + = Threads::split_range (dof.begin_active(), + dof.end(), n_threads); + typedef std_cxx1x::tuple&, + const hp::DoFHandler&, + const hp::QCollection&> Commons; -template -void MatrixCreator::create_laplace_matrix (const Mapping &mapping, - const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - const Function &rhs, - Vector &rhs_vector, - const Function * const coefficient) -{ - Assert (matrix.m() == dof.n_dofs(), - ExcDimensionMismatch (matrix.m(), dof.n_dofs())); - Assert (matrix.n() == dof.n_dofs(), - ExcDimensionMismatch (matrix.n(), dof.n_dofs())); - - hp::FECollection fe_collection (dof.get_fe()); - hp::QCollection q_collection (q); - hp::MappingCollection mapping_collection (mapping); - internal::MatrixCreator::AssemblerData::Scratch - assembler_data (fe_collection, - update_gradients | update_values | - update_JxW_values | update_quadrature_points, - coefficient, &rhs, - q_collection, mapping_collection); - internal::MatrixCreator::AssemblerData::CopyData copy_data; - copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(), - assembler_data.fe_collection.max_dofs_per_cell()); - copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell()); - copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell()); - - void (*copy_local_to_global) (const internal::MatrixCreator::AssemblerData::CopyData &, - SparseMatrix *, - Vector *) - = &internal::MatrixCreator:: - copy_local_to_global, Vector >; - - WorkStream::run (dof.begin_active(), - static_cast::active_cell_iterator>(dof.end()), - &internal::MatrixCreator::laplace_assembler::active_cell_iterator>, - std_cxx1x::bind (copy_local_to_global, - std_cxx1x::_1, &matrix, &rhs_vector), - assembler_data, - copy_data); -} + // mutex to synchronise access to + // the matrix + Threads::ThreadMutex mutex; + // then assemble in parallel + typedef void (*create_boundary_mass_matrix_1_t) + (Commons, + SparseMatrix &matrix, + const typename FunctionMap::type &boundary_functions, + Vector &rhs_vector, + std::vector &dof_to_boundary_mapping, + const Function * const coefficient, + const std::vector& component_mapping, + const MatrixCreator::internal::IteratorRange > range, + Threads::ThreadMutex &mutex); + create_boundary_mass_matrix_1_t p = &create_boundary_mass_matrix_1; +//TODO: Use WorkStream here + for (unsigned int thread=0; thread -void MatrixCreator::create_laplace_matrix (const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - const Function &rhs, - Vector &rhs_vector, - const Function * const coefficient) -{ - Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); - create_laplace_matrix(StaticMappingQ1::mapping, dof, q, - matrix, rhs, rhs_vector, coefficient); -} -template -void MatrixCreator::create_laplace_matrix (const hp::MappingCollection &mapping, - const hp::DoFHandler &dof, - const hp::QCollection &q, - SparseMatrix &matrix, - const Function * const coefficient) -{ - Assert (matrix.m() == dof.n_dofs(), - ExcDimensionMismatch (matrix.m(), dof.n_dofs())); - Assert (matrix.n() == dof.n_dofs(), - ExcDimensionMismatch (matrix.n(), dof.n_dofs())); - - internal::MatrixCreator::AssemblerData::Scratch - assembler_data (dof.get_fe(), - update_gradients | update_JxW_values | - (coefficient != 0 ? update_quadrature_points : UpdateFlags(0)), - coefficient, /*rhs_function=*/0, - q, mapping); - internal::MatrixCreator::AssemblerData::CopyData copy_data; - copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(), - assembler_data.fe_collection.max_dofs_per_cell()); - copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell()); - copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell()); - - void (*copy_local_to_global) (const internal::MatrixCreator::AssemblerData::CopyData &, - SparseMatrix *, - Vector *) - = &internal::MatrixCreator:: - copy_local_to_global, Vector >; - - WorkStream::run (dof.begin_active(), - static_cast::active_cell_iterator>(dof.end()), - &internal::MatrixCreator::laplace_assembler::active_cell_iterator>, - std_cxx1x::bind (copy_local_to_global, - std_cxx1x::_1, &matrix, (Vector*)0), - assembler_data, - copy_data); -} + template + void create_boundary_mass_matrix (const hp::DoFHandler &dof, + const hp::QCollection &q, + SparseMatrix &matrix, + const typename FunctionMap::type &rhs, + Vector &rhs_vector, + std::vector &dof_to_boundary_mapping, + const Function * const a, + std::vector component_mapping) + { + Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); + create_boundary_mass_matrix(hp::StaticMappingQ1::mapping_collection, dof, q, + matrix,rhs, rhs_vector, dof_to_boundary_mapping, a, component_mapping); + } -template -void MatrixCreator::create_laplace_matrix (const hp::DoFHandler &dof, - const hp::QCollection &q, - SparseMatrix &matrix, - const Function * const coefficient) -{ - Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); - create_laplace_matrix(hp::StaticMappingQ1::mapping_collection, dof, q, matrix, coefficient); -} + template + void create_laplace_matrix (const Mapping &mapping, + const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const Function * const coefficient) + { + Assert (matrix.m() == dof.n_dofs(), + ExcDimensionMismatch (matrix.m(), dof.n_dofs())); + Assert (matrix.n() == dof.n_dofs(), + ExcDimensionMismatch (matrix.n(), dof.n_dofs())); + + hp::FECollection fe_collection (dof.get_fe()); + hp::QCollection q_collection (q); + hp::MappingCollection mapping_collection (mapping); + MatrixCreator::internal::AssemblerData::Scratch + assembler_data (fe_collection, + update_gradients | update_JxW_values | + (coefficient != 0 ? update_quadrature_points : UpdateFlags(0)), + coefficient, /*rhs_function=*/0, + q_collection, mapping_collection); + MatrixCreator::internal::AssemblerData::CopyData copy_data; + copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(), + assembler_data.fe_collection.max_dofs_per_cell()); + copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell()); + copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell()); + + void (*copy_local_to_global) (const MatrixCreator::internal::AssemblerData::CopyData &, + SparseMatrix *, + Vector *) + = &MatrixCreator::internal:: + copy_local_to_global, Vector >; + + WorkStream::run (dof.begin_active(), + static_cast::active_cell_iterator>(dof.end()), + &MatrixCreator::internal::laplace_assembler::active_cell_iterator>, + std_cxx1x::bind (copy_local_to_global, + std_cxx1x::_1, &matrix, (Vector*)0), + assembler_data, + copy_data); + } -template -void MatrixCreator::create_laplace_matrix (const hp::MappingCollection &mapping, - const hp::DoFHandler &dof, - const hp::QCollection &q, - SparseMatrix &matrix, - const Function &rhs, - Vector &rhs_vector, - const Function * const coefficient) -{ - Assert (matrix.m() == dof.n_dofs(), - ExcDimensionMismatch (matrix.m(), dof.n_dofs())); - Assert (matrix.n() == dof.n_dofs(), - ExcDimensionMismatch (matrix.n(), dof.n_dofs())); - - internal::MatrixCreator::AssemblerData::Scratch - assembler_data (dof.get_fe(), - update_gradients | update_values | - update_JxW_values | update_quadrature_points, - coefficient, &rhs, - q, mapping); - internal::MatrixCreator::AssemblerData::CopyData copy_data; - copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(), - assembler_data.fe_collection.max_dofs_per_cell()); - copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell()); - copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell()); - - void (*copy_local_to_global) (const internal::MatrixCreator::AssemblerData::CopyData &, - SparseMatrix *, - Vector *) - = &internal::MatrixCreator:: - copy_local_to_global, Vector >; - - WorkStream::run (dof.begin_active(), - static_cast::active_cell_iterator>(dof.end()), - &internal::MatrixCreator::laplace_assembler::active_cell_iterator>, - std_cxx1x::bind (copy_local_to_global, - std_cxx1x::_1, &matrix, &rhs_vector), - assembler_data, - copy_data); -} + template + void create_laplace_matrix (const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const Function * const coefficient) + { + Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); + create_laplace_matrix(StaticMappingQ1::mapping, dof, q, matrix, coefficient); + } -template -void MatrixCreator::create_laplace_matrix (const hp::DoFHandler &dof, - const hp::QCollection &q, - SparseMatrix &matrix, - const Function &rhs, - Vector &rhs_vector, - const Function * const coefficient) -{ - Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); - create_laplace_matrix(hp::StaticMappingQ1::mapping_collection, dof, q, - matrix, rhs, rhs_vector, coefficient); -} + template + void create_laplace_matrix (const Mapping &mapping, + const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const Function &rhs, + Vector &rhs_vector, + const Function * const coefficient) + { + Assert (matrix.m() == dof.n_dofs(), + ExcDimensionMismatch (matrix.m(), dof.n_dofs())); + Assert (matrix.n() == dof.n_dofs(), + ExcDimensionMismatch (matrix.n(), dof.n_dofs())); + + hp::FECollection fe_collection (dof.get_fe()); + hp::QCollection q_collection (q); + hp::MappingCollection mapping_collection (mapping); + MatrixCreator::internal::AssemblerData::Scratch + assembler_data (fe_collection, + update_gradients | update_values | + update_JxW_values | update_quadrature_points, + coefficient, &rhs, + q_collection, mapping_collection); + MatrixCreator::internal::AssemblerData::CopyData copy_data; + copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(), + assembler_data.fe_collection.max_dofs_per_cell()); + copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell()); + copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell()); + + void (*copy_local_to_global) (const MatrixCreator::internal::AssemblerData::CopyData &, + SparseMatrix *, + Vector *) + = &MatrixCreator::internal:: + copy_local_to_global, Vector >; + + WorkStream::run (dof.begin_active(), + static_cast::active_cell_iterator>(dof.end()), + &MatrixCreator::internal::laplace_assembler::active_cell_iterator>, + std_cxx1x::bind (copy_local_to_global, + std_cxx1x::_1, &matrix, &rhs_vector), + assembler_data, + copy_data); + } + + + + template + void create_laplace_matrix (const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const Function &rhs, + Vector &rhs_vector, + const Function * const coefficient) + { + Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); + create_laplace_matrix(StaticMappingQ1::mapping, dof, q, + matrix, rhs, rhs_vector, coefficient); + } + template + void create_laplace_matrix (const hp::MappingCollection &mapping, + const hp::DoFHandler &dof, + const hp::QCollection &q, + SparseMatrix &matrix, + const Function * const coefficient) + { + Assert (matrix.m() == dof.n_dofs(), + ExcDimensionMismatch (matrix.m(), dof.n_dofs())); + Assert (matrix.n() == dof.n_dofs(), + ExcDimensionMismatch (matrix.n(), dof.n_dofs())); + + MatrixCreator::internal::AssemblerData::Scratch + assembler_data (dof.get_fe(), + update_gradients | update_JxW_values | + (coefficient != 0 ? update_quadrature_points : UpdateFlags(0)), + coefficient, /*rhs_function=*/0, + q, mapping); + MatrixCreator::internal::AssemblerData::CopyData copy_data; + copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(), + assembler_data.fe_collection.max_dofs_per_cell()); + copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell()); + copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell()); + + void (*copy_local_to_global) (const MatrixCreator::internal::AssemblerData::CopyData &, + SparseMatrix *, + Vector *) + = &MatrixCreator::internal:: + copy_local_to_global, Vector >; + + WorkStream::run (dof.begin_active(), + static_cast::active_cell_iterator>(dof.end()), + &MatrixCreator::internal::laplace_assembler::active_cell_iterator>, + std_cxx1x::bind (copy_local_to_global, + std_cxx1x::_1, &matrix, (Vector*)0), + assembler_data, + copy_data); + } + + + + template + void create_laplace_matrix (const hp::DoFHandler &dof, + const hp::QCollection &q, + SparseMatrix &matrix, + const Function * const coefficient) + { + Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); + create_laplace_matrix(hp::StaticMappingQ1::mapping_collection, dof, q, matrix, coefficient); + } + + + + template + void create_laplace_matrix (const hp::MappingCollection &mapping, + const hp::DoFHandler &dof, + const hp::QCollection &q, + SparseMatrix &matrix, + const Function &rhs, + Vector &rhs_vector, + const Function * const coefficient) + { + Assert (matrix.m() == dof.n_dofs(), + ExcDimensionMismatch (matrix.m(), dof.n_dofs())); + Assert (matrix.n() == dof.n_dofs(), + ExcDimensionMismatch (matrix.n(), dof.n_dofs())); + + MatrixCreator::internal::AssemblerData::Scratch + assembler_data (dof.get_fe(), + update_gradients | update_values | + update_JxW_values | update_quadrature_points, + coefficient, &rhs, + q, mapping); + MatrixCreator::internal::AssemblerData::CopyData copy_data; + copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(), + assembler_data.fe_collection.max_dofs_per_cell()); + copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell()); + copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell()); + + void (*copy_local_to_global) (const MatrixCreator::internal::AssemblerData::CopyData &, + SparseMatrix *, + Vector *) + = &MatrixCreator::internal:: + copy_local_to_global, Vector >; + + WorkStream::run (dof.begin_active(), + static_cast::active_cell_iterator>(dof.end()), + &MatrixCreator::internal::laplace_assembler::active_cell_iterator>, + std_cxx1x::bind (copy_local_to_global, + std_cxx1x::_1, &matrix, &rhs_vector), + assembler_data, + copy_data); + } + + + + template + void create_laplace_matrix (const hp::DoFHandler &dof, + const hp::QCollection &q, + SparseMatrix &matrix, + const Function &rhs, + Vector &rhs_vector, + const Function * const coefficient) + { + Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); + create_laplace_matrix(hp::StaticMappingQ1::mapping_collection, dof, q, + matrix, rhs, rhs_vector, coefficient); + } + +} // namespace MatrixCreator + + //TODO:[WB] I don't think that the optimized storage of diagonals is needed (GK) template @@ -1901,7 +1957,7 @@ MatrixTools::apply_boundary_values (const std::map &boundar break; } - + std::map::const_iterator dof = boundary_values.begin(), endd = boundary_values.end(); const SparsityPattern &sparsity = matrix.get_sparsity_pattern(); @@ -1910,10 +1966,10 @@ MatrixTools::apply_boundary_values (const std::map &boundar for (; dof != endd; ++dof) { Assert (dof->first < n_dofs, ExcInternalError()); - + const unsigned int dof_number = dof->first; // for each boundary dof: - + // set entries of this line // to zero except for the diagonal // entry. Note that the diagonal @@ -1968,7 +2024,7 @@ MatrixTools::apply_boundary_values (const std::map &boundar // of this line for the Gauss // elimination step const number diagonal_entry = matrix.diag_element(dof_number); - + // we have to loop over all // rows of the matrix which // have a nonzero entry in @@ -2012,11 +2068,11 @@ MatrixTools::apply_boundary_values (const std::map &boundar const unsigned int global_entry = (p - &sparsity_colnums[sparsity_rowstart[0]]); - + // correct right hand side right_hand_side(row) -= matrix.global_entry(global_entry) / diagonal_entry * new_rhs; - + // set matrix entry to zero matrix.global_entry(global_entry) = 0.; } @@ -2038,14 +2094,14 @@ MatrixTools::apply_boundary_values (const std::map &boundar const bool eliminate_columns) { const unsigned int blocks = matrix.n_block_rows(); - + Assert (matrix.n() == right_hand_side.size(), ExcDimensionMismatch(matrix.n(), right_hand_side.size())); Assert (matrix.n() == solution.size(), ExcDimensionMismatch(matrix.n(), solution.size())); Assert (matrix.n_block_rows() == matrix.n_block_cols(), ExcNotQuadratic()); - Assert (matrix.get_sparsity_pattern().get_row_indices() == + Assert (matrix.get_sparsity_pattern().get_row_indices() == matrix.get_sparsity_pattern().get_column_indices(), ExcNotQuadratic()); Assert (matrix.get_sparsity_pattern().get_column_indices() == @@ -2059,7 +2115,7 @@ MatrixTools::apply_boundary_values (const std::map &boundar Assert (matrix.block(i,i).get_sparsity_pattern().optimize_diagonal(), SparsityPattern::ExcDiagonalNotOptimized()); - + // if no boundary values are to be applied // simply return if (boundary_values.size() == 0) @@ -2080,7 +2136,7 @@ MatrixTools::apply_boundary_values (const std::map &boundar for (unsigned int i=0; i &boundar // blocks? if so, use 1.0 instead if (first_nonzero_diagonal_entry == 0) first_nonzero_diagonal_entry = 1; - + std::map::const_iterator dof = boundary_values.begin(), endd = boundary_values.end(); const BlockSparsityPattern & sparsity_pattern = matrix.get_sparsity_pattern(); - + // pointer to the mapping between // global and block indices. since // the row and column mappings are @@ -2108,7 +2164,7 @@ MatrixTools::apply_boundary_values (const std::map &boundar // one of them const BlockIndices & index_mapping = sparsity_pattern.get_column_indices(); - + // now loop over all boundary dofs for (; dof != endd; ++dof) { @@ -2122,7 +2178,7 @@ MatrixTools::apply_boundary_values (const std::map &boundar block_index = index_mapping.global_to_local (dof_number); // for each boundary dof: - + // set entries of this line // to zero except for the diagonal // entry. Note that the diagonal @@ -2149,16 +2205,16 @@ MatrixTools::apply_boundary_values (const std::map &boundar // matrix so the diagonal // element is the first of // this row. - const unsigned int + const unsigned int last = local_sparsity.get_rowstart_indices()[block_index.second+1], first = (block_col == block_index.first ? local_sparsity.get_rowstart_indices()[block_index.second]+1 : local_sparsity.get_rowstart_indices()[block_index.second]); - + for (unsigned int j=first; j &boundar number new_rhs; if (matrix.block(block_index.first, block_index.first) .diag_element(block_index.second) != 0.0) - new_rhs = dof->second * + new_rhs = dof->second * matrix.block(block_index.first, block_index.first) .diag_element(block_index.second); else @@ -2204,10 +2260,10 @@ MatrixTools::apply_boundary_values (const std::map &boundar // store the only nonzero entry // of this line for the Gauss // elimination step - const number diagonal_entry + const number diagonal_entry = matrix.block(block_index.first,block_index.first) .diag_element(block_index.second); - + // we have to loop over all // rows of the matrix which // have a nonzero entry in @@ -2239,7 +2295,7 @@ MatrixTools::apply_boundary_values (const std::map &boundar = sparsity_pattern.block (block_row, block_index.first); const SparsityPattern &transpose_sparsity = sparsity_pattern.block (block_index.first, block_row); - + // traverse the row of // the transpose block // to find the @@ -2253,7 +2309,7 @@ MatrixTools::apply_boundary_values (const std::map &boundar transpose_sparsity.get_rowstart_indices()[block_index.second]+1 : transpose_sparsity.get_rowstart_indices()[block_index.second]), last = transpose_sparsity.get_rowstart_indices()[block_index.second+1]; - + for (unsigned int j=first; j &boundar (p != &this_sparsity.get_column_numbers() [this_sparsity.get_rowstart_indices()[row+1]]), ExcInternalError()); - + const unsigned int global_entry = (p - @@ -2336,7 +2392,7 @@ MatrixTools::apply_boundary_values (const std::map &boundar right_hand_side.block(block_row)(row) -= matrix.block(block_row,block_index.first).global_entry(global_entry) / diagonal_entry * new_rhs; - + // set matrix entry to zero matrix.block(block_row,block_index.first).global_entry(global_entry) = 0.; } @@ -2380,7 +2436,7 @@ namespace PETScWrappers ExcInternalError()); Assert (local_range == solution.local_range(), ExcInternalError()); - + // we have to read and write from this // matrix (in this order). this will only @@ -2501,7 +2557,7 @@ apply_boundary_values (const std::map &boundary_values, right_hand_side, eliminate_columns); // compress the matrix once we're done - matrix.compress (); + matrix.compress (); } @@ -2719,7 +2775,7 @@ namespace TrilinosWrappers ExcDimensionMismatch(matrix.n(), solution.size())); Assert (matrix.n_block_rows() == matrix.n_block_cols(), ExcNotQuadratic()); - + const unsigned int n_blocks = matrix.n_block_rows(); matrix.compress(); @@ -2746,7 +2802,7 @@ namespace TrilinosWrappers std::pair (index,dof->second)); } } - + // Now call the non-block variants on // the diagonal subblocks and the // solution/rhs. @@ -2756,8 +2812,8 @@ namespace TrilinosWrappers solution.block(block), right_hand_side.block(block), eliminate_columns); - - // Finally, we need to do something + + // Finally, we need to do something // about the off-diagonal matrices. This // is luckily not difficult. Just clear // the whole row. @@ -2765,7 +2821,7 @@ namespace TrilinosWrappers { const std::pair local_range = matrix.block(block_m,0).local_range(); - + std::vector constrained_rows; for (std::map::const_iterator dof = block_boundary_values[block_m].begin(); @@ -2774,7 +2830,7 @@ namespace TrilinosWrappers if ((dof->first >= local_range.first) && (dof->first < local_range.second)) constrained_rows.push_back (dof->first); - + for (unsigned int block_n=0; block_n &boundary_values, TrilinosWrappers::BlockVector &right_hand_side, const bool eliminate_columns) { - TrilinosWrappers::apply_block_boundary_values (boundary_values, matrix, - solution, right_hand_side, + TrilinosWrappers::apply_block_boundary_values (boundary_values, matrix, + solution, right_hand_side, eliminate_columns); } @@ -2840,8 +2896,8 @@ apply_boundary_values (const std::map &boundary_values, TrilinosWrappers::MPI::BlockVector &right_hand_side, const bool eliminate_columns) { - TrilinosWrappers::apply_block_boundary_values (boundary_values, matrix, - solution, right_hand_side, + TrilinosWrappers::apply_block_boundary_values (boundary_values, matrix, + solution, right_hand_side, eliminate_columns); } @@ -2943,7 +2999,7 @@ local_apply_boundary_values (const std::map &boundary_value } else local_matrix(i,i) = std::fabs(local_matrix(i,i)); - + // and replace rhs entry by correct // value local_rhs(i) = local_matrix(i,i) * boundary_value->second;