From: Daniel Arndt Date: Thu, 17 Nov 2016 18:16:57 +0000 (+0100) Subject: Remove project_generic from public view X-Git-Tag: v8.5.0-rc1~372^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=ba4997d304d6511d35716846b205469057039790;p=dealii.git Remove project_generic from public view --- diff --git a/include/deal.II/numerics/vector_tools.h b/include/deal.II/numerics/vector_tools.h index f7a0996289..24c1fe87ea 100644 --- a/include/deal.II/numerics/vector_tools.h +++ b/include/deal.II/numerics/vector_tools.h @@ -933,33 +933,6 @@ namespace VectorTools const std_cxx11::function< VectorizedArray (const unsigned int, const unsigned int)> func, VectorType &vec_result); - /** - * Implementation for the project() function with finite elements - * are supported by the MatrixFree class for arbitrary number of - * components and degree of the FiniteElement. - * - * This function should be used if you have more than four components - * or the degree of your FiniteElement is higher than eight. For all the - * other cases project() is already instantiated. - * - * The first two template arguments have to be specified explicitly. - * @p vec_result is expected to not have any ghost entries. - * @p project_to_boundary_first and @p enforce_zero_boundary are not yet - * implemented. - */ - template - void project_generic (const Mapping &mapping, - const DoFHandler &dof, - const ConstraintMatrix &constraints, - const Quadrature &quadrature, - const Function &function, - VectorType &vec_result, - const bool enforce_zero_boundary = false, - const Quadrature &q_boundary = (dim > 1 ? - QGauss(2) : - Quadrature(0)), - const bool project_to_boundary_first = false); - /** * Compute Dirichlet boundary conditions. This function makes up a map of * degrees of freedom subject to Dirichlet boundary conditions and the diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h index 23ed9d56a4..03266f2ac2 100644 --- a/include/deal.II/numerics/vector_tools.templates.h +++ b/include/deal.II/numerics/vector_tools.templates.h @@ -895,14 +895,98 @@ namespace VectorTools - /** - * Helper interface. After figuring out the number of components - * in project_parallel, we determine the degree of the FiniteElement - * and call project_generic with the appropriate template arguments. + /* + * Generic implementation for project using a MatrixFree implementation + * with arbitrary number of components and degree of the FiniteElement. */ + template + void project_generic (const Mapping &mapping, + const DoFHandler &dof, + const ConstraintMatrix &constraints, + const Quadrature &quadrature, + const Function &function, + VectorType &vec_result, + const bool enforce_zero_boundary, + const Quadrature &q_boundary, + const bool project_to_boundary_first) + { + Assert (project_to_boundary_first == false, ExcNotImplemented()); + Assert (enforce_zero_boundary == false, ExcNotImplemented()); + (void) enforce_zero_boundary; + (void) project_to_boundary_first; + (void) q_boundary; + + const IndexSet locally_owned_dofs = dof.locally_owned_dofs(); + IndexSet locally_relevant_dofs; + DoFTools::extract_locally_relevant_dofs(dof, locally_relevant_dofs); + + typedef typename VectorType::value_type number; + Assert (dof.get_fe().n_components() == function.n_components, + ExcDimensionMismatch(dof.get_fe().n_components(), + function.n_components)); + Assert (vec_result.size() == dof.n_dofs(), + ExcDimensionMismatch (vec_result.size(), dof.n_dofs())); + Assert (dof.get_fe().degree == fe_degree, + ExcDimensionMismatch(fe_degree, dof.get_fe().degree)); + Assert (dof.get_fe().n_components() == components, + ExcDimensionMismatch(components, dof.get_fe().n_components())); + + // set up mass matrix and right hand side + typename MatrixFree::AdditionalData additional_data; + additional_data.tasks_parallel_scheme = + MatrixFree::AdditionalData::partition_color; + if (const parallel::Triangulation *parallel_tria = + dynamic_cast*>(&dof.get_tria())) + additional_data.mpi_communicator = parallel_tria->get_communicator(); + additional_data.mapping_update_flags = (update_values | update_JxW_values); + MatrixFree matrix_free; + matrix_free.reinit (mapping, dof, constraints, + QGauss<1>(fe_degree+2), additional_data); + typedef MatrixFreeOperators::MassOperator MatrixType; + MatrixType mass_matrix; + mass_matrix.initialize(matrix_free); + mass_matrix.compute_diagonal(); + + typedef LinearAlgebra::distributed::Vector LocalVectorType; + LocalVectorType vec, rhs, inhomogeneities; + matrix_free.initialize_dof_vector(vec); + matrix_free.initialize_dof_vector(rhs); + matrix_free.initialize_dof_vector(inhomogeneities); + constraints.distribute(inhomogeneities); + inhomogeneities*=-1.; + + create_right_hand_side (mapping, dof, quadrature, function, rhs, constraints); + mass_matrix.vmult_add(rhs, inhomogeneities); + + // now invert the matrix + // Allow for a maximum of 5*n steps to reduce the residual by 10^-12. n + // steps may not be sufficient, since roundoff errors may accumulate for + // 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); + PreconditionJacobi preconditioner; + preconditioner.initialize(mass_matrix, 1.); + cg.solve (mass_matrix, vec, rhs, preconditioner); + vec+=inhomogeneities; + + constraints.distribute (vec); + IndexSet::ElementIterator it = locally_owned_dofs.begin(); + for (; it!=locally_owned_dofs.end(); ++it) + vec_result(*it) = vec(*it); + vec_result.compress(VectorOperation::insert); + } + + + + /** + * Helper interface. After figuring out the number of components + * in project_parallel, we determine the degree of the FiniteElement + * and call project_generic with the appropriate template arguments. + */ template void project_parallel_helper (const Mapping &mapping, - const dealii::DoFHandler &dof, + const DoFHandler &dof, const ConstraintMatrix &constraints, const Quadrature &quadrature, const Function &function, @@ -962,9 +1046,7 @@ namespace VectorTools break; default: - Assert(false, - ExcMessage("Call project_generic yourself specifying the " - "appropriate template arguments!")); + Assert(false, ExcNotImplemented()); } } @@ -974,7 +1056,7 @@ namespace VectorTools // Used to determine the number of components. template void project_parallel (const Mapping &mapping, - const dealii::DoFHandler &dof, + const DoFHandler &dof, const ConstraintMatrix &constraints, const Quadrature &quadrature, const Function &function, @@ -1010,9 +1092,7 @@ namespace VectorTools break; default: - Assert(false, - ExcMessage("Call project_generic yourself specifying the " - "appropriate template arguments!")); + Assert(false, ExcNotImplemented()); } } @@ -1020,7 +1100,7 @@ namespace VectorTools template void project_matrix_free (const Mapping &mapping, - const dealii::DoFHandler &dof, + const DoFHandler &dof, const ConstraintMatrix &constraints, const Quadrature &quadrature, const Function &function, @@ -1084,9 +1164,6 @@ namespace VectorTools const std_cxx11::function< typename VectorType::value_type (const typename DoFHandler::active_cell_iterator &, const unsigned int)> func, VectorType &vec_result) { - const parallel::Triangulation *parallel_tria = - dynamic_cast*>(&dof.get_tria()); - typedef typename VectorType::value_type Number; Assert (dof.get_fe().n_components() == 1, ExcDimensionMismatch(dof.get_fe().n_components(), @@ -1100,7 +1177,8 @@ namespace VectorTools typename MatrixFree::AdditionalData additional_data; additional_data.tasks_parallel_scheme = MatrixFree::AdditionalData::partition_color; - if (parallel_tria !=0) + if (const parallel::Triangulation *parallel_tria = + dynamic_cast*>(&dof.get_tria())) additional_data.mpi_communicator = parallel_tria->get_communicator(); additional_data.mapping_update_flags = (update_values | update_JxW_values); MatrixFree matrix_free; @@ -1356,87 +1434,6 @@ namespace VectorTools - template - void project_generic (const Mapping &mapping, - const DoFHandler &dof, - const ConstraintMatrix &constraints, - const Quadrature &quadrature, - const Function &function, - VectorType &vec_result, - const bool enforce_zero_boundary, - const Quadrature &q_boundary, - const bool project_to_boundary_first) - { - const parallel::Triangulation *parallel_tria = - dynamic_cast*>(&dof.get_tria()); - Assert (project_to_boundary_first == false, ExcNotImplemented()); - Assert (enforce_zero_boundary == false, ExcNotImplemented()); - (void) enforce_zero_boundary; - (void) project_to_boundary_first; - (void) q_boundary; - - const IndexSet locally_owned_dofs = dof.locally_owned_dofs(); - IndexSet locally_relevant_dofs; - DoFTools::extract_locally_relevant_dofs(dof, locally_relevant_dofs); - - typedef typename VectorType::value_type number; - Assert (dof.get_fe().n_components() == function.n_components, - ExcDimensionMismatch(dof.get_fe().n_components(), - function.n_components)); - Assert (vec_result.size() == dof.n_dofs(), - ExcDimensionMismatch (vec_result.size(), dof.n_dofs())); - Assert (dof.get_fe().degree == fe_degree, - ExcDimensionMismatch(fe_degree, dof.get_fe().degree)); - Assert (dof.get_fe().n_components() == components, - ExcDimensionMismatch(components, dof.get_fe().n_components())); - - // set up mass matrix and right hand side - typename MatrixFree::AdditionalData additional_data; - additional_data.tasks_parallel_scheme = - MatrixFree::AdditionalData::partition_color; - if (parallel_tria) - additional_data.mpi_communicator = parallel_tria->get_communicator(); - additional_data.mapping_update_flags = (update_values | update_JxW_values); - MatrixFree matrix_free; - matrix_free.reinit (mapping, dof, constraints, - QGauss<1>(fe_degree+2), additional_data); - typedef MatrixFreeOperators::MassOperator MatrixType; - MatrixType mass_matrix; - mass_matrix.initialize(matrix_free); - mass_matrix.compute_diagonal(); - - typedef LinearAlgebra::distributed::Vector LocalVectorType; - LocalVectorType vec, rhs, inhomogeneities; - matrix_free.initialize_dof_vector(vec); - matrix_free.initialize_dof_vector(rhs); - matrix_free.initialize_dof_vector(inhomogeneities); - constraints.distribute(inhomogeneities); - inhomogeneities*=-1.; - - create_right_hand_side (mapping, dof, quadrature, function, rhs, constraints); - mass_matrix.vmult_add(rhs, inhomogeneities); - - // now invert the matrix - // Allow for a maximum of 5*n steps to reduce the residual by 10^-12. n - // steps may not be sufficient, since roundoff errors may accumulate for - // 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); - PreconditionJacobi preconditioner; - preconditioner.initialize(mass_matrix, 1.); - cg.solve (mass_matrix, vec, rhs, preconditioner); - vec+=inhomogeneities; - - constraints.distribute (vec); - IndexSet::ElementIterator it = locally_owned_dofs.begin(); - for (; it!=locally_owned_dofs.end(); ++it) - vec_result(*it) = vec(*it); - vec_result.compress(VectorOperation::insert); - } - - - template void project (const Mapping &mapping, const DoFHandler &dof,