From: Denis Davydov Date: Tue, 23 Feb 2016 15:31:12 +0000 (+0100) Subject: continue templating VectorTools; introduce templates to MatrixTools and DoFTools X-Git-Tag: v8.5.0-rc1~1279^2~3 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=fb64b4cd6523a4984fbd9a3ad50b1096dc4a0837;p=dealii.git continue templating VectorTools; introduce templates to MatrixTools and DoFTools (i) drop the typedef FunctionMap from (hp::)DoFHandler and template n_boundary_dofs(). Fix a few places where it was used. (ii) template make_boundary_sparsity_pattern(...,FunctionMap,...) (iii) new templates in MatrixCreator:: create_mass_matrix(), create_boundary_mass_matrix(), apply_boundary_values() and local_apply_boundary_values(). (iv) number templates in VectorTools:: interpolate_boundary_values(), project_boundary_values(), (v) instantiate ConstraintMatrix::distribute_local_to_global for SparseMatrix --- diff --git a/include/deal.II/dofs/dof_handler.h b/include/deal.II/dofs/dof_handler.h index 727372449b..1e98e5cb14 100644 --- a/include/deal.II/dofs/dof_handler.h +++ b/include/deal.II/dofs/dof_handler.h @@ -258,11 +258,6 @@ public: typedef typename LevelSelector::face_iterator level_face_iterator; - /** - * Alias the @p FunctionMap type declared elsewhere. - */ - typedef typename dealii::FunctionMap::type FunctionMap; - /** * Make the dimension available in function templates. */ @@ -706,9 +701,12 @@ public: * reason that a @p map rather than a @p set is used is the same as * described in the section on the @p make_boundary_sparsity_pattern * function. + * + * The type of boundary_ids equals typename FunctionMap::type . */ + template types::global_dof_index - n_boundary_dofs (const FunctionMap &boundary_ids) const; + n_boundary_dofs (const std::map*> &boundary_ids) const; /** * Same function, but with different data type of the argument, which is @@ -1098,7 +1096,6 @@ private: #ifndef DOXYGEN template <> types::global_dof_index DoFHandler<1>::n_boundary_dofs () const; -template <> types::global_dof_index DoFHandler<1>::n_boundary_dofs (const FunctionMap &) const; template <> types::global_dof_index DoFHandler<1>::n_boundary_dofs (const std::set &) const; template <> void DoFHandler<1>::renumber_dofs(unsigned int,const std::vector &new_numbers); diff --git a/include/deal.II/dofs/dof_tools.h b/include/deal.II/dofs/dof_tools.h index b1eef9f9e1..4a9253d003 100644 --- a/include/deal.II/dofs/dof_tools.h +++ b/include/deal.II/dofs/dof_tools.h @@ -681,17 +681,17 @@ namespace DoFTools * This function could have been written by passing a @p set of boundary_id * numbers. However, most of the functions throughout deal.II dealing with * boundary indicators take a mapping of boundary indicators and the - * corresponding boundary function, i.e., a FunctionMap argument. + * corresponding boundary function, i.e., a std::map*> argument. * Correspondingly, this function does the same, though the actual boundary * function is ignored here. (Consequently, if you don't have any such * boundary functions, just create a map with the boundary indicators you * want and set the function pointers to null pointers). */ - template + template void make_boundary_sparsity_pattern (const DoFHandlerType &dof, - const typename FunctionMap::type &boundary_ids, + const std::map*> &boundary_ids, const std::vector &dof_to_boundary_mapping, SparsityPatternType &sparsity); diff --git a/include/deal.II/hp/dof_handler.h b/include/deal.II/hp/dof_handler.h index dd2e24e14f..0577b0f170 100644 --- a/include/deal.II/hp/dof_handler.h +++ b/include/deal.II/hp/dof_handler.h @@ -202,11 +202,6 @@ namespace hp typedef typename LevelSelector::face_iterator level_face_iterator; - /** - * Alias the @p FunctionMap type declared elsewhere. - */ - typedef typename FunctionMap::type FunctionMap; - /** * Make the dimension available in function templates. */ @@ -526,9 +521,12 @@ namespace hp * reason that a @p map rather than a @p set is used is the same as * described in the section on the @p make_boundary_sparsity_pattern * function. + * + * The type of @p boundary_ids equals typename FunctionMap::type. */ + template types::global_dof_index - n_boundary_dofs (const FunctionMap &boundary_ids) const; + n_boundary_dofs (const std::map*> &boundary_ids) const; /** * Same function, but with different data type of the argument, which is diff --git a/include/deal.II/numerics/matrix_tools.h b/include/deal.II/numerics/matrix_tools.h index 4784417ad8..9c6b0ea26d 100644 --- a/include/deal.II/numerics/matrix_tools.h +++ b/include/deal.II/numerics/matrix_tools.h @@ -53,6 +53,8 @@ namespace hp #ifdef DEAL_II_WITH_PETSC +# include + namespace PETScWrappers { class SparseMatrix; @@ -68,6 +70,9 @@ namespace PETScWrappers #endif #ifdef DEAL_II_WITH_TRILINOS +// TODO: move this to a stand-alone header and include from there? +# define TrilinosScalar double + namespace TrilinosWrappers { class SparseMatrix; @@ -175,8 +180,8 @@ namespace TrilinosWrappers * * The create_boundary_mass_matrix() creates the matrix with entries $m_{ij} = * \int_{\Gamma} \phi_i \phi_j dx$, where $\Gamma$ is the union of boundary - * parts with indicators contained in a FunctionMap passed to the function - * (i.e. if you want to set up the mass matrix for the parts of the boundary + * parts with indicators contained in a std::map*> + * passed to the function (i.e. if you want to set up the mass matrix for the parts of the boundary * with indicators zero and 2, you pass the function a map of unsigned * chars as parameter @p boundary_functions containing the keys zero and * 2). The size of the matrix is equal to the number of degrees of freedom @@ -232,23 +237,23 @@ namespace MatrixCreator * * See the general documentation of this namespace for more information. */ - template + template void create_mass_matrix (const Mapping &mapping, const DoFHandler &dof, const Quadrature &q, SparseMatrix &matrix, - const Function *const a = 0, + const Function *const a = 0, const ConstraintMatrix &constraints = ConstraintMatrix()); /** * Calls the create_mass_matrix() function, see above, with * mapping=MappingQGeneric@(1). */ - template + template void create_mass_matrix (const DoFHandler &dof, const Quadrature &q, SparseMatrix &matrix, - const Function *const a = 0, + const Function *const a = 0, const ConstraintMatrix &constraints = ConstraintMatrix()); /** @@ -267,73 +272,73 @@ namespace MatrixCreator * * See the general documentation of this namespace for more information. */ - template + template void create_mass_matrix (const Mapping &mapping, const DoFHandler &dof, const Quadrature &q, SparseMatrix &matrix, - const Function &rhs, + const Function &rhs, Vector &rhs_vector, - const Function *const a = 0, + const Function *const a = 0, const ConstraintMatrix &constraints = ConstraintMatrix()); /** * Calls the create_mass_matrix() function, see above, with * mapping=MappingQGeneric@(1). */ - template + template void create_mass_matrix (const DoFHandler &dof, const Quadrature &q, SparseMatrix &matrix, - const Function &rhs, + const Function &rhs, Vector &rhs_vector, - const Function *const a = 0, + const Function *const a = 0, const ConstraintMatrix &constraints = ConstraintMatrix()); /** * Same function as above, but for hp objects. */ - template + template void create_mass_matrix (const hp::MappingCollection &mapping, const hp::DoFHandler &dof, const hp::QCollection &q, SparseMatrix &matrix, - const Function *const a = 0, + const Function *const a = 0, const ConstraintMatrix &constraints = ConstraintMatrix()); /** * Same function as above, but for hp objects. */ - template + template void create_mass_matrix (const hp::DoFHandler &dof, const hp::QCollection &q, SparseMatrix &matrix, - const Function *const a = 0, + const Function *const a = 0, const ConstraintMatrix &constraints = ConstraintMatrix()); /** * Same function as above, but for hp objects. */ - template + template void create_mass_matrix (const hp::MappingCollection &mapping, const hp::DoFHandler &dof, const hp::QCollection &q, SparseMatrix &matrix, - const Function &rhs, + const Function &rhs, Vector &rhs_vector, - const Function *const a = 0, + const Function *const a = 0, const ConstraintMatrix &constraints = ConstraintMatrix()); /** * Same function as above, but for hp objects. */ - template + template void create_mass_matrix (const hp::DoFHandler &dof, const hp::QCollection &q, SparseMatrix &matrix, - const Function &rhs, + const Function &rhs, Vector &rhs_vector, - const Function *const a = 0, + const Function *const a = 0, const ConstraintMatrix &constraints = ConstraintMatrix()); @@ -359,15 +364,15 @@ namespace MatrixCreator * @todo This function does not work for finite elements with cell-dependent * shape functions. */ - template + 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, + SparseMatrix &matrix, + const std::map*> &boundary_functions, + Vector &rhs_vector, std::vector &dof_to_boundary_mapping, - const Function *const weight = 0, + const Function *const weight = 0, std::vector component_mapping = std::vector()); @@ -375,41 +380,41 @@ namespace MatrixCreator * Calls the create_boundary_mass_matrix() function, see above, with * mapping=MappingQGeneric@(1). */ - template + template void create_boundary_mass_matrix (const DoFHandler &dof, const Quadrature &q, - SparseMatrix &matrix, - const typename FunctionMap::type &boundary_functions, - Vector &rhs_vector, + SparseMatrix &matrix, + const std::map*> &boundary_functions, + Vector &rhs_vector, std::vector &dof_to_boundary_mapping, - const Function *const a = 0, + const Function *const a = 0, std::vector component_mapping = std::vector()); /** * Same function as above, but for hp objects. */ - template + 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, + SparseMatrix &matrix, + const std::map*> &boundary_functions, + Vector &rhs_vector, std::vector &dof_to_boundary_mapping, - const Function *const a = 0, + const Function *const a = 0, std::vector component_mapping = std::vector()); /** * Same function as above, but for hp objects. */ - template + template void create_boundary_mass_matrix (const hp::DoFHandler &dof, const hp::QCollection &q, - SparseMatrix &matrix, - const typename FunctionMap::type &boundary_functions, - Vector &rhs_vector, + SparseMatrix &matrix, + const std::map*> &boundary_functions, + Vector &rhs_vector, std::vector &dof_to_boundary_mapping, - const Function *const a = 0, + const Function *const a = 0, std::vector component_mapping = std::vector()); /** @@ -757,7 +762,7 @@ namespace MatrixTools */ template void - apply_boundary_values (const std::map &boundary_values, + apply_boundary_values (const std::map &boundary_values, SparseMatrix &matrix, Vector &solution, Vector &right_hand_side, @@ -770,7 +775,7 @@ namespace MatrixTools */ template void - apply_boundary_values (const std::map &boundary_values, + apply_boundary_values (const std::map &boundary_values, BlockSparseMatrix &matrix, BlockVector &solution, BlockVector &right_hand_side, @@ -797,7 +802,7 @@ namespace MatrixTools * This function is used in step-17 and step-18. */ void - apply_boundary_values (const std::map &boundary_values, + apply_boundary_values (const std::map &boundary_values, PETScWrappers::SparseMatrix &matrix, PETScWrappers::Vector &solution, PETScWrappers::Vector &right_hand_side, @@ -823,7 +828,7 @@ namespace MatrixTools * rows. */ void - apply_boundary_values (const std::map &boundary_values, + apply_boundary_values (const std::map &boundary_values, PETScWrappers::MPI::SparseMatrix &matrix, PETScWrappers::MPI::Vector &solution, PETScWrappers::MPI::Vector &right_hand_side, @@ -833,7 +838,7 @@ namespace MatrixTools * Same function, but for parallel PETSc matrices and a non-parallel vector. */ void - apply_boundary_values (const std::map &boundary_values, + apply_boundary_values (const std::map &boundary_values, PETScWrappers::MPI::SparseMatrix &matrix, PETScWrappers::Vector &solution, PETScWrappers::MPI::Vector &right_hand_side, @@ -843,7 +848,7 @@ namespace MatrixTools * Same as above but for BlockSparseMatrix. */ void - apply_boundary_values (const std::map &boundary_values, + apply_boundary_values (const std::map &boundary_values, PETScWrappers::MPI::BlockSparseMatrix &matrix, PETScWrappers::MPI::BlockVector &solution, PETScWrappers::MPI::BlockVector &right_hand_side, @@ -870,7 +875,7 @@ namespace MatrixTools * name in this namespace.) */ void - apply_boundary_values (const std::map &boundary_values, + apply_boundary_values (const std::map &boundary_values, TrilinosWrappers::SparseMatrix &matrix, TrilinosWrappers::Vector &solution, TrilinosWrappers::Vector &right_hand_side, @@ -881,7 +886,7 @@ namespace MatrixTools * structures. */ void - apply_boundary_values (const std::map &boundary_values, + apply_boundary_values (const std::map &boundary_values, TrilinosWrappers::BlockSparseMatrix &matrix, TrilinosWrappers::BlockVector &solution, TrilinosWrappers::BlockVector &right_hand_side, @@ -907,7 +912,7 @@ namespace MatrixTools * rows. */ void - apply_boundary_values (const std::map &boundary_values, + apply_boundary_values (const std::map &boundary_values, TrilinosWrappers::SparseMatrix &matrix, TrilinosWrappers::MPI::Vector &solution, TrilinosWrappers::MPI::Vector &right_hand_side, @@ -918,7 +923,7 @@ namespace MatrixTools * structures. */ void - apply_boundary_values (const std::map &boundary_values, + apply_boundary_values (const std::map &boundary_values, TrilinosWrappers::BlockSparseMatrix &matrix, TrilinosWrappers::MPI::BlockVector &solution, TrilinosWrappers::MPI::BlockVector &right_hand_side, @@ -943,11 +948,12 @@ namespace MatrixTools * * @dealiiVideoLecture{21.6,21.65} */ + template void - local_apply_boundary_values (const std::map &boundary_values, + local_apply_boundary_values (const std::map &boundary_values, const std::vector &local_dof_indices, - FullMatrix &local_matrix, - Vector &local_rhs, + FullMatrix &local_matrix, + Vector &local_rhs, const bool eliminate_columns); /** diff --git a/include/deal.II/numerics/vector_tools.h b/include/deal.II/numerics/vector_tools.h index 64f046f0c4..3203daec37 100644 --- a/include/deal.II/numerics/vector_tools.h +++ b/include/deal.II/numerics/vector_tools.h @@ -32,7 +32,6 @@ DEAL_II_NAMESPACE_OPEN template class Function; -template struct FunctionMap; template class Quadrature; template class QGauss; @@ -111,7 +110,8 @@ class ConstraintMatrix; * MatrixTools::apply_boundary_values() function. The projection of the trace * of the function to the boundary is done with the * VectorTools::project_boundary_values() (see below) function, which is - * called with a map of boundary functions FunctionMap in which all boundary + * called with a map of boundary functions + * std::map*> in which all boundary * indicators from zero to numbers::internal_face_boundary_id-1 * (numbers::internal_face_boundary_id is used for other purposes, see the * Triangulation class documentation) point to the function to be projected. @@ -217,7 +217,8 @@ class ConstraintMatrix; * trace of the function to the boundary. * * The projection takes place on all boundary parts with boundary indicators - * listed in the map (FunctioMap::FunctionMap) of boundary functions. These + * listed in the map (std::map*>) + * of boundary functions. These * boundary parts may or may not be continuous. For these boundary parts, the * mass matrix is assembled using the * MatrixTools::create_boundary_mass_matrix() function, as well as the @@ -698,26 +699,26 @@ namespace VectorTools * * See the general documentation of this namespace for more information. */ - template + template void interpolate_boundary_values (const Mapping &mapping, const DoFHandlerType &dof, - const typename FunctionMap::type &function_map, - std::map &boundary_values, + const std::map*> &function_map, + std::map &boundary_values, const ComponentMask &component_mask = ComponentMask()); /** * Like the previous function, but take a mapping collection to go with the * hp::DoFHandler object. */ - template + template void interpolate_boundary_values (const hp::MappingCollection &mapping, const hp::DoFHandler &dof, - const typename FunctionMap::type &function_map, - std::map &boundary_values, + const std::map*> &function_map, + std::map &boundary_values, const ComponentMask &component_mask = ComponentMask()); /** @@ -729,14 +730,14 @@ namespace VectorTools * @see * @ref GlossBoundaryIndicator "Glossary entry on boundary indicators" */ - template + template void interpolate_boundary_values (const Mapping &mapping, const DoFHandlerType &dof, const types::boundary_id boundary_component, - const Function &boundary_function, - std::map &boundary_values, + const Function &boundary_function, + std::map &boundary_values, const ComponentMask &component_mask = ComponentMask()); /** @@ -748,13 +749,13 @@ namespace VectorTools * @see * @ref GlossBoundaryIndicator "Glossary entry on boundary indicators" */ - template + template void interpolate_boundary_values (const DoFHandlerType &dof, const types::boundary_id boundary_component, - const Function &boundary_function, - std::map &boundary_values, + const Function &boundary_function, + std::map &boundary_values, const ComponentMask &component_mask = ComponentMask()); @@ -764,12 +765,12 @@ namespace VectorTools * apply as for the previous function, in particular about the use of the * component mask and the requires size of the function object. */ - template + template void interpolate_boundary_values (const DoFHandlerType &dof, - const typename FunctionMap::type &function_map, - std::map &boundary_values, + const std::map*> &function_map, + std::map &boundary_values, const ComponentMask &component_mask = ComponentMask()); @@ -834,14 +835,14 @@ namespace VectorTools * * @ingroup constraints */ - template + template void interpolate_boundary_values - (const Mapping &mapping, - const DoFHandlerType &dof, - const typename FunctionMap::type &function_map, - ConstraintMatrix &constraints, - const ComponentMask &component_mask = ComponentMask()); + (const Mapping &mapping, + const DoFHandlerType &dof, + const std::map*> &function_map, + ConstraintMatrix &constraints, + const ComponentMask &component_mask = ComponentMask()); /** * Same function as above, but taking only one pair of boundary indicator @@ -854,13 +855,13 @@ namespace VectorTools * @see * @ref GlossBoundaryIndicator "Glossary entry on boundary indicators" */ - template + template void interpolate_boundary_values (const Mapping &mapping, const DoFHandlerType &dof, const types::boundary_id boundary_component, - const Function &boundary_function, + const Function &boundary_function, ConstraintMatrix &constraints, const ComponentMask &component_mask = ComponentMask()); @@ -875,12 +876,12 @@ namespace VectorTools * @see * @ref GlossBoundaryIndicator "Glossary entry on boundary indicators" */ - template + template void interpolate_boundary_values (const DoFHandlerType &dof, const types::boundary_id boundary_component, - const Function &boundary_function, + const Function &boundary_function, ConstraintMatrix &constraints, const ComponentMask &component_mask = ComponentMask()); @@ -893,11 +894,11 @@ namespace VectorTools * * @ingroup constraints */ - template + template void interpolate_boundary_values (const DoFHandlerType &dof, - const typename FunctionMap::type &function_map, + const std::map*> &function_map, ConstraintMatrix &constraints, const ComponentMask &component_mask = ComponentMask()); @@ -952,45 +953,45 @@ namespace VectorTools * component number in @p boundary_functions that should be used for this * component in @p dof. By default, no remapping is applied. */ - template + template void project_boundary_values (const Mapping &mapping, const DoFHandler &dof, - const typename FunctionMap::type &boundary_functions, + const std::map*> &boundary_functions, const Quadrature &q, - std::map &boundary_values, + std::map &boundary_values, std::vector component_mapping = std::vector()); /** * Calls the project_boundary_values() function, see above, with * mapping=MappingQGeneric@(1). */ - template + template void project_boundary_values (const DoFHandler &dof, - const typename FunctionMap::type &boundary_function, + const std::map*> &boundary_function, const Quadrature &q, - std::map &boundary_values, + std::map &boundary_values, std::vector component_mapping = std::vector()); /** * Same as above, but for objects of type hp::DoFHandler */ - template + template void project_boundary_values (const hp::MappingCollection &mapping, const hp::DoFHandler &dof, - const typename FunctionMap::type &boundary_functions, + const std::map*> &boundary_functions, const hp::QCollection &q, - std::map &boundary_values, + std::map &boundary_values, std::vector component_mapping = std::vector()); /** * Calls the project_boundary_values() function, see above, with * mapping=MappingQGeneric@(1). */ - template + template void project_boundary_values (const hp::DoFHandler &dof, - const typename FunctionMap::type &boundary_function, + const std::map*> &boundary_function, const hp::QCollection &q, - std::map &boundary_values, + std::map &boundary_values, std::vector component_mapping = std::vector()); /** @@ -1031,10 +1032,10 @@ namespace VectorTools * * @ingroup constraints */ - template + template void project_boundary_values (const Mapping &mapping, const DoFHandler &dof, - const typename FunctionMap::type &boundary_functions, + const std::map*> &boundary_functions, const Quadrature &q, ConstraintMatrix &constraints, std::vector component_mapping = std::vector()); @@ -1045,9 +1046,9 @@ namespace VectorTools * * @ingroup constraints */ - template + template void project_boundary_values (const DoFHandler &dof, - const typename FunctionMap::type &boundary_function, + const std::map*> &boundary_function, const Quadrature &q, ConstraintMatrix &constraints, std::vector component_mapping = std::vector()); diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h index f2a40b952b..68d9cc1983 100644 --- a/include/deal.II/numerics/vector_tools.templates.h +++ b/include/deal.II/numerics/vector_tools.templates.h @@ -500,10 +500,10 @@ namespace VectorTools * mapping here because the function we evaluate for the DoFs is zero in * the mapped locations as well as in the original, unmapped locations */ - template + template void interpolate_zero_boundary_values (const DoFHandlerType &dof_handler, - std::map &boundary_values) + std::map &boundary_values) { const unsigned int dim = DoFHandlerType::dimension; @@ -690,15 +690,16 @@ namespace VectorTools * Compute the boundary values to be used in the project() functions. */ template class DoFHandlerType, - template class M_or_MC, template class Q_or_QC> + template class M_or_MC, template class Q_or_QC, + typename number> void project_compute_b_v (const M_or_MC &mapping, const DoFHandlerType &dof, - const Function &function, + const Function &function, const bool enforce_zero_boundary, const Q_or_QC &q_boundary, const bool project_to_boundary_first, - std::map &boundary_values) + std::map &boundary_values) { if (enforce_zero_boundary == true) // no need to project boundary @@ -722,7 +723,7 @@ namespace VectorTools const std::vector used_boundary_ids = dof.get_triangulation().get_boundary_ids(); - typename FunctionMap::type boundary_functions; + std::map*> boundary_functions; for (unsigned int i=0; i bool constraints_and_b_v_are_compatible (const ConstraintMatrix &constraints, - std::map &boundary_values) + std::map &boundary_values) { - for (std::map::iterator it=boundary_values.begin(); + for (typename std::map::iterator it=boundary_values.begin(); it != boundary_values.end(); ++it) if (constraints.is_constrained(it->first)) //TODO: This looks wrong -- shouldn't it be ==0 in the first condition and && ? @@ -762,12 +763,13 @@ namespace VectorTools const DoFHandlerType &dof, const ConstraintMatrix &constraints, const Q_or_QC &quadrature, - const Function &function, + const Function &function, VectorType &vec_result, const bool enforce_zero_boundary, const Q_or_QC &q_boundary, const bool project_to_boundary_first) { + typedef typename VectorType::value_type number; Assert (dof.get_fe().n_components() == function.n_components, ExcDimensionMismatch(dof.get_fe().n_components(), function.n_components)); @@ -775,16 +777,16 @@ namespace VectorTools ExcDimensionMismatch (vec_result.size(), dof.n_dofs())); // make up boundary values - std::map boundary_values; + std::map boundary_values; project_compute_b_v(mapping, dof, function, enforce_zero_boundary, q_boundary, project_to_boundary_first, boundary_values); // check if constraints are compatible (see below) const bool constraints_are_compatible = - constraints_and_b_v_are_compatible(constraints, boundary_values); + constraints_and_b_v_are_compatible(constraints, boundary_values); // set up mass matrix and right hand side - Vector vec (dof.n_dofs()); + Vector vec (dof.n_dofs()); SparsityPattern sparsity; { DynamicSparsityPattern dsp (dof.n_dofs(), dof.n_dofs()); @@ -793,8 +795,8 @@ namespace VectorTools sparsity.copy_from (dsp); } - SparseMatrix mass_matrix (sparsity); - Vector tmp (mass_matrix.n()); + SparseMatrix mass_matrix (sparsity); + Vector tmp (mass_matrix.n()); // If the constraint matrix does not conflict with the given boundary // values (i.e., it either does not contain boundary values or it contains @@ -803,7 +805,7 @@ namespace VectorTools // interpolate the boundary values and then condense the matrix and vector if (constraints_are_compatible) { - const Function *dummy = 0; + const Function *dummy = 0; MatrixCreator::create_mass_matrix (mapping, dof, quadrature, mass_matrix, function, tmp, dummy, constraints); @@ -827,10 +829,10 @@ namespace VectorTools // steps may not be sufficient, since roundoff errors may accumulate for // badly conditioned matrices ReductionControl control(5*tmp.size(), 0., 1e-12, false, false); - GrowingVectorMemory<> memory; - SolverCG<> cg(control,memory); + GrowingVectorMemory > memory; + SolverCG > cg(control,memory); - PreconditionSSOR<> prec; + PreconditionSSOR > prec; prec.initialize(mass_matrix, 1.2); cg.solve (mass_matrix, vec, tmp, prec); @@ -851,7 +853,7 @@ namespace VectorTools const DoFHandler &dof, const ConstraintMatrix &constraints, const Quadrature &quadrature, - const Function &function, + const Function &function, VectorType &vec_result, const bool enforce_zero_boundary, const Quadrature &q_boundary, @@ -868,7 +870,7 @@ namespace VectorTools void project (const DoFHandler &dof, const ConstraintMatrix &constraints, const Quadrature &quadrature, - const Function &function, + const Function &function, VectorType &vec, const bool enforce_zero_boundary, const Quadrature &q_boundary, @@ -885,7 +887,7 @@ namespace VectorTools const hp::DoFHandler &dof, const ConstraintMatrix &constraints, const hp::QCollection &quadrature, - const Function &function, + const Function &function, VectorType &vec_result, const bool enforce_zero_boundary, const hp::QCollection &q_boundary, @@ -902,7 +904,7 @@ namespace VectorTools void project (const hp::DoFHandler &dof, const ConstraintMatrix &constraints, const hp::QCollection &quadrature, - const Function &function, + const Function &function, VectorType &vec, const bool enforce_zero_boundary, const hp::QCollection &q_boundary, @@ -1665,14 +1667,14 @@ namespace VectorTools // faces are points and it is far // easier to simply work on // individual vertices - template class M_or_MC> + template class M_or_MC> static inline void do_interpolate_boundary_values (const M_or_MC &, - const DoFHandlerType &dof, - const typename FunctionMap::type &function_map, - std::map &boundary_values, - const ComponentMask &component_mask, + const DoFHandlerType &dof, + const std::map*> &function_map, + std::map &boundary_values, + const ComponentMask &component_mask, const dealii::internal::int2type<1>) { const unsigned int dim = DoFHandlerType::dimension; @@ -1697,7 +1699,7 @@ namespace VectorTools && (function_map.find(cell->face(direction)->boundary_id()) != function_map.end())) { - const Function &boundary_function + const Function &boundary_function = *function_map.find(cell->face(direction)->boundary_id())->second; // get the FE corresponding to this @@ -1725,7 +1727,7 @@ namespace VectorTools // of the number of // components of the // function - Vector function_values (fe.n_components()); + Vector function_values (fe.n_components()); if (fe.n_components() == 1) function_values(0) = boundary_function.value (cell->vertex(direction)); @@ -1748,15 +1750,15 @@ namespace VectorTools // dim_, it is clearly less specialized than the 1d function above and // whenever possible (i.e., if dim==1), the function template above // will be used - template class M_or_MC, int dim_> + template class M_or_MC, int dim_> static inline void do_interpolate_boundary_values - (const M_or_MC &mapping, - const DoFHandlerType &dof, - const typename FunctionMap::type &function_map, - std::map &boundary_values, - const ComponentMask &component_mask, + (const M_or_MC &mapping, + const DoFHandlerType &dof, + const std::map*> &function_map, + std::map &boundary_values, + const ComponentMask &component_mask, const dealii::internal::int2type) { const unsigned int dim = DoFHandlerType::dimension; @@ -1780,7 +1782,7 @@ namespace VectorTools const unsigned int n_components = DoFTools::n_components(dof); const bool fe_is_system = (n_components != 1); - for (typename FunctionMap::type::const_iterator i=function_map.begin(); + for (typename std::map*>::const_iterator i=function_map.begin(); i!=function_map.end(); ++i) Assert (n_components == i->second->n_components, ExcDimensionMismatch(n_components, i->second->n_components)); @@ -1795,8 +1797,8 @@ namespace VectorTools // array to store the values of the boundary function at the boundary // points. have two arrays for scalar and vector functions to use the // more efficient one respectively - std::vector dof_values_scalar; - std::vector > dof_values_system; + std::vector dof_values_scalar; + std::vector > dof_values_system; dof_values_scalar.reserve (DoFTools::max_dofs_per_face (dof)); dof_values_system.reserve (DoFTools::max_dofs_per_face (dof)); @@ -1909,7 +1911,7 @@ namespace VectorTools // allocating temporary if possible if (dof_values_system.size() < fe.dofs_per_face) dof_values_system.resize (fe.dofs_per_face, - Vector(fe.n_components())); + Vector(fe.n_components())); else dof_values_system.resize (fe.dofs_per_face); @@ -1991,15 +1993,14 @@ namespace VectorTools - template + template void - interpolate_boundary_values - (const Mapping &mapping, - const DoFHandlerType &dof, - const typename FunctionMap::type &function_map, - std::map &boundary_values, - const ComponentMask &component_mask_) + (const Mapping &mapping, + const DoFHandlerType &dof, + const std::map*> &function_map, + std::map &boundary_values, + const ComponentMask &component_mask_) { do_interpolate_boundary_values (mapping, dof, function_map, boundary_values, component_mask_, @@ -2008,31 +2009,31 @@ namespace VectorTools - template + template void interpolate_boundary_values (const Mapping &mapping, const DoFHandlerType &dof, const types::boundary_id boundary_component, - const Function &boundary_function, - std::map &boundary_values, + const Function &boundary_function, + std::map &boundary_values, const ComponentMask &component_mask) { - typename FunctionMap::type function_map; + std::map*> function_map; function_map[boundary_component] = &boundary_function; interpolate_boundary_values (mapping, dof, function_map, boundary_values, component_mask); } - template + template void interpolate_boundary_values - (const hp::MappingCollection &mapping, - const hp::DoFHandler &dof, - const typename FunctionMap::type &function_map, - std::map &boundary_values, - const ComponentMask &component_mask_) + (const hp::MappingCollection &mapping, + const hp::DoFHandler &dof, + const std::map*> &function_map, + std::map &boundary_values, + const ComponentMask &component_mask_) { do_interpolate_boundary_values (mapping, dof, function_map, boundary_values, component_mask_, @@ -2041,14 +2042,14 @@ namespace VectorTools - template + template void interpolate_boundary_values - (const DoFHandlerType &dof, - const types::boundary_id boundary_component, - const Function &boundary_function, - std::map &boundary_values, - const ComponentMask &component_mask) + (const DoFHandlerType &dof, + const types::boundary_id boundary_component, + const Function &boundary_function, + std::map &boundary_values, + const ComponentMask &component_mask) { interpolate_boundary_values(StaticMappingQ1::mapping, dof, boundary_component, @@ -2057,13 +2058,13 @@ namespace VectorTools - template + template void interpolate_boundary_values - (const DoFHandlerType &dof, - const typename FunctionMap::type &function_map, - std::map &boundary_values, - const ComponentMask &component_mask) + (const DoFHandlerType &dof, + const std::map*> &function_map, + std::map &boundary_values, + const ComponentMask &component_mask) { interpolate_boundary_values (StaticMappingQ1::mapping, @@ -2078,19 +2079,19 @@ namespace VectorTools - template + template void interpolate_boundary_values - (const Mapping &mapping, - const DoFHandlerType &dof, - const typename FunctionMap::type &function_map, - ConstraintMatrix &constraints, - const ComponentMask &component_mask_) + (const Mapping &mapping, + const DoFHandlerType &dof, + const std::map*> &function_map, + ConstraintMatrix &constraints, + const ComponentMask &component_mask_) { - std::map boundary_values; + std::map boundary_values; interpolate_boundary_values (mapping, dof, function_map, boundary_values, component_mask_); - std::map::const_iterator boundary_value = + typename std::map::const_iterator boundary_value = boundary_values.begin(); for ( ; boundary_value !=boundary_values.end(); ++boundary_value) { @@ -2107,17 +2108,17 @@ namespace VectorTools - template + template void interpolate_boundary_values (const Mapping &mapping, const DoFHandlerType &dof, const types::boundary_id boundary_component, - const Function &boundary_function, + const Function &boundary_function, ConstraintMatrix &constraints, const ComponentMask &component_mask) { - typename FunctionMap::type function_map; + std::map*> function_map; function_map[boundary_component] = &boundary_function; interpolate_boundary_values (mapping, dof, function_map, constraints, component_mask); @@ -2125,14 +2126,14 @@ namespace VectorTools - template + template void interpolate_boundary_values - (const DoFHandlerType &dof, - const types::boundary_id boundary_component, - const Function &boundary_function, - ConstraintMatrix &constraints, - const ComponentMask &component_mask) + (const DoFHandlerType &dof, + const types::boundary_id boundary_component, + const Function &boundary_function, + ConstraintMatrix &constraints, + const ComponentMask &component_mask) { interpolate_boundary_values (StaticMappingQ1::mapping, @@ -2142,13 +2143,13 @@ namespace VectorTools - template + template void interpolate_boundary_values - (const DoFHandlerType &dof, - const typename FunctionMap::type &function_map, - ConstraintMatrix &constraints, - const ComponentMask &component_mask) + (const DoFHandlerType &dof, + const std::map*> &function_map, + ConstraintMatrix &constraints, + const ComponentMask &component_mask) { interpolate_boundary_values (StaticMappingQ1::mapping, @@ -2165,14 +2166,14 @@ namespace VectorTools namespace { template class DoFHandlerType, - template class M_or_MC, template class Q_or_QC> + template class M_or_MC, template class Q_or_QC, typename number> void - do_project_boundary_values (const M_or_MC &mapping, - const DoFHandlerType &dof, - const typename FunctionMap::type &boundary_functions, - const Q_or_QC &q, - std::map &boundary_values, - std::vector component_mapping) + do_project_boundary_values (const M_or_MC &mapping, + const DoFHandlerType &dof, + const std::map*> &boundary_functions, + const Q_or_QC &q, + std::map &boundary_values, + std::vector component_mapping) { // in 1d, projection onto the 0d end points == interpolation if (dim == 1) @@ -2204,7 +2205,7 @@ namespace VectorTools std::vector dof_to_boundary_mapping; std::set selected_boundary_components; - for (typename FunctionMap::type::const_iterator i=boundary_functions.begin(); + for (typename std::map*>::const_iterator i=boundary_functions.begin(); i!=boundary_functions.end(); ++i) selected_boundary_components.insert (i->first); @@ -2216,8 +2217,8 @@ namespace VectorTools return; // set up sparsity structure - DynamicSparsityPattern dsp(dof.n_boundary_dofs (boundary_functions), - dof.n_boundary_dofs (boundary_functions)); + DynamicSparsityPattern dsp(dof.n_boundary_dofs(boundary_functions), + dof.n_boundary_dofs(boundary_functions)); DoFTools::make_boundary_sparsity_pattern (dof, boundary_functions, dof_to_boundary_mapping, @@ -2266,13 +2267,13 @@ namespace VectorTools // make mass matrix and right hand side - SparseMatrix mass_matrix(sparsity); - Vector rhs(sparsity.n_rows()); + SparseMatrix mass_matrix(sparsity); + Vector rhs(sparsity.n_rows()); MatrixCreator::create_boundary_mass_matrix (mapping, dof, q, mass_matrix, boundary_functions, - rhs, dof_to_boundary_mapping, (const Function *) 0, + rhs, dof_to_boundary_mapping, (const Function *) 0, component_mapping); // For certain weird elements, @@ -2287,11 +2288,11 @@ namespace VectorTools //TODO: Maybe we should figure out if the element really needs this - FilteredMatrix > filtered_mass_matrix(mass_matrix, true); - FilteredMatrix > filtered_precondition; + FilteredMatrix > filtered_mass_matrix(mass_matrix, true); + FilteredMatrix > filtered_precondition; std::vector excluded_dofs(mass_matrix.m(), false); - double max_element = 0.; + number max_element = 0.; for (unsigned int i=0; i max_element) max_element = mass_matrix.diag_element(i); @@ -2305,11 +2306,11 @@ namespace VectorTools excluded_dofs[i] = true; } - Vector boundary_projection (rhs.size()); + Vector boundary_projection (rhs.size()); // cannot reduce residual in a useful way if we are close to the square // root of the minimal double value - if (rhs.norm_sqr() < 1e28 * std::numeric_limits::min()) + if (rhs.norm_sqr() < 1e28 * std::numeric_limits::min()) boundary_projection = 0; else { @@ -2317,10 +2318,10 @@ namespace VectorTools // steps may not be sufficient, since roundoff errors may accumulate for // badly conditioned matrices ReductionControl control(5*rhs.size(), 0., 1.e-12, false, false); - GrowingVectorMemory<> memory; - SolverCG<> cg(control,memory); + GrowingVectorMemory > memory; + SolverCG > cg(control,memory); - PreconditionSSOR<> prec; + PreconditionSSOR > prec; prec.initialize(mass_matrix, 1.2); filtered_precondition.initialize(prec, true); // solve @@ -2347,14 +2348,14 @@ namespace VectorTools } } - template + template void - project_boundary_values (const Mapping &mapping, - const DoFHandler &dof, - const typename FunctionMap::type &boundary_functions, - const Quadrature &q, - std::map &boundary_values, - std::vector component_mapping) + project_boundary_values (const Mapping &mapping, + const DoFHandler &dof, + const std::map*> &boundary_functions, + const Quadrature &q, + std::map &boundary_values, + std::vector component_mapping) { do_project_boundary_values(mapping, dof, boundary_functions, q, boundary_values, component_mapping); @@ -2362,13 +2363,13 @@ namespace VectorTools - template + template void - project_boundary_values (const DoFHandler &dof, - const typename FunctionMap::type &boundary_functions, - const Quadrature &q, - std::map &boundary_values, - std::vector component_mapping) + project_boundary_values (const DoFHandler &dof, + const std::map*> &boundary_functions, + const Quadrature &q, + std::map &boundary_values, + std::vector component_mapping) { project_boundary_values(StaticMappingQ1::mapping, dof, boundary_functions, q, boundary_values, component_mapping); @@ -2376,13 +2377,13 @@ namespace VectorTools - template - void project_boundary_values (const hp::MappingCollection &mapping, - const hp::DoFHandler &dof, - const typename FunctionMap::type &boundary_functions, - const hp::QCollection &q, - std::map &boundary_values, - std::vector component_mapping) + template + void project_boundary_values (const hp::MappingCollection &mapping, + const hp::DoFHandler &dof, + const std::map*> &boundary_functions, + const hp::QCollection &q, + std::map &boundary_values, + std::vector component_mapping) { do_project_boundary_values (mapping, dof, boundary_functions, @@ -2392,12 +2393,12 @@ namespace VectorTools - template - void project_boundary_values (const hp::DoFHandler &dof, - const typename FunctionMap::type &boundary_function, - const hp::QCollection &q, - std::map &boundary_values, - std::vector component_mapping) + template + void project_boundary_values (const hp::DoFHandler &dof, + const std::map*> &boundary_function, + const hp::QCollection &q, + std::map &boundary_values, + std::vector component_mapping) { project_boundary_values (hp::StaticMappingQ1::mapping_collection, dof, boundary_function, @@ -2410,19 +2411,19 @@ namespace VectorTools - template + template void - project_boundary_values (const Mapping &mapping, - const DoFHandler &dof, - const typename FunctionMap::type &boundary_functions, - const Quadrature &q, - ConstraintMatrix &constraints, - std::vector component_mapping) + project_boundary_values (const Mapping &mapping, + const DoFHandler &dof, + const std::map*> &boundary_functions, + const Quadrature &q, + ConstraintMatrix &constraints, + std::vector component_mapping) { - std::map boundary_values; + std::map boundary_values; project_boundary_values (mapping, dof, boundary_functions, q, boundary_values, component_mapping); - std::map::const_iterator boundary_value = + typename std::map::const_iterator boundary_value = boundary_values.begin(); for ( ; boundary_value !=boundary_values.end(); ++boundary_value) { @@ -2437,13 +2438,13 @@ namespace VectorTools - template + template void - project_boundary_values (const DoFHandler &dof, - const typename FunctionMap::type &boundary_functions, - const Quadrature &q, - ConstraintMatrix &constraints, - std::vector component_mapping) + project_boundary_values (const DoFHandler &dof, + const std::map*> &boundary_functions, + const Quadrature &q, + ConstraintMatrix &constraints, + std::vector component_mapping) { project_boundary_values(StaticMappingQ1::mapping, dof, boundary_functions, q, constraints, component_mapping); diff --git a/source/dofs/dof_handler.cc b/source/dofs/dof_handler.cc index 4627df3861..9b992decb3 100644 --- a/source/dofs/dof_handler.cc +++ b/source/dofs/dof_handler.cc @@ -1004,12 +1004,13 @@ types::global_dof_index DoFHandler<1>::n_boundary_dofs () const template <> -types::global_dof_index DoFHandler<1>::n_boundary_dofs (const FunctionMap &boundary_ids) const +template +types::global_dof_index DoFHandler<1>::n_boundary_dofs (const std::map*> &boundary_ids) const { // check that only boundary // indicators 0 and 1 are allowed // in 1d - for (FunctionMap::const_iterator i=boundary_ids.begin(); + for (typename std::map*>::const_iterator i=boundary_ids.begin(); i!=boundary_ids.end(); ++i) Assert ((i->first == 0) || (i->first == 1), ExcInvalidBoundaryIndicator()); @@ -1043,12 +1044,13 @@ types::global_dof_index DoFHandler<1,2>::n_boundary_dofs () const template <> -types::global_dof_index DoFHandler<1,2>::n_boundary_dofs (const FunctionMap &boundary_ids) const +template +types::global_dof_index DoFHandler<1,2>::n_boundary_dofs (const std::map*> &boundary_ids) const { // check that only boundary // indicators 0 and 1 are allowed // in 1d - for (FunctionMap::const_iterator i=boundary_ids.begin(); + for (typename std::map*>::const_iterator i=boundary_ids.begin(); i!=boundary_ids.end(); ++i) Assert ((i->first == 0) || (i->first == 1), ExcInvalidBoundaryIndicator()); @@ -1113,8 +1115,9 @@ types::global_dof_index DoFHandler::n_boundary_dofs () const template +template types::global_dof_index -DoFHandler::n_boundary_dofs (const FunctionMap &boundary_ids) const +DoFHandler::n_boundary_dofs (const std::map*> &boundary_ids) const { Assert (boundary_ids.find(numbers::internal_face_boundary_id) == boundary_ids.end(), ExcInvalidBoundaryIndicator()); diff --git a/source/dofs/dof_handler.inst.in b/source/dofs/dof_handler.inst.in index 8780b8cb08..17dcab09ff 100644 --- a/source/dofs/dof_handler.inst.in +++ b/source/dofs/dof_handler.inst.in @@ -14,6 +14,15 @@ // --------------------------------------------------------------------- +for (scalar: REAL_SCALARS; deal_II_dimension : DIMENSIONS) + { + template types::global_dof_index DoFHandler::n_boundary_dofs (const std::map*> &boundary_ids) const; + +#if deal_II_dimension < 3 + template types::global_dof_index DoFHandler::n_boundary_dofs (const std::map*> &boundary_ids) const; +#endif + + } for (deal_II_dimension : DIMENSIONS) { diff --git a/source/dofs/dof_tools_sparsity.cc b/source/dofs/dof_tools_sparsity.cc index 733f2189cd..eb98b99c72 100644 --- a/source/dofs/dof_tools_sparsity.cc +++ b/source/dofs/dof_tools_sparsity.cc @@ -338,7 +338,7 @@ namespace DoFTools { // there are only 2 boundary indicators in 1d, so it is no // performance problem to call the other function - typename DoFHandlerType::FunctionMap boundary_ids; + std::map*> boundary_ids; boundary_ids[0] = 0; boundary_ids[1] = 0; make_boundary_sparsity_pattern @@ -397,10 +397,10 @@ namespace DoFTools - template + template void make_boundary_sparsity_pattern (const DoFHandlerType &dof, - const typename FunctionMap::type &boundary_ids, + const std::map*> &boundary_ids, const std::vector &dof_to_boundary_mapping, SparsityPatternType &sparsity) { diff --git a/source/dofs/dof_tools_sparsity.inst.in b/source/dofs/dof_tools_sparsity.inst.in index dc2f35246c..2aa9dabccf 100644 --- a/source/dofs/dof_tools_sparsity.inst.in +++ b/source/dofs/dof_tools_sparsity.inst.in @@ -13,7 +13,47 @@ // // --------------------------------------------------------------------- +for (scalar: REAL_SCALARS; SP : SPARSITY_PATTERNS; deal_II_dimension : DIMENSIONS) + { + template void + DoFTools::make_boundary_sparsity_pattern,SP,scalar> + (const DoFHandler& dof, + const std::map*> &boundary_ids, + const std::vector &dof_to_boundary_mapping, + SP &sparsity); + + template void + DoFTools::make_boundary_sparsity_pattern,SP,scalar> + (const hp::DoFHandler& dof, + const std::map*> &boundary_ids, + const std::vector &dof_to_boundary_mapping, + SP &sparsity); + +#if deal_II_dimension < 3 + template void + DoFTools::make_boundary_sparsity_pattern,SP,scalar> + (const hp::DoFHandler& dof, + const std::map*> &boundary_ids, + const std::vector &dof_to_boundary_mapping, + SP &sparsity); + + template void + DoFTools::make_boundary_sparsity_pattern,SP,scalar> + (const DoFHandler& dof, + const std::map*> &boundary_ids, + const std::vector &dof_to_boundary_mapping, + SP &sparsity); + //template void + //DoFTools::make_boundary_sparsity_pattern,SP,scalar> + //(const hp::DoFHandler& dof, + // const std::map*> &boundary_ids, + // const std::vector &dof_to_boundary_mapping, + // SP &sparsity); + +#endif + + } for (SP : SPARSITY_PATTERNS; deal_II_dimension : DIMENSIONS) { @@ -75,29 +115,6 @@ for (SP : SPARSITY_PATTERNS; deal_II_dimension : DIMENSIONS) const std::vector &, SP &); - template void - DoFTools::make_boundary_sparsity_pattern,SP> - (const DoFHandler& dof, - const FunctionMap::type &boundary_ids, - const std::vector &dof_to_boundary_mapping, - SP &sparsity); - - template void - DoFTools::make_boundary_sparsity_pattern,SP> - (const hp::DoFHandler& dof, - const FunctionMap::type &boundary_ids, - const std::vector &dof_to_boundary_mapping, - SP &sparsity); - -#if deal_II_dimension < 3 - template void - DoFTools::make_boundary_sparsity_pattern,SP> - (const hp::DoFHandler& dof, - const FunctionMap::type &boundary_ids, - const std::vector &dof_to_boundary_mapping, - SP &sparsity); - #endif - template void DoFTools::make_flux_sparsity_pattern,SP> (const DoFHandler &dof, @@ -201,20 +218,6 @@ for (SP : SPARSITY_PATTERNS; deal_II_dimension : DIMENSIONS) // const std::vector &, // SP &); - template void - DoFTools::make_boundary_sparsity_pattern,SP> - (const DoFHandler& dof, - const FunctionMap::type &boundary_ids, - const std::vector &dof_to_boundary_mapping, - SP &sparsity); - - //template void - //DoFTools::make_boundary_sparsity_pattern,SP> - //(const hp::DoFHandler& dof, - // const FunctionMap::type &boundary_ids, - // const std::vector &dof_to_boundary_mapping, - // SP &sparsity); - #endif @@ -279,16 +282,16 @@ for (SP : SPARSITY_PATTERNS; deal_II_dimension : DIMENSIONS) SP &); template void - DoFTools::make_boundary_sparsity_pattern,SP> + DoFTools::make_boundary_sparsity_pattern,SP,double> (const DoFHandler<1,3>& dof, - const FunctionMap<3>::type &boundary_ids, + const std::map*> &boundary_ids, const std::vector &dof_to_boundary_mapping, SP &sparsity); template void - DoFTools::make_boundary_sparsity_pattern,SP> + DoFTools::make_boundary_sparsity_pattern,SP,double> (const hp::DoFHandler<1,3>& dof, - const FunctionMap<3>::type &boundary_ids, + const std::map*> &boundary_ids, const std::vector &dof_to_boundary_mapping, SP &sparsity); diff --git a/source/hp/dof_handler.cc b/source/hp/dof_handler.cc index 3279b07d34..da115ac68c 100644 --- a/source/hp/dof_handler.cc +++ b/source/hp/dof_handler.cc @@ -1832,14 +1832,15 @@ namespace hp template <> - types::global_dof_index DoFHandler<1>::n_boundary_dofs (const FunctionMap &boundary_ids) const + template + types::global_dof_index DoFHandler<1>::n_boundary_dofs (const std::map*> &boundary_ids) const { Assert (finite_elements != 0, ExcNoFESelected()); // check that only boundary // indicators 0 and 1 are allowed // in 1d - for (FunctionMap::const_iterator i=boundary_ids.begin(); + for (typename std::map*>::const_iterator i=boundary_ids.begin(); i!=boundary_ids.end(); ++i) Assert ((i->first == 0) || (i->first == 1), ExcInvalidBoundaryIndicator()); @@ -1916,7 +1917,8 @@ namespace hp } template <> - types::global_dof_index DoFHandler<1,2>::n_boundary_dofs (const FunctionMap &) const + template + types::global_dof_index DoFHandler<1,2>::n_boundary_dofs (const std::map*> &) const { Assert(false,ExcNotImplemented()); return 0; @@ -1939,7 +1941,8 @@ namespace hp } template <> - types::global_dof_index DoFHandler<1,3>::n_boundary_dofs (const FunctionMap &) const + template + types::global_dof_index DoFHandler<1,3>::n_boundary_dofs (const std::map*> &) const { Assert(false,ExcNotImplemented()); return 0; @@ -1992,8 +1995,9 @@ namespace hp template + template types::global_dof_index - DoFHandler::n_boundary_dofs (const FunctionMap &boundary_ids) const + DoFHandler::n_boundary_dofs (const std::map*> &boundary_ids) const { Assert (finite_elements != 0, ExcNoFESelected()); Assert (boundary_ids.find(numbers::internal_face_boundary_id) == boundary_ids.end(), @@ -2073,7 +2077,8 @@ namespace hp template <> - types::global_dof_index DoFHandler<2,3>::n_boundary_dofs (const FunctionMap &) const + template + types::global_dof_index DoFHandler<2,3>::n_boundary_dofs (const std::map*> &) const { Assert(false,ExcNotImplemented()); return 0; diff --git a/source/hp/dof_handler.inst.in b/source/hp/dof_handler.inst.in index a872a73729..dff1f33058 100644 --- a/source/hp/dof_handler.inst.in +++ b/source/hp/dof_handler.inst.in @@ -13,6 +13,18 @@ // // --------------------------------------------------------------------- +for (scalar: REAL_SCALARS; deal_II_dimension : DIMENSIONS) + { + namespace hp + \{ + template types::global_dof_index DoFHandler::n_boundary_dofs (const std::map*> &boundary_ids) const; + +#if deal_II_dimension < 3 + template types::global_dof_index DoFHandler::n_boundary_dofs (const std::map*> &boundary_ids) const; +#endif + \} + } + for (deal_II_dimension : DIMENSIONS) { diff --git a/source/lac/constraint_matrix.cc b/source/lac/constraint_matrix.cc index 214d602358..cfea9cf2cf 100644 --- a/source/lac/constraint_matrix.cc +++ b/source/lac/constraint_matrix.cc @@ -1283,6 +1283,7 @@ PARALLEL_VECTOR_FUNCTIONS(TrilinosWrappers::MPI::BlockVector); bool , \ internal::bool2type) const +MATRIX_FUNCTIONS(SparseMatrix); MATRIX_FUNCTIONS(SparseMatrix); MATRIX_FUNCTIONS(SparseMatrix); MATRIX_FUNCTIONS(FullMatrix); diff --git a/source/lac/constraint_matrix.inst.in b/source/lac/constraint_matrix.inst.in index 78c03f8d5f..4859a11632 100644 --- a/source/lac/constraint_matrix.inst.in +++ b/source/lac/constraint_matrix.inst.in @@ -18,7 +18,7 @@ for (S: REAL_SCALARS; T : DEAL_II_VEC_TEMPLATES) template void ConstraintMatrix::condense >(const T &, T &) const; template void ConstraintMatrix::condense >(T &vec) const; template void ConstraintMatrix::distribute_local_to_global > ( - const Vector&, const std::vector &, T &, const FullMatrix&) const; + const Vector&, const std::vector &, T &, const FullMatrix&) const; template void ConstraintMatrix::set_zero >(T &) const; } diff --git a/source/numerics/matrix_tools.cc b/source/numerics/matrix_tools.cc index ddd5358422..07a54505c5 100644 --- a/source/numerics/matrix_tools.cc +++ b/source/numerics/matrix_tools.cc @@ -137,13 +137,14 @@ namespace MatrixCreator namespace AssemblerData { template + int spacedim, + typename number> struct Scratch { Scratch (const ::dealii::hp::FECollection &fe, const UpdateFlags update_flags, - const Function *coefficient, - const Function *rhs_function, + const Function *coefficient, + const Function *rhs_function, const ::dealii::hp::QCollection &quadrature, const ::dealii::hp::MappingCollection &mapping) : @@ -156,10 +157,10 @@ namespace MatrixCreator update_flags), coefficient_values(quadrature_collection.max_n_quadrature_points()), coefficient_vector_values (quadrature_collection.max_n_quadrature_points(), - dealii::Vector (fe_collection.n_components())), + dealii::Vector (fe_collection.n_components())), rhs_values(quadrature_collection.max_n_quadrature_points()), rhs_vector_values(quadrature_collection.max_n_quadrature_points(), - dealii::Vector (fe_collection.n_components())), + dealii::Vector (fe_collection.n_components())), coefficient (coefficient), rhs_function (rhs_function), update_flags (update_flags) @@ -189,13 +190,13 @@ namespace MatrixCreator ::dealii::hp::FEValues x_fe_values; - std::vector coefficient_values; - std::vector > coefficient_vector_values; - std::vector rhs_values; - std::vector > rhs_vector_values; + std::vector coefficient_values; + std::vector > coefficient_vector_values; + std::vector rhs_values; + std::vector > rhs_vector_values; - const Function *coefficient; - const Function *rhs_function; + const Function *coefficient; + const Function *rhs_function; const UpdateFlags update_flags; }; @@ -217,7 +218,7 @@ namespace MatrixCreator typename CellIterator, typename number> void mass_assembler (const CellIterator &cell, - MatrixCreator::internal::AssemblerData::Scratch &data, + MatrixCreator::internal::AssemblerData::Scratch &data, MatrixCreator::internal::AssemblerData::CopyData ©_data) { data.x_fe_values.reinit (cell); @@ -255,7 +256,7 @@ namespace MatrixCreator else { data.rhs_vector_values.resize (n_q_points, - dealii::Vector(n_components)); + dealii::Vector(n_components)); data.rhs_function->vector_value_list (fe_values.get_quadrature_points(), data.rhs_vector_values); } @@ -273,14 +274,14 @@ namespace MatrixCreator else { data.coefficient_vector_values.resize (n_q_points, - dealii::Vector(n_components)); + dealii::Vector(n_components)); data.coefficient->vector_value_list (fe_values.get_quadrature_points(), data.coefficient_vector_values); } } - double add_data; + number add_data; const std::vector &JxW = fe_values.get_JxW_values(); for (unsigned int i=0; i void laplace_assembler (const CellIterator &cell, - MatrixCreator::internal::AssemblerData::Scratch &data, + MatrixCreator::internal::AssemblerData::Scratch &data, MatrixCreator::internal::AssemblerData::CopyData ©_data) { data.x_fe_values.reinit (cell); @@ -646,12 +647,12 @@ namespace MatrixCreator namespace MatrixCreator { - template + template void create_mass_matrix (const Mapping &mapping, const DoFHandler &dof, const Quadrature &q, SparseMatrix &matrix, - const Function *const coefficient, + const Function *const coefficient, const ConstraintMatrix &constraints) { Assert (matrix.m() == dof.n_dofs(), @@ -662,7 +663,7 @@ namespace MatrixCreator hp::FECollection fe_collection (dof.get_fe()); hp::QCollection q_collection (q); hp::MappingCollection mapping_collection (mapping); - MatrixCreator::internal::AssemblerData::Scratch + MatrixCreator::internal::AssemblerData::Scratch assembler_data (fe_collection, update_values | update_JxW_values | (coefficient != 0 ? update_quadrature_points : UpdateFlags(0)), @@ -688,11 +689,11 @@ namespace MatrixCreator - template + template void create_mass_matrix (const DoFHandler &dof, const Quadrature &q, SparseMatrix &matrix, - const Function *const coefficient, + const Function *const coefficient, const ConstraintMatrix &constraints) { create_mass_matrix(StaticMappingQ1::mapping, dof, @@ -701,14 +702,14 @@ namespace MatrixCreator - template + template void create_mass_matrix (const Mapping &mapping, const DoFHandler &dof, const Quadrature &q, SparseMatrix &matrix, - const Function &rhs, + const Function &rhs, Vector &rhs_vector, - const Function *const coefficient, + const Function *const coefficient, const ConstraintMatrix &constraints) { Assert (matrix.m() == dof.n_dofs(), @@ -719,7 +720,7 @@ namespace MatrixCreator hp::FECollection fe_collection (dof.get_fe()); hp::QCollection q_collection (q); hp::MappingCollection mapping_collection (mapping); - MatrixCreator::internal::AssemblerData::Scratch + MatrixCreator::internal::AssemblerData::Scratch assembler_data (fe_collection, update_values | update_JxW_values | update_quadrature_points, @@ -744,13 +745,13 @@ namespace MatrixCreator - template + template void create_mass_matrix (const DoFHandler &dof, const Quadrature &q, SparseMatrix &matrix, - const Function &rhs, + const Function &rhs, Vector &rhs_vector, - const Function *const coefficient, + const Function *const coefficient, const ConstraintMatrix &constraints) { create_mass_matrix(StaticMappingQ1::mapping, @@ -760,12 +761,12 @@ namespace MatrixCreator - template + template void create_mass_matrix (const hp::MappingCollection &mapping, const hp::DoFHandler &dof, const hp::QCollection &q, SparseMatrix &matrix, - const Function *const coefficient, + const Function *const coefficient, const ConstraintMatrix &constraints) { Assert (matrix.m() == dof.n_dofs(), @@ -773,7 +774,7 @@ namespace MatrixCreator Assert (matrix.n() == dof.n_dofs(), ExcDimensionMismatch (matrix.n(), dof.n_dofs())); - MatrixCreator::internal::AssemblerData::Scratch + MatrixCreator::internal::AssemblerData::Scratch assembler_data (dof.get_fe(), update_values | update_JxW_values | (coefficient != 0 ? update_quadrature_points : UpdateFlags(0)), @@ -798,11 +799,11 @@ namespace MatrixCreator - template + template void create_mass_matrix (const hp::DoFHandler &dof, const hp::QCollection &q, SparseMatrix &matrix, - const Function *const coefficient, + const Function *const coefficient, const ConstraintMatrix &constraints) { create_mass_matrix(hp::StaticMappingQ1::mapping_collection, @@ -811,14 +812,14 @@ namespace MatrixCreator - template + template void create_mass_matrix (const hp::MappingCollection &mapping, const hp::DoFHandler &dof, const hp::QCollection &q, SparseMatrix &matrix, - const Function &rhs, + const Function &rhs, Vector &rhs_vector, - const Function *const coefficient, + const Function *const coefficient, const ConstraintMatrix &constraints) { Assert (matrix.m() == dof.n_dofs(), @@ -826,7 +827,7 @@ namespace MatrixCreator Assert (matrix.n() == dof.n_dofs(), ExcDimensionMismatch (matrix.n(), dof.n_dofs())); - MatrixCreator::internal::AssemblerData::Scratch + MatrixCreator::internal::AssemblerData::Scratch assembler_data (dof.get_fe(), update_values | update_JxW_values | update_quadrature_points, @@ -851,13 +852,13 @@ namespace MatrixCreator - template + template void create_mass_matrix (const hp::DoFHandler &dof, const hp::QCollection &q, SparseMatrix &matrix, - const Function &rhs, + const Function &rhs, Vector &rhs_vector, - const Function *const coefficient, + const Function *const coefficient, const ConstraintMatrix &constraints) { create_mass_matrix(hp::StaticMappingQ1::mapping_collection, dof, q, @@ -868,7 +869,7 @@ namespace MatrixCreator namespace { - template + template void create_boundary_mass_matrix_1 (typename DoFHandler::active_cell_iterator const &cell, MatrixCreator::internal::AssemblerBoundary::Scratch const &, @@ -877,8 +878,8 @@ namespace MatrixCreator Mapping const &mapping, FiniteElement const &fe, Quadrature const &q, - typename FunctionMap::type const &boundary_functions, - Function const *const coefficient, + std::map*> const &boundary_functions, + Function const *const coefficient, std::vector const &component_mapping) { @@ -904,15 +905,15 @@ namespace MatrixCreator // 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)); + 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 rhs_values_scalar (fe_values.n_quadrature_points); + std::vector > rhs_values_system (fe_values.n_quadrature_points, + Vector(n_function_components)); copy_data.dofs.resize(copy_data.dofs_per_cell); cell->get_dof_indices (copy_data.dofs); @@ -1071,13 +1072,13 @@ namespace MatrixCreator } } - template + template void copy_boundary_mass_matrix_1(MatrixCreator::internal::AssemblerBoundary::CopyData > const ©_data, - typename FunctionMap::type const &boundary_functions, + std::map*> const &boundary_functions, std::vector const &dof_to_boundary_mapping, - SparseMatrix &matrix, - Vector &rhs_vector) + SparseMatrix &matrix, + Vector &rhs_vector) { // now transfer cell matrix and vector to the whole boundary matrix // @@ -1139,34 +1140,67 @@ namespace MatrixCreator template <> void - create_boundary_mass_matrix_1<1,3> (DoFHandler<1,3>::active_cell_iterator const &/*cell*/, - MatrixCreator::internal::AssemblerBoundary::Scratch const &, - MatrixCreator::internal::AssemblerBoundary::CopyData > &/*copy_data*/, - Mapping<1,3> const &, - FiniteElement<1,3> const &, - Quadrature<0> const &, - FunctionMap<3>::type const &/*boundary_functions*/, - Function<3> const *const /*coefficient*/, - std::vector const &/*component_mapping*/) + create_boundary_mass_matrix_1<1,3,float> (DoFHandler<1,3>::active_cell_iterator const &/*cell*/, + MatrixCreator::internal::AssemblerBoundary::Scratch const &, + MatrixCreator::internal::AssemblerBoundary::CopyData > &/*copy_data*/, + Mapping<1,3> const &, + FiniteElement<1,3> const &, + Quadrature<0> const &, + std::map*> const &/*boundary_functions*/, + Function<3,float> const *const /*coefficient*/, + std::vector const &/*component_mapping*/) + { + Assert(false,ExcNotImplemented()); + } + + template <> + void + create_boundary_mass_matrix_1<1,3,double> (DoFHandler<1,3>::active_cell_iterator const &/*cell*/, + MatrixCreator::internal::AssemblerBoundary::Scratch const &, + MatrixCreator::internal::AssemblerBoundary::CopyData > &/*copy_data*/, + Mapping<1,3> const &, + FiniteElement<1,3> const &, + Quadrature<0> const &, + std::map*> const &/*boundary_functions*/, + Function<3,double> const *const /*coefficient*/, + std::vector const &/*component_mapping*/) + { + Assert(false,ExcNotImplemented()); + } + + template <> + void + create_boundary_mass_matrix_1<1,3,long double> (DoFHandler<1,3>::active_cell_iterator const &/*cell*/, + MatrixCreator::internal::AssemblerBoundary::Scratch const &, + MatrixCreator::internal::AssemblerBoundary::CopyData > &/*copy_data*/, + Mapping<1,3> const &, + FiniteElement<1,3> const &, + Quadrature<0> const &, + std::map*> const &/*boundary_functions*/, + Function<3,long double> const *const /*coefficient*/, + std::vector const &/*component_mapping*/) { Assert(false,ExcNotImplemented()); } + } - template + 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, + SparseMatrix &matrix, + const std::map*> &boundary_functions, + Vector &rhs_vector, std::vector &dof_to_boundary_mapping, - const Function *const coefficient, + const Function *const coefficient, std::vector component_mapping) { // what would that be in 1d? the @@ -1181,7 +1215,7 @@ namespace MatrixCreator const FiniteElement &fe = dof.get_fe(); const unsigned int n_components = fe.n_components(); - Assert (matrix.n() == dof.n_boundary_dofs(boundary_functions), + Assert (matrix.n() == dof.n_boundary_dofs (boundary_functions), ExcInternalError()); Assert (matrix.n() == matrix.m(), ExcInternalError()); Assert (matrix.n() == rhs_vector.size(), ExcInternalError()); @@ -1208,14 +1242,14 @@ namespace MatrixCreator static_cast::active_cell_iterator const &,MatrixCreator::internal::AssemblerBoundary::Scratch const &, MatrixCreator::internal::AssemblerBoundary::CopyData > &)> > - (std_cxx11::bind(&create_boundary_mass_matrix_1,std_cxx11::_1,std_cxx11::_2, + (std_cxx11::bind(&create_boundary_mass_matrix_1,std_cxx11::_1,std_cxx11::_2, std_cxx11::_3, std_cxx11::cref(mapping),std_cxx11::cref(fe),std_cxx11::cref(q), std_cxx11::cref(boundary_functions),coefficient, std_cxx11::cref(component_mapping))), static_cast > const &)> > (std_cxx11::bind( - ©_boundary_mass_matrix_1, + ©_boundary_mass_matrix_1, std_cxx11::_1, std_cxx11::cref(boundary_functions), std_cxx11::cref(dof_to_boundary_mapping), @@ -1230,7 +1264,7 @@ namespace MatrixCreator namespace { - template + template void create_hp_boundary_mass_matrix_1 (typename hp::DoFHandler::active_cell_iterator const &cell, @@ -1240,8 +1274,8 @@ namespace MatrixCreator hp::MappingCollection const &mapping, hp::FECollection const &fe_collection, hp::QCollection const &q, - const typename FunctionMap::type &boundary_functions, - Function const *const coefficient, + const std::map*> &boundary_functions, + Function const *const coefficient, std::vector const &component_mapping) { const unsigned int n_components = fe_collection.n_components(); @@ -1264,11 +1298,11 @@ namespace MatrixCreator // 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 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_on_face_vector (dofs_per_face); @@ -1301,7 +1335,7 @@ namespace MatrixCreator // FE has several components { rhs_values_system.resize (fe_values.n_quadrature_points, - Vector(n_function_components)); + Vector(n_function_components)); boundary_functions.find(cell->face(face)->boundary_id()) ->second->vector_value_list (fe_values.get_quadrature_points(), rhs_values_system); @@ -1337,7 +1371,7 @@ namespace MatrixCreator else { coefficient_vector_values.resize (fe_values.n_quadrature_points, - Vector(n_components)); + Vector(n_components)); coefficient->vector_value_list (fe_values.get_quadrature_points(), coefficient_vector_values); for (unsigned int point=0; point + template void copy_hp_boundary_mass_matrix_1(MatrixCreator::internal::AssemblerBoundary ::CopyData > const ©_data, - typename FunctionMap::type const &boundary_functions, + std::map*> const &boundary_functions, std::vector const &dof_to_boundary_mapping, - SparseMatrix &matrix, - Vector &rhs_vector) + SparseMatrix &matrix, + Vector &rhs_vector) { // now transfer cell matrix and vector to the whole boundary matrix // @@ -1556,14 +1590,14 @@ namespace MatrixCreator - template + template void create_boundary_mass_matrix (const DoFHandler &dof, const Quadrature &q, - SparseMatrix &matrix, - const typename FunctionMap::type &rhs, - Vector &rhs_vector, + SparseMatrix &matrix, + const std::map*> &rhs, + Vector &rhs_vector, std::vector &dof_to_boundary_mapping, - const Function *const a, + const Function *const a, std::vector component_mapping) { create_boundary_mass_matrix(StaticMappingQ1::mapping, dof, q, @@ -1572,16 +1606,16 @@ namespace MatrixCreator - template + 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, + SparseMatrix &matrix, + const std::map*> &boundary_functions, + Vector &rhs_vector, std::vector &dof_to_boundary_mapping, - const Function *const coefficient, + const Function *const coefficient, std::vector component_mapping) { // what would that be in 1d? the @@ -1596,7 +1630,7 @@ namespace MatrixCreator 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), + Assert (matrix.n() == dof.n_boundary_dofs (boundary_functions), ExcInternalError()); Assert (matrix.n() == matrix.m(), ExcInternalError()); Assert (matrix.n() == rhs_vector.size(), ExcInternalError()); @@ -1623,14 +1657,14 @@ namespace MatrixCreator static_cast::active_cell_iterator const &,MatrixCreator::internal::AssemblerBoundary::Scratch const &, MatrixCreator::internal::AssemblerBoundary::CopyData > &)> > - (std_cxx11::bind( &create_hp_boundary_mass_matrix_1,std_cxx11::_1,std_cxx11::_2, + (std_cxx11::bind( &create_hp_boundary_mass_matrix_1,std_cxx11::_1,std_cxx11::_2, std_cxx11::_3, std_cxx11::cref(mapping),std_cxx11::cref(fe_collection),std_cxx11::cref(q), std_cxx11::cref(boundary_functions),coefficient, std_cxx11::cref(component_mapping))), static_cast > const &)> > ( - std_cxx11::bind( ©_hp_boundary_mass_matrix_1, + std_cxx11::bind( ©_hp_boundary_mass_matrix_1, std_cxx11::_1, std_cxx11::cref(boundary_functions), std_cxx11::cref(dof_to_boundary_mapping), @@ -1643,14 +1677,14 @@ namespace MatrixCreator - template + template void create_boundary_mass_matrix (const hp::DoFHandler &dof, const hp::QCollection &q, - SparseMatrix &matrix, - const typename FunctionMap::type &rhs, - Vector &rhs_vector, + SparseMatrix &matrix, + const std::map*> &rhs, + Vector &rhs_vector, std::vector &dof_to_boundary_mapping, - const Function *const a, + const Function *const a, std::vector component_mapping) { create_boundary_mass_matrix(hp::StaticMappingQ1::mapping_collection, dof, q, @@ -1675,7 +1709,7 @@ namespace MatrixCreator hp::FECollection fe_collection (dof.get_fe()); hp::QCollection q_collection (q); hp::MappingCollection mapping_collection (mapping); - MatrixCreator::internal::AssemblerData::Scratch + MatrixCreator::internal::AssemblerData::Scratch assembler_data (fe_collection, update_gradients | update_JxW_values | (coefficient != 0 ? update_quadrature_points : UpdateFlags(0)), @@ -1733,7 +1767,7 @@ namespace MatrixCreator hp::FECollection fe_collection (dof.get_fe()); hp::QCollection q_collection (q); hp::MappingCollection mapping_collection (mapping); - MatrixCreator::internal::AssemblerData::Scratch + MatrixCreator::internal::AssemblerData::Scratch assembler_data (fe_collection, update_gradients | update_values | update_JxW_values | update_quadrature_points, @@ -1788,7 +1822,7 @@ namespace MatrixCreator Assert (matrix.n() == dof.n_dofs(), ExcDimensionMismatch (matrix.n(), dof.n_dofs())); - MatrixCreator::internal::AssemblerData::Scratch + MatrixCreator::internal::AssemblerData::Scratch assembler_data (dof.get_fe(), update_gradients | update_JxW_values | (coefficient != 0 ? update_quadrature_points : UpdateFlags(0)), @@ -1843,7 +1877,7 @@ namespace MatrixCreator Assert (matrix.n() == dof.n_dofs(), ExcDimensionMismatch (matrix.n(), dof.n_dofs())); - MatrixCreator::internal::AssemblerData::Scratch + MatrixCreator::internal::AssemblerData::Scratch assembler_data (dof.get_fe(), update_gradients | update_values | update_JxW_values | update_quadrature_points, @@ -1901,7 +1935,7 @@ namespace MatrixTools //TODO:[WB] I don't think that the optimized storage of diagonals is needed (GK) template void - apply_boundary_values (const std::map &boundary_values, + apply_boundary_values (const std::map &boundary_values, SparseMatrix &matrix, Vector &solution, Vector &right_hand_side, @@ -1937,8 +1971,8 @@ namespace MatrixTools } - std::map::const_iterator dof = boundary_values.begin(), - endd = boundary_values.end(); + typename std::map::const_iterator dof = boundary_values.begin(), + endd = boundary_values.end(); for (; dof != endd; ++dof) { Assert (dof->first < n_dofs, ExcInternalError()); @@ -2059,7 +2093,7 @@ namespace MatrixTools template void - apply_boundary_values (const std::map &boundary_values, + apply_boundary_values (const std::map &boundary_values, BlockSparseMatrix &matrix, BlockVector &solution, BlockVector &right_hand_side, @@ -2119,8 +2153,8 @@ namespace MatrixTools first_nonzero_diagonal_entry = 1; - std::map::const_iterator dof = boundary_values.begin(), - endd = boundary_values.end(); + typename std::map::const_iterator dof = boundary_values.begin(), + endd = boundary_values.end(); const BlockSparsityPattern & sparsity_pattern = matrix.get_sparsity_pattern(); @@ -2324,7 +2358,7 @@ namespace MatrixTools { template void - apply_boundary_values (const std::map &boundary_values, + apply_boundary_values (const std::map &boundary_values, PETScMatrix &matrix, PETScVector &solution, PETScVector &right_hand_side, @@ -2425,7 +2459,7 @@ namespace MatrixTools void - apply_boundary_values (const std::map &boundary_values, + apply_boundary_values (const std::map &boundary_values, PETScWrappers::SparseMatrix &matrix, PETScWrappers::Vector &solution, PETScWrappers::Vector &right_hand_side, @@ -2440,7 +2474,7 @@ namespace MatrixTools void - apply_boundary_values (const std::map &boundary_values, + apply_boundary_values (const std::map &boundary_values, PETScWrappers::MPI::SparseMatrix &matrix, PETScWrappers::MPI::Vector &solution, PETScWrappers::MPI::Vector &right_hand_side, @@ -2454,7 +2488,7 @@ namespace MatrixTools void - apply_boundary_values (const std::map &boundary_values, + apply_boundary_values (const std::map &boundary_values, PETScWrappers::MPI::BlockSparseMatrix &matrix, PETScWrappers::MPI::BlockVector &solution, PETScWrappers::MPI::BlockVector &right_hand_side, @@ -2473,11 +2507,11 @@ namespace MatrixTools // into blocks for the boundary values. // To this end, generate a vector of // maps with the respective indices. - std::vector > block_boundary_values(n_blocks); + std::vector > block_boundary_values(n_blocks); { int block = 0; dealii::types::global_dof_index offset = 0; - for (std::map::const_iterator + for (std::map::const_iterator dof = boundary_values.begin(); dof != boundary_values.end(); ++dof) @@ -2512,7 +2546,7 @@ namespace MatrixTools = matrix.block(block_m,0).local_range(); std::vector constrained_rows; - for (std::map::const_iterator + for (std::map::const_iterator dof = block_boundary_values[block_m].begin(); dof != block_boundary_values[block_m].end(); ++dof) @@ -2538,7 +2572,7 @@ namespace MatrixTools { template void - apply_boundary_values (const std::map &boundary_values, + apply_boundary_values (const std::map &boundary_values, TrilinosMatrix &matrix, TrilinosVector &solution, TrilinosVector &right_hand_side, @@ -2579,7 +2613,7 @@ namespace MatrixTools // figure out which rows of the matrix we // have to eliminate on this processor std::vector constrained_rows; - for (std::map::const_iterator + for (std::map::const_iterator dof = boundary_values.begin(); dof != boundary_values.end(); ++dof) @@ -2598,7 +2632,7 @@ namespace MatrixTools std::vector indices; std::vector solution_values; - for (std::map::const_iterator + for (std::map::const_iterator dof = boundary_values.begin(); dof != boundary_values.end(); ++dof) @@ -2635,7 +2669,7 @@ namespace MatrixTools template void - apply_block_boundary_values (const std::map &boundary_values, + apply_block_boundary_values (const std::map &boundary_values, TrilinosMatrix &matrix, TrilinosBlockVector &solution, TrilinosBlockVector &right_hand_side, @@ -2656,11 +2690,11 @@ namespace MatrixTools // into blocks for the boundary values. // To this end, generate a vector of // maps with the respective indices. - std::vector > block_boundary_values(n_blocks); + std::vector > block_boundary_values(n_blocks); { int block=0; types::global_dof_index offset = 0; - for (std::map::const_iterator + for (std::map::const_iterator dof = boundary_values.begin(); dof != boundary_values.end(); ++dof) @@ -2672,7 +2706,7 @@ namespace MatrixTools } const types::global_dof_index index = dof->first - offset; block_boundary_values[block].insert( - std::pair (index,dof->second)); + std::pair (index,dof->second)); } } @@ -2696,7 +2730,7 @@ namespace MatrixTools = matrix.block(block_m,0).local_range(); std::vector constrained_rows; - for (std::map::const_iterator + for (std::map::const_iterator dof = block_boundary_values[block_m].begin(); dof != block_boundary_values[block_m].end(); ++dof) @@ -2716,7 +2750,7 @@ namespace MatrixTools void - apply_boundary_values (const std::map &boundary_values, + apply_boundary_values (const std::map &boundary_values, TrilinosWrappers::SparseMatrix &matrix, TrilinosWrappers::Vector &solution, TrilinosWrappers::Vector &right_hand_side, @@ -2731,7 +2765,7 @@ namespace MatrixTools void - apply_boundary_values (const std::map &boundary_values, + apply_boundary_values (const std::map &boundary_values, TrilinosWrappers::SparseMatrix &matrix, TrilinosWrappers::MPI::Vector &solution, TrilinosWrappers::MPI::Vector &right_hand_side, @@ -2746,7 +2780,7 @@ namespace MatrixTools void - apply_boundary_values (const std::map &boundary_values, + apply_boundary_values (const std::map &boundary_values, TrilinosWrappers::BlockSparseMatrix &matrix, TrilinosWrappers::BlockVector &solution, TrilinosWrappers::BlockVector &right_hand_side, @@ -2760,7 +2794,7 @@ namespace MatrixTools void - apply_boundary_values (const std::map &boundary_values, + apply_boundary_values (const std::map &boundary_values, TrilinosWrappers::BlockSparseMatrix &matrix, TrilinosWrappers::MPI::BlockVector &solution, TrilinosWrappers::MPI::BlockVector &right_hand_side, @@ -2774,12 +2808,12 @@ namespace MatrixTools #endif - + template void - local_apply_boundary_values (const std::map &boundary_values, + local_apply_boundary_values (const std::map &boundary_values, const std::vector &local_dof_indices, - FullMatrix &local_matrix, - Vector &local_rhs, + FullMatrix &local_matrix, + Vector &local_rhs, const bool eliminate_columns) { Assert (local_dof_indices.size() == local_matrix.m(), @@ -2820,11 +2854,11 @@ namespace MatrixTools // // we only compute this value lazily the // first time we need it. - double average_diagonal = 0; + number average_diagonal = 0; const unsigned int n_local_dofs = local_dof_indices.size(); for (unsigned int i=0; i::const_iterator + const typename std::map::const_iterator boundary_value = boundary_values.find (local_dof_indices[i]); if (boundary_value != boundary_values.end()) { @@ -2903,14 +2937,31 @@ namespace MatrixTools Vector &solution, Vector &right_hand_side, const bool eliminate_columns); + + template + void + apply_boundary_values (const std::map &boundary_values, + SparseMatrix &matrix, + Vector &solution, + Vector &right_hand_side, + const bool eliminate_columns); + template void - apply_boundary_values (const std::map &boundary_values, + apply_boundary_values (const std::map &boundary_values, SparseMatrix &matrix, Vector &solution, Vector &right_hand_side, const bool eliminate_columns); + template + void + apply_boundary_values (const std::map &boundary_values, + BlockSparseMatrix &matrix, + BlockVector &solution, + BlockVector &right_hand_side, + const bool eliminate_columns); + template void apply_boundary_values (const std::map &boundary_values, @@ -2920,7 +2971,7 @@ namespace MatrixTools const bool eliminate_columns); template void - apply_boundary_values (const std::map &boundary_values, + apply_boundary_values (const std::map &boundary_values, BlockSparseMatrix &matrix, BlockVector &solution, BlockVector &right_hand_side, diff --git a/source/numerics/matrix_tools.inst.in b/source/numerics/matrix_tools.inst.in index a2e9346dc5..8787c14630 100644 --- a/source/numerics/matrix_tools.inst.in +++ b/source/numerics/matrix_tools.inst.in @@ -13,184 +13,170 @@ // // --------------------------------------------------------------------- -for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) +for (number: REAL_SCALARS) { -#if deal_II_dimension <= deal_II_space_dimension + template + void MatrixTools::local_apply_boundary_values + (const std::map &boundary_values, + const std::vector &local_dof_indices, + FullMatrix &local_matrix, + Vector &local_rhs, + const bool eliminate_columns); -// non-hp version of create_mass_matrix + } + +for (scalar: REAL_SCALARS; deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) + { +#if deal_II_dimension <= deal_II_space_dimension + // non-hp version of create_mass_matrix template - void MatrixCreator::create_mass_matrix + void MatrixCreator::create_mass_matrix (const Mapping &mapping, const DoFHandler &dof, const Quadrature &q, - SparseMatrix &matrix, - const Function * const coefficient, + SparseMatrix &matrix, + const Function * const coefficient, const ConstraintMatrix &constraints); template - void MatrixCreator::create_mass_matrix + void MatrixCreator::create_mass_matrix (const DoFHandler &dof, const Quadrature &q, - SparseMatrix &matrix, - const Function * const coefficient, + SparseMatrix &matrix, + const Function * const coefficient, const ConstraintMatrix &constraints); template - void MatrixCreator::create_mass_matrix + 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, + SparseMatrix &matrix, + const Function &rhs, + Vector &rhs_vector, + const Function * const coefficient, const ConstraintMatrix &constraints); template - void MatrixCreator::create_mass_matrix + void MatrixCreator::create_mass_matrix (const DoFHandler &dof, const Quadrature &q, - SparseMatrix &matrix, - const Function &rhs, - Vector &rhs_vector, - const Function * const coefficient, + SparseMatrix &matrix, + const Function &rhs, + Vector &rhs_vector, + const Function * const coefficient, const ConstraintMatrix &constraints); +#endif + } +for (scalar: REAL_SCALARS; deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) + { +#if deal_II_dimension <= deal_II_space_dimension template - void MatrixCreator::create_boundary_mass_matrix + void MatrixCreator::create_boundary_mass_matrix (const DoFHandler &dof, const Quadrature &q, - SparseMatrix &matrix, - const FunctionMap::type &rhs, - Vector &rhs_vector, + SparseMatrix &matrix, + const std::map*> &rhs, + Vector &rhs_vector, std::vector &dof_to_boundary_mapping, - const Function * const a, + const Function * const a, std::vector); template - void MatrixCreator::create_boundary_mass_matrix + void MatrixCreator::create_boundary_mass_matrix (const Mapping &, const DoFHandler &dof, const Quadrature &q, - SparseMatrix &matrix, - const FunctionMap::type &rhs, - Vector &rhs_vector, + SparseMatrix &matrix, + const std::map*> &rhs, + Vector &rhs_vector, std::vector &dof_to_boundary_mapping, - const Function * const a, + const Function * const a, std::vector); template void - MatrixCreator::create_boundary_mass_matrix + MatrixCreator::create_boundary_mass_matrix (const hp::MappingCollection&, const hp::DoFHandler&, const hp::QCollection&, - SparseMatrix&, - const FunctionMap::type&, - Vector&, + SparseMatrix&, + const std::map*>&, + Vector&, std::vector&, - const Function * const, + const Function * const, std::vector); template - void MatrixCreator::create_boundary_mass_matrix + void MatrixCreator::create_boundary_mass_matrix (const hp::DoFHandler&, const hp::QCollection&, - SparseMatrix&, - const FunctionMap::type&, - Vector&, + SparseMatrix&, + const std::map*>&, + Vector&, std::vector&, - const Function * const, + const Function * const, std::vector); - -// same for float - template - void MatrixCreator::create_mass_matrix - (const Mapping &mapping, - const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - const Function * const coefficient, - const ConstraintMatrix &constraints); - template - void MatrixCreator::create_mass_matrix - (const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - const Function * const coefficient, - const ConstraintMatrix &constraints); - 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, - const ConstraintMatrix &constraints); - template - void MatrixCreator::create_mass_matrix - (const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - const Function &rhs, - Vector &rhs_vector, - const Function * const coefficient, - const ConstraintMatrix &constraints); - #endif } //TODO[SP]: replace by // where applicable and move to codimension cases above also when applicable -for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) +for (scalar : REAL_SCALARS; deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) { // hp versions of functions #if deal_II_dimension <= deal_II_space_dimension template - void MatrixCreator::create_mass_matrix + void MatrixCreator::create_mass_matrix (const hp::MappingCollection &mapping, const hp::DoFHandler &dof, const hp::QCollection &q, - SparseMatrix &matrix, - const Function * const coefficient, + SparseMatrix &matrix, + const Function * const coefficient, const ConstraintMatrix &constraints); template - void MatrixCreator::create_mass_matrix + 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, + SparseMatrix &matrix, + const Function &rhs, + Vector &rhs_vector, + const Function * const coefficient, const ConstraintMatrix &constraints); - #endif #if deal_II_dimension == deal_II_space_dimension template - void MatrixCreator::create_mass_matrix + void MatrixCreator::create_mass_matrix (const hp::DoFHandler &dof, const hp::QCollection &q, - SparseMatrix &matrix, - const Function * const coefficient, + SparseMatrix &matrix, + const Function * const coefficient, const ConstraintMatrix &constraints); template - void MatrixCreator::create_mass_matrix + 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, + SparseMatrix &matrix, + const Function &rhs, + Vector &rhs_vector, + const Function * const coefficient, const ConstraintMatrix &constraints); +#endif + } + + +for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) + { +#if deal_II_dimension == deal_II_space_dimension + // non-hp versions of create_laplace_matrix template void MatrixCreator::create_laplace_matrix @@ -264,55 +250,5 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS const ConstraintMatrix &constraints); #endif - -// same for float -#if deal_II_dimension <= deal_II_space_dimension - - template - void MatrixCreator::create_mass_matrix - (const hp::MappingCollection &mapping, - const hp::DoFHandler &dof, - const hp::QCollection &q, - SparseMatrix &matrix, - const Function * const coefficient, - const ConstraintMatrix &constraints); - - 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, - const ConstraintMatrix &constraints); - - -#endif - -#if deal_II_dimension == deal_II_space_dimension - - - template - void MatrixCreator::create_mass_matrix - (const hp::DoFHandler &dof, - const hp::QCollection &q, - SparseMatrix &matrix, - const Function * const coefficient, - const ConstraintMatrix &constraints); - - 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, - const ConstraintMatrix &constraints); - -#endif - } diff --git a/source/numerics/vector_tools_boundary.inst.in b/source/numerics/vector_tools_boundary.inst.in index c23601adae..b5fef1b0e5 100644 --- a/source/numerics/vector_tools_boundary.inst.in +++ b/source/numerics/vector_tools_boundary.inst.in @@ -22,7 +22,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS; void interpolate_boundary_values (const Mapping &, const DH &, - const FunctionMap::type &, + const std::map*> &, std::map &, const ComponentMask &); @@ -46,7 +46,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS; template void interpolate_boundary_values (const DH &, - const FunctionMap::type &, + const std::map*> &, std::map &, const ComponentMask &); @@ -54,7 +54,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS; void interpolate_boundary_values ( const Mapping &, const DH &, - const FunctionMap::type &, + const std::map*> &, ConstraintMatrix &, const ComponentMask &); @@ -78,7 +78,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS; template void interpolate_boundary_values ( const DH &, - const FunctionMap::type &, + const std::map*> &, ConstraintMatrix &, const ComponentMask &); \} @@ -94,14 +94,14 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) void project_boundary_values (const Mapping &, const DoFHandler &, - const FunctionMap::type &, + const std::map*> &, const Quadrature&, std::map&, std::vector); template void project_boundary_values (const DoFHandler &, - const FunctionMap::type &, + const std::map*> &, const Quadrature&, std::map&, std::vector); @@ -110,14 +110,14 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) void project_boundary_values (const Mapping &, const DoFHandler &, - const FunctionMap::type &, + const std::map*> &, const Quadrature&, ConstraintMatrix&, std::vector); template void project_boundary_values (const DoFHandler &, - const FunctionMap::type &, + const std::map*> &, const Quadrature&, ConstraintMatrix&, std::vector); @@ -125,7 +125,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) template void project_boundary_values (const hp::DoFHandler &, - const FunctionMap::type &, + const std::map*> &, const hp::QCollection &, std::map &, std::vector); diff --git a/source/numerics/vector_tools_project.inst.in b/source/numerics/vector_tools_project.inst.in index 2afeedbb2e..ffbdda0908 100644 --- a/source/numerics/vector_tools_project.inst.in +++ b/source/numerics/vector_tools_project.inst.in @@ -25,7 +25,7 @@ for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimens const DoFHandler &, const ConstraintMatrix &, const Quadrature &, - const Function &, + const Function &, VEC &, const bool, const Quadrature &, @@ -36,7 +36,7 @@ for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimens (const DoFHandler &, const ConstraintMatrix &, const Quadrature &, - const Function &, + const Function &, VEC &, const bool, const Quadrature &, @@ -48,7 +48,7 @@ for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimens const hp::DoFHandler &, const ConstraintMatrix &, const hp::QCollection &, - const Function &, + const Function &, VEC &, const bool, const hp::QCollection &, @@ -59,7 +59,7 @@ for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimens (const hp::DoFHandler &, const ConstraintMatrix &, const hp::QCollection &, - const Function &, + const Function &, VEC &, const bool, const hp::QCollection &,