From: Fabian Castelli Date: Wed, 6 Jul 2022 14:38:22 +0000 (+0200) Subject: Clean up of VectorTools::project functions X-Git-Tag: v9.5.0-rc1~1054^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=900e3ebaaedf7f49c122c0a63fb666d5d9485cb4;p=dealii.git Clean up of VectorTools::project functions --- diff --git a/include/deal.II/numerics/vector_tools_project.templates.h b/include/deal.II/numerics/vector_tools_project.templates.h index bd8497e8b6..f99d25cd06 100644 --- a/include/deal.II/numerics/vector_tools_project.templates.h +++ b/include/deal.II/numerics/vector_tools_project.templates.h @@ -52,11 +52,11 @@ 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 DoFHandler & dof_handler, - std::map &boundary_values) + std::map &boundary_values) { // loop over all boundary faces // to get all dof indices of @@ -79,11 +79,8 @@ namespace VectorTools // that is actually wholly on // the boundary, not only by // one line or one vertex - typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); std::vector face_dof_indices; - for (; cell != endc; ++cell) + for (const auto &cell : dof_handler.active_cell_iterators()) for (auto f : GeometryInfo::face_indices()) if (cell->at_boundary(f)) { @@ -102,6 +99,8 @@ namespace VectorTools } } + + /** * Compute the boundary values to be used in the project() functions. */ @@ -111,16 +110,16 @@ namespace VectorTools class M_or_MC, template class Q_or_QC, - typename number> + typename Number> void project_compute_b_v( const M_or_MC & mapping, const DoFHandler & 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 @@ -143,7 +142,7 @@ namespace VectorTools const std::vector used_boundary_ids = dof.get_triangulation().get_boundary_ids(); - std::map *> + std::map *> boundary_functions; for (const auto used_boundary_id : used_boundary_ids) boundary_functions[used_boundary_id] = &function; @@ -153,6 +152,7 @@ namespace VectorTools } + /* * MatrixFree implementation of project() for an arbitrary number of * components of the FiniteElement. @@ -180,12 +180,8 @@ namespace VectorTools (void)q_boundary; AssertDimension(dof.get_fe_collection().size(), 1); - - Assert(dof.get_fe(0).n_components() == function.n_components, - ExcDimensionMismatch(dof.get_fe(0).n_components(), - function.n_components)); - Assert(dof.get_fe(0).n_components() == components, - ExcDimensionMismatch(components, dof.get_fe(0).n_components())); + AssertDimension(dof.get_fe(0).n_components(), function.n_components); + AssertDimension(dof.get_fe(0).n_components(), components); Quadrature quadrature_mf; @@ -329,8 +325,10 @@ namespace VectorTools - // Helper interface for the matrix-free implementation of project(). - // Used to determine the number of components. + /** + * Helper interface for the matrix-free implementation of project(). Used + * to determine the number of components. + */ template void project_matrix_free_component( @@ -422,8 +420,7 @@ namespace VectorTools const Quadrature &q_boundary, const bool project_to_boundary_first) { - Assert(vec_result.size() == dof.n_dofs(), - ExcDimensionMismatch(vec_result.size(), dof.n_dofs())); + AssertDimension(vec_result.size(), dof.n_dofs()); LinearAlgebra::distributed::Vector work_result; @@ -437,11 +434,9 @@ namespace VectorTools q_boundary, project_to_boundary_first); - const IndexSet & locally_owned_dofs = dof.locally_owned_dofs(); - IndexSet::ElementIterator it = locally_owned_dofs.begin(); - for (; it != locally_owned_dofs.end(); ++it) - ::dealii::internal::ElementAccess::set(work_result(*it), - *it, + for (const auto i : dof.locally_owned_dofs()) + ::dealii::internal::ElementAccess::set(work_result(i), + i, vec_result); vec_result.compress(VectorOperation::insert); } @@ -452,11 +447,11 @@ namespace VectorTools * Return whether the boundary values try to constrain a degree of freedom * that is already constrained to something else */ - template + template bool constraints_and_b_v_are_compatible( - const AffineConstraints & constraints, - std::map &boundary_values) + const AffineConstraints & constraints, + std::map &boundary_values) { for (const auto &boundary_value : boundary_values) if (constraints.is_constrained(boundary_value.first)) @@ -495,15 +490,12 @@ namespace VectorTools const Q_or_QC &q_boundary, const bool project_to_boundary_first) { - using number = typename VectorType::value_type; - Assert(dof.get_fe(0).n_components() == function.n_components, - ExcDimensionMismatch(dof.get_fe(0).n_components(), - function.n_components)); - Assert(vec_result.size() == dof.n_dofs(), - ExcDimensionMismatch(vec_result.size(), dof.n_dofs())); + using Number = typename VectorType::value_type; + AssertDimension(dof.get_fe(0).n_components(), function.n_components); + AssertDimension(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, @@ -514,11 +506,11 @@ namespace VectorTools // check if constraints are compatible (see below) const bool constraints_are_compatible = - constraints_and_b_v_are_compatible(constraints, + 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()); @@ -529,8 +521,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 constraints object does not conflict with the given boundary // values (i.e., it either does not contain boundary values or it contains @@ -539,7 +531,7 @@ namespace VectorTools // interpolate the boundary values and then condense the matrix and vector if (constraints_are_compatible) { - const Function *dummy = nullptr; + const Function *dummy = nullptr; MatrixCreator::create_mass_matrix(mapping, dof, quadrature, @@ -566,10 +558,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); @@ -584,6 +576,8 @@ namespace VectorTools vec_result); } + + template void project_parallel( @@ -596,11 +590,10 @@ namespace VectorTools const unsigned int)> & func, VectorType & vec_result) { - using Number = typename VectorType::value_type; - Assert(dof.get_fe(0).n_components() == 1, - ExcDimensionMismatch(dof.get_fe(0).n_components(), 1)); - Assert(vec_result.size() == dof.n_dofs(), - ExcDimensionMismatch(vec_result.size(), dof.n_dofs())); + using Number = typename VectorType::value_type; + using LocalVectorType = LinearAlgebra::distributed::Vector; + AssertDimension(dof.get_fe(0).n_components(), 1); + AssertDimension(vec_result.size(), dof.n_dofs()); // set up mass matrix and right hand side typename MatrixFree::AdditionalData additional_data; @@ -615,13 +608,12 @@ namespace VectorTools constraints, QGauss<1>(dof.get_fe().degree + 2), additional_data); - using MatrixType = MatrixFreeOperators:: - MassOperator>; + using MatrixType = + MatrixFreeOperators::MassOperator; MatrixType mass_matrix; mass_matrix.initialize(matrix_free); mass_matrix.compute_diagonal(); - using LocalVectorType = LinearAlgebra::distributed::Vector; LocalVectorType vec, rhs, inhomogeneities; matrix_free->initialize_dof_vector(vec); matrix_free->initialize_dof_vector(rhs); @@ -640,10 +632,7 @@ namespace VectorTools const unsigned int n_q_points = quadrature.size(); Vector cell_rhs(dofs_per_cell); std::vector local_dof_indices(dofs_per_cell); - typename DoFHandler::active_cell_iterator - cell = dof.begin_active(), - endc = dof.end(); - for (; cell != endc; ++cell) + for (const auto &cell : dof.active_cell_iterators()) if (cell->is_locally_owned()) { cell_rhs = 0; @@ -672,7 +661,7 @@ namespace VectorTools // badly conditioned matrices. This behavior can be observed, e.g. for // FE_Q_Hierarchical for degree higher than three. ReductionControl control(5 * rhs.size(), 0., 1e-12, false, false); - SolverCG> cg(control); + SolverCG cg(control); typename PreconditionJacobi::AdditionalData data(0.8); PreconditionJacobi preconditioner; preconditioner.initialize(mass_matrix, data); @@ -707,19 +696,17 @@ namespace VectorTools const DoFHandler &dof = matrix_free->get_dof_handler(fe_component); - using Number = typename VectorType::value_type; - Assert(dof.get_fe(0).n_components() == 1, - ExcDimensionMismatch(dof.get_fe(0).n_components(), 1)); - Assert(vec_result.size() == dof.n_dofs(), - ExcDimensionMismatch(vec_result.size(), dof.n_dofs())); + using Number = typename VectorType::value_type; + using LocalVectorType = LinearAlgebra::distributed::Vector; + AssertDimension(dof.get_fe(0).n_components(), 1); + AssertDimension(vec_result.size(), dof.n_dofs()); - using MatrixType = MatrixFreeOperators:: - MassOperator>; + using MatrixType = + MatrixFreeOperators::MassOperator; MatrixType mass_matrix; mass_matrix.initialize(matrix_free, {fe_component}); mass_matrix.compute_diagonal(); - using LocalVectorType = LinearAlgebra::distributed::Vector; LocalVectorType vec, rhs, inhomogeneities; matrix_free->initialize_dof_vector(vec, fe_component); matrix_free->initialize_dof_vector(rhs, fe_component); @@ -753,7 +740,7 @@ namespace VectorTools // badly conditioned matrices. This behavior can be observed, e.g. for // FE_Q_Hierarchical for degree higher than three. ReductionControl control(5 * rhs.size(), 0., 1e-12, false, false); - SolverCG> cg(control); + SolverCG cg(control); typename PreconditionJacobi::AdditionalData data(0.8); PreconditionJacobi preconditioner; preconditioner.initialize(mass_matrix, data); @@ -771,11 +758,13 @@ namespace VectorTools vec_result.compress(VectorOperation::insert); } + + /** * Specialization of project() for the case dim==spacedim and with correct - * number types for MatrixFree support. Check if we actually can use the - * MatrixFree implementation or need to use the matrix based one - * nonetheless based on the number of components. + * number types for MatrixFree support. Check if we actually can use the + * MatrixFree implementation or need to use the matrix based one nonetheless + * based on the number of components. */ template std::enable_if_t< @@ -834,6 +823,8 @@ namespace VectorTools } } + + /** * Specialization of project() for complex numbers or `dim < spacedim`, * for which we are sure that we cannot use the MatrixFree implementation. @@ -867,8 +858,11 @@ namespace VectorTools q_boundary, project_to_boundary_first); } + } // namespace internal + + template void project(const Mapping & mapping, @@ -1006,6 +1000,7 @@ namespace VectorTools } + template void project(const DoFHandler & dof, @@ -1027,6 +1022,7 @@ namespace VectorTools q_boundary, project_to_boundary_first); } + } // namespace VectorTools DEAL_II_NAMESPACE_CLOSE