From: Matthias Maier Date: Sat, 26 May 2018 01:21:29 +0000 (-0500) Subject: RCM: Update numerics, part 3: vector_tools X-Git-Tag: v9.1.0-rc1~1067^2~19 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=90549e7b2caa7bde953d4944aae4d4348c764394;p=dealii.git RCM: Update numerics, part 3: vector_tools --- diff --git a/include/deal.II/numerics/vector_tools.h b/include/deal.II/numerics/vector_tools.h index b6ad27195b..8cc2624ff6 100644 --- a/include/deal.II/numerics/vector_tools.h +++ b/include/deal.II/numerics/vector_tools.h @@ -58,7 +58,8 @@ namespace hp template class QCollection; } -class ConstraintMatrix; +template +class AffineConstraints; // TODO: Move documentation of functions to the functions! @@ -731,8 +732,8 @@ namespace VectorTools * * For this case (continuous elements on grids with hanging nodes), please * use the interpolate_to_different_mesh function with an additional - * ConstraintMatrix argument, see below, or make the field conforming - * yourself by calling the @p ConstraintsMatrix::distribute function of your + * AffineConstraints argument, see below, or make the field conforming + * yourself by calling the @p AffineConstraints::distribute function of your * hanging node constraints object. * * @note: This function works with parallel::distributed::Triangulation, but @@ -768,11 +769,12 @@ namespace VectorTools typename VectorType, template class DoFHandlerType> void - interpolate_to_different_mesh(const DoFHandlerType &dof1, - const VectorType & u1, - const DoFHandlerType &dof2, - const ConstraintMatrix &constraints, - VectorType & u2); + interpolate_to_different_mesh( + const DoFHandlerType & dof1, + const VectorType & u1, + const DoFHandlerType & dof2, + const AffineConstraints &constraints, + VectorType & u2); /** * The same function as above, but takes an InterGridMap object directly as @@ -787,10 +789,10 @@ namespace VectorTools template class DoFHandlerType> void interpolate_to_different_mesh( - const InterGridMap> &intergridmap, - const VectorType & u1, - const ConstraintMatrix & constraints, - VectorType & u2); + const InterGridMap> & intergridmap, + const VectorType & u1, + const AffineConstraints &constraints, + VectorType & u2); /** * Compute the projection of @p function to the finite element space. @@ -845,10 +847,10 @@ namespace VectorTools */ template void - project(const Mapping & mapping, - const DoFHandler &dof, - const ConstraintMatrix & constraints, - const Quadrature & quadrature, + project(const Mapping & mapping, + const DoFHandler & dof, + const AffineConstraints &constraints, + const Quadrature & quadrature, const Function &function, VectorType & vec, const bool enforce_zero_boundary = false, @@ -863,9 +865,9 @@ namespace VectorTools */ template void - project(const DoFHandler &dof, - const ConstraintMatrix & constraints, - const Quadrature & quadrature, + project(const DoFHandler & dof, + const AffineConstraints &constraints, + const Quadrature & quadrature, const Function &function, VectorType & vec, const bool enforce_zero_boundary = false, @@ -880,10 +882,10 @@ namespace VectorTools */ template void - project(const hp::MappingCollection &mapping, - const hp::DoFHandler & dof, - const ConstraintMatrix & constraints, - const hp::QCollection & quadrature, + project(const hp::MappingCollection & mapping, + const hp::DoFHandler & dof, + const AffineConstraints &constraints, + const hp::QCollection & quadrature, const Function &function, VectorType & vec, const bool enforce_zero_boundary = false, @@ -897,9 +899,9 @@ namespace VectorTools */ template void - project(const hp::DoFHandler &dof, - const ConstraintMatrix & constraints, - const hp::QCollection & quadrature, + project(const hp::DoFHandler & dof, + const AffineConstraints &constraints, + const hp::QCollection & quadrature, const Function &function, VectorType & vec, const bool enforce_zero_boundary = false, @@ -930,14 +932,14 @@ namespace VectorTools */ template void - project(const Mapping & mapping, - const DoFHandler &dof, - const ConstraintMatrix & constraints, - const Quadrature & quadrature, + project(const Mapping & mapping, + const DoFHandler & dof, + const AffineConstraints &constraints, + const Quadrature & quadrature, const std::function::active_cell_iterator &, - const unsigned int)> & func, - VectorType & vec_result); + const unsigned int)> & func, + VectorType & vec_result); /** * The same as above for projection of scalar-valued MatrixFree quadrature @@ -967,7 +969,7 @@ namespace VectorTools void project(std::shared_ptr< const MatrixFree> data, - const ConstraintMatrix & constraints, + const AffineConstraints &constraints, const unsigned int n_q_points_1d, const std::function( const unsigned int, @@ -983,7 +985,7 @@ namespace VectorTools void project(std::shared_ptr< const MatrixFree> data, - const ConstraintMatrix & constraints, + const AffineConstraints &constraints, const std::function( const unsigned int, const unsigned int)> & func, @@ -1145,7 +1147,7 @@ namespace VectorTools * * @note When combining adaptively refined meshes with hanging node * constraints and boundary conditions like from the current function within - * one ConstraintMatrix object, the hanging node constraints should always + * one AffineConstraints object, the hanging node constraints should always * be set first, and then the boundary conditions since boundary conditions * are not set in the second operation on degrees of freedom that are * already constrained. This makes sure that the discretization remains @@ -1201,9 +1203,9 @@ namespace VectorTools const Mapping & mapping, const DoFHandlerType &dof, const std::map *> - & function_map, - ConstraintMatrix & constraints, - const ComponentMask &component_mask = ComponentMask()); + & function_map, + AffineConstraints &constraints, + const ComponentMask & component_mask = ComponentMask()); /** * Same function as above, but taking only one pair of boundary indicator @@ -1226,7 +1228,7 @@ namespace VectorTools const DoFHandlerType &dof, const types::boundary_id boundary_component, const Function & boundary_function, - ConstraintMatrix & constraints, + AffineConstraints & constraints, const ComponentMask & component_mask = ComponentMask()); /** @@ -1249,7 +1251,7 @@ namespace VectorTools const DoFHandlerType &dof, const types::boundary_id boundary_component, const Function & boundary_function, - ConstraintMatrix & constraints, + AffineConstraints & constraints, const ComponentMask & component_mask = ComponentMask()); @@ -1269,9 +1271,9 @@ namespace VectorTools interpolate_boundary_values( const DoFHandlerType &dof, const std::map *> - & function_map, - ConstraintMatrix & constraints, - const ComponentMask &component_mask = ComponentMask()); + & function_map, + AffineConstraints &constraints, + const ComponentMask & component_mask = ComponentMask()); /** @@ -1405,7 +1407,7 @@ namespace VectorTools * * @note When combining adaptively refined meshes with hanging node * constraints and boundary conditions like from the current function within - * one ConstraintMatrix object, the hanging node constraints should always + * one AffineConstraints object, the hanging node constraints should always * be set first, and then the boundary conditions since boundary conditions * are not set in the second operation on degrees of freedom that are * already constrained. This makes sure that the discretization remains @@ -1437,7 +1439,7 @@ namespace VectorTools const std::map *> & boundary_functions, const Quadrature &q, - ConstraintMatrix & constraints, + AffineConstraints &constraints, std::vector component_mapping = std::vector()); /** @@ -1453,7 +1455,7 @@ namespace VectorTools const std::map *> & boundary_function, const Quadrature &q, - ConstraintMatrix & constraints, + AffineConstraints &constraints, std::vector component_mapping = std::vector()); @@ -1518,7 +1520,7 @@ namespace VectorTools const unsigned int first_vector_component, const Function &boundary_function, const types::boundary_id boundary_component, - ConstraintMatrix & constraints, + AffineConstraints & constraints, const Mapping & mapping = StaticMappingQ1::mapping); /** @@ -1536,7 +1538,7 @@ namespace VectorTools const unsigned int first_vector_component, const Function & boundary_function, const types::boundary_id boundary_component, - ConstraintMatrix & constraints, + AffineConstraints & constraints, const hp::MappingCollection &mapping_collection = hp::StaticMappingQ1::mapping_collection); @@ -1642,7 +1644,7 @@ namespace VectorTools const unsigned int first_vector_component, const Function &boundary_function, const types::boundary_id boundary_component, - ConstraintMatrix & constraints, + AffineConstraints & constraints, const Mapping & mapping = StaticMappingQ1::mapping); @@ -1659,7 +1661,7 @@ namespace VectorTools const unsigned int first_vector_component, const Function & boundary_function, const types::boundary_id boundary_component, - ConstraintMatrix & constraints, + AffineConstraints & constraints, const hp::MappingCollection &mapping_collection = hp::StaticMappingQ1::mapping_collection); @@ -1675,7 +1677,7 @@ namespace VectorTools * elements. Thus it throws an exception, if it is called with other finite * elements. * - * If the ConstraintMatrix @p constraints contained values or other + * If the AffineConstraints object @p constraints contained values or other * constraints before, the new ones are added or the old ones overwritten, * if a node of the boundary part to be used was already in the list of * constraints. This is handled by using inhomogeneous constraints. Please @@ -1717,7 +1719,7 @@ namespace VectorTools const unsigned int first_vector_component, const Function &boundary_function, const types::boundary_id boundary_component, - ConstraintMatrix & constraints, + AffineConstraints & constraints, const Mapping & mapping = StaticMappingQ1::mapping); /** @@ -1735,7 +1737,7 @@ namespace VectorTools const unsigned int first_vector_component, const Function & boundary_function, const types::boundary_id boundary_component, - ConstraintMatrix & constraints, + AffineConstraints & constraints, const hp::MappingCollection &mapping_collection = hp::StaticMappingQ1::mapping_collection); @@ -1973,7 +1975,7 @@ namespace VectorTools const unsigned int first_vector_component, const std::set & boundary_ids, typename FunctionMap::type &function_map, - ConstraintMatrix & constraints, + AffineConstraints & constraints, const Mapping &mapping = StaticMappingQ1::mapping); /** @@ -1994,7 +1996,7 @@ namespace VectorTools const DoFHandlerType &dof_handler, const unsigned int first_vector_component, const std::set & boundary_ids, - ConstraintMatrix & constraints, + AffineConstraints & constraints, const Mapping &mapping = StaticMappingQ1::mapping); /** @@ -2020,7 +2022,7 @@ namespace VectorTools const unsigned int first_vector_component, const std::set & boundary_ids, typename FunctionMap::type &function_map, - ConstraintMatrix & constraints, + AffineConstraints & constraints, const Mapping &mapping = StaticMappingQ1::mapping); /** @@ -2037,7 +2039,7 @@ namespace VectorTools const DoFHandlerType &dof_handler, const unsigned int first_vector_component, const std::set & boundary_ids, - ConstraintMatrix & constraints, + AffineConstraints & constraints, const Mapping &mapping = StaticMappingQ1::mapping); @@ -2061,7 +2063,8 @@ namespace VectorTools const Quadrature & q, const Function &rhs, VectorType & rhs_vector, - const ConstraintMatrix &constraints = ConstraintMatrix()); + const AffineConstraints & constraints = + AffineConstraints()); /** * Call the create_right_hand_side() function, see above, with @@ -2074,7 +2077,8 @@ namespace VectorTools const Quadrature & q, const Function &rhs, VectorType & rhs_vector, - const ConstraintMatrix &constraints = ConstraintMatrix()); + const AffineConstraints & constraints = + AffineConstraints()); /** * Like the previous set of functions, but for hp objects. @@ -2087,7 +2091,8 @@ namespace VectorTools const hp::QCollection & q, const Function &rhs, VectorType & rhs_vector, - const ConstraintMatrix &constraints = ConstraintMatrix()); + const AffineConstraints & constraints = + AffineConstraints()); /** * Like the previous set of functions, but for hp objects. @@ -2099,7 +2104,8 @@ namespace VectorTools const hp::QCollection & q, const Function &rhs, VectorType & rhs_vector, - const ConstraintMatrix &constraints = ConstraintMatrix()); + const AffineConstraints & constraints = + AffineConstraints()); /** * Create a right hand side vector for a point source at point @p p. In diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h index 67b1d22127..3260c64c43 100644 --- a/include/deal.II/numerics/vector_tools.templates.h +++ b/include/deal.II/numerics/vector_tools.templates.h @@ -55,8 +55,8 @@ #include #include +#include #include -#include #include #include #include @@ -747,7 +747,7 @@ namespace VectorTools InterGridMap> intergridmap; intergridmap.make_mapping(dof1, dof2); - ConstraintMatrix dummy; + AffineConstraints dummy; dummy.close(); interpolate_to_different_mesh(intergridmap, u1, dummy, u2); @@ -760,11 +760,12 @@ namespace VectorTools typename VectorType, template class DoFHandlerType> void - interpolate_to_different_mesh(const DoFHandlerType &dof1, - const VectorType & u1, - const DoFHandlerType &dof2, - const ConstraintMatrix &constraints, - VectorType & u2) + interpolate_to_different_mesh( + const DoFHandlerType & dof1, + const VectorType & u1, + const DoFHandlerType & dof2, + const AffineConstraints &constraints, + VectorType & u2) { Assert(GridTools::have_same_coarse_mesh(dof1, dof2), ExcMessage("The two DoF handlers must represent triangulations that " @@ -803,10 +804,10 @@ namespace VectorTools template class DoFHandlerType> void interpolate_to_different_mesh( - const InterGridMap> &intergridmap, - const VectorType & u1, - const ConstraintMatrix & constraints, - VectorType & u2) + const InterGridMap> & intergridmap, + const VectorType & u1, + const AffineConstraints &constraints, + VectorType & u2) { const DoFHandlerType &dof1 = intergridmap.get_source_grid(); const DoFHandlerType &dof2 = @@ -938,7 +939,7 @@ namespace VectorTools template bool constraints_and_b_v_are_compatible( - const ConstraintMatrix & constraints, + const AffineConstraints & constraints, std::map &boundary_values) { for (typename std::map::iterator it = @@ -1000,7 +1001,7 @@ namespace VectorTools do_project( const M_or_MC & mapping, const DoFHandlerType & dof, - const ConstraintMatrix & constraints, + const AffineConstraints & constraints, const Q_or_QC & quadrature, const Function &function, VectorType & vec_result, @@ -1099,7 +1100,7 @@ namespace VectorTools project_matrix_free( const Mapping & mapping, const DoFHandler &dof, - const ConstraintMatrix & constraints, + const AffineConstraints &constraints, const Quadrature & quadrature, const Function< spacedim, @@ -1206,7 +1207,7 @@ namespace VectorTools project_matrix_free_degree( const Mapping & mapping, const DoFHandler &dof, - const ConstraintMatrix & constraints, + const AffineConstraints &constraints, const Quadrature & quadrature, const Function< spacedim, @@ -1277,7 +1278,7 @@ namespace VectorTools project_matrix_free_component( const Mapping & mapping, const DoFHandler &dof, - const ConstraintMatrix & constraints, + const AffineConstraints &constraints, const Quadrature & quadrature, const Function< spacedim, @@ -1355,7 +1356,7 @@ namespace VectorTools project_matrix_free_copy_vector( const Mapping & mapping, const DoFHandler & dof, - const ConstraintMatrix & constraints, + const AffineConstraints & constraints, const Quadrature & quadrature, const Function &function, VectorType & vec_result, @@ -1395,15 +1396,16 @@ namespace VectorTools */ template void - project(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) + project( + const Mapping & mapping, + const DoFHandler & dof, + const AffineConstraints &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) { // If we can, use the matrix-free implementation bool use_matrix_free = @@ -1450,14 +1452,14 @@ namespace VectorTools template void project_parallel( - const Mapping & mapping, - const DoFHandler &dof, - const ConstraintMatrix & constraints, - const Quadrature & quadrature, + const Mapping & mapping, + const DoFHandler & dof, + const AffineConstraints &constraints, + const Quadrature & quadrature, const std::function::active_cell_iterator &, - const unsigned int)> & func, - VectorType & vec_result) + const unsigned int)> & func, + VectorType & vec_result) { typedef typename VectorType::value_type Number; Assert(dof.get_fe(0).n_components() == 1, @@ -1567,13 +1569,13 @@ namespace VectorTools void project_parallel( std::shared_ptr> - matrix_free, - const ConstraintMatrix &constraints, + matrix_free, + const AffineConstraints &constraints, const std::function( const unsigned int, - const unsigned int)> &func, - VectorType & vec_result, - const unsigned int fe_component) + const unsigned int)> & func, + VectorType & vec_result, + const unsigned int fe_component) { const DoFHandler &dof = matrix_free->get_dof_handler(fe_component); @@ -1655,14 +1657,14 @@ namespace VectorTools template void - project(const Mapping & mapping, - const DoFHandler &dof, - const ConstraintMatrix & constraints, - const Quadrature & quadrature, + project(const Mapping & mapping, + const DoFHandler & dof, + const AffineConstraints &constraints, + const Quadrature & quadrature, const std::function::active_cell_iterator &, - const unsigned int)> & func, - VectorType & vec_result) + const unsigned int)> & func, + VectorType & vec_result) { switch (dof.get_fe().degree) { @@ -1690,7 +1692,7 @@ namespace VectorTools void project(std::shared_ptr< const MatrixFree> matrix_free, - const ConstraintMatrix & constraints, + const AffineConstraints &constraints, const unsigned int n_q_points_1d, const std::function( const unsigned int, @@ -1731,7 +1733,7 @@ namespace VectorTools void project(std::shared_ptr< const MatrixFree> matrix_free, - const ConstraintMatrix & constraints, + const AffineConstraints &constraints, const std::function( const unsigned int, const unsigned int)> & func, @@ -1750,10 +1752,10 @@ namespace VectorTools template void - project(const Mapping & mapping, - const DoFHandler &dof, - const ConstraintMatrix & constraints, - const Quadrature & quadrature, + project(const Mapping & mapping, + const DoFHandler & dof, + const AffineConstraints &constraints, + const Quadrature & quadrature, const Function &function, VectorType & vec_result, const bool enforce_zero_boundary, @@ -1803,9 +1805,9 @@ namespace VectorTools template void - project(const DoFHandler &dof, - const ConstraintMatrix & constraints, - const Quadrature & quadrature, + project(const DoFHandler & dof, + const AffineConstraints &constraints, + const Quadrature & quadrature, const Function &function, VectorType & vec, const bool enforce_zero_boundary, @@ -1833,10 +1835,10 @@ namespace VectorTools template void - project(const hp::MappingCollection &mapping, - const hp::DoFHandler & dof, - const ConstraintMatrix & constraints, - const hp::QCollection & quadrature, + project(const hp::MappingCollection & mapping, + const hp::DoFHandler & dof, + const AffineConstraints &constraints, + const hp::QCollection & quadrature, const Function &function, VectorType & vec_result, const bool enforce_zero_boundary, @@ -1861,9 +1863,9 @@ namespace VectorTools template void - project(const hp::DoFHandler &dof, - const ConstraintMatrix & constraints, - const hp::QCollection & quadrature, + project(const hp::DoFHandler & dof, + const AffineConstraints &constraints, + const hp::QCollection & quadrature, const Function &function, VectorType & vec, const bool enforce_zero_boundary, @@ -1891,7 +1893,7 @@ namespace VectorTools const Quadrature & quadrature, const Function &rhs_function, VectorType & rhs_vector, - const ConstraintMatrix & constraints) + const AffineConstraints & constraints) { typedef typename VectorType::value_type Number; @@ -2007,7 +2009,7 @@ namespace VectorTools const Quadrature & quadrature, const Function &rhs_function, VectorType & rhs_vector, - const ConstraintMatrix & constraints) + const AffineConstraints & constraints) { create_right_hand_side(StaticMappingQ1::mapping, dof_handler, @@ -2027,7 +2029,7 @@ namespace VectorTools const hp::QCollection & quadrature, const Function &rhs_function, VectorType & rhs_vector, - const ConstraintMatrix & constraints) + const AffineConstraints & constraints) { typedef typename VectorType::value_type Number; @@ -2160,7 +2162,7 @@ namespace VectorTools const hp::QCollection & quadrature, const Function &rhs_function, VectorType & rhs_vector, - const ConstraintMatrix & constraints) + const AffineConstraints & constraints) { create_right_hand_side( hp::StaticMappingQ1::mapping_collection, @@ -3130,7 +3132,8 @@ namespace VectorTools - // ----------- interpolate_boundary_values for ConstraintMatrix -------------- + // ----------- interpolate_boundary_values for AffineConstraints + // -------------- @@ -3143,9 +3146,9 @@ namespace VectorTools const Mapping & mapping, const DoFHandlerType &dof, const std::map *> - & function_map, - ConstraintMatrix & constraints, - const ComponentMask &component_mask_) + & function_map, + AffineConstraints &constraints, + const ComponentMask & component_mask_) { std::map boundary_values; interpolate_boundary_values( @@ -3176,7 +3179,7 @@ namespace VectorTools const DoFHandlerType &dof, const types::boundary_id boundary_component, const Function & boundary_function, - ConstraintMatrix & constraints, + AffineConstraints & constraints, const ComponentMask & component_mask) { std::map *> @@ -3197,7 +3200,7 @@ namespace VectorTools const DoFHandlerType &dof, const types::boundary_id boundary_component, const Function & boundary_function, - ConstraintMatrix & constraints, + AffineConstraints & constraints, const ComponentMask & component_mask) { interpolate_boundary_values(StaticMappingQ1::mapping, @@ -3218,9 +3221,9 @@ namespace VectorTools interpolate_boundary_values( const DoFHandlerType &dof, const std::map *> - & function_map, - ConstraintMatrix & constraints, - const ComponentMask &component_mask) + & function_map, + AffineConstraints &constraints, + const ComponentMask & component_mask) { interpolate_boundary_values(StaticMappingQ1::mapping, dof, @@ -3479,9 +3482,8 @@ namespace VectorTools // make mass matrix and right hand side - SparseMatrix::real_type> - mass_matrix(sparsity); - Vector rhs(sparsity.n_rows()); + SparseMatrix mass_matrix(sparsity); + Vector rhs(sparsity.n_rows()); MatrixCreator::create_boundary_mass_matrix( @@ -3492,9 +3494,7 @@ namespace VectorTools boundary_functions, rhs, dof_to_boundary_mapping, - (const Function::real_type> - *)nullptr, + (const Function *)nullptr, component_mapping); // For certain weird elements, @@ -3509,11 +3509,9 @@ namespace VectorTools // TODO: Maybe we should figure out if the element really needs this - FilteredMatrix::real_type>> - filtered_mass_matrix(mass_matrix, true); - FilteredMatrix::real_type>> - filtered_precondition; - std::vector excluded_dofs(mass_matrix.m(), false); + FilteredMatrix> filtered_mass_matrix(mass_matrix, true); + FilteredMatrix> filtered_precondition; + std::vector excluded_dofs(mass_matrix.m(), false); // we assemble mass matrix with unit weight, // thus it will be real-valued irrespectively of the underlying algebra @@ -3643,7 +3641,7 @@ namespace VectorTools } - // ---- implementation for project_boundary_values with ConstraintMatrix ---- + // ---- implementation for project_boundary_values with AffineConstraints ---- @@ -3655,7 +3653,7 @@ namespace VectorTools const std::map *> & boundary_functions, const Quadrature &q, - ConstraintMatrix & constraints, + AffineConstraints &constraints, std::vector component_mapping) { std::map boundary_values; @@ -3683,7 +3681,7 @@ namespace VectorTools const std::map *> & boundary_functions, const Quadrature &q, - ConstraintMatrix & constraints, + AffineConstraints &constraints, std::vector component_mapping) { project_boundary_values(StaticMappingQ1::mapping, @@ -3768,69 +3766,42 @@ namespace VectorTools void add_constraint(const VectorDoFTuple &dof_indices, const Tensor<1, dim> & constraining_vector, - ConstraintMatrix & constraints, + AffineConstraints &constraints, const double inhomogeneity = 0) { - // choose the DoF that has the - // largest component in the - // constraining_vector as the - // one to be constrained as - // this makes the process - // stable in cases where the - // constraining_vector has the - // form n=(1,0) or n=(0,1) + // choose the DoF that has the largest component in the + // constraining_vector as the one to be constrained as this makes the + // process stable in cases where the constraining_vector has the form + // n=(1,0) or n=(0,1) // - // we get constraints of the form - // x0 = a_1*x1 + a2*x2 + ... - // if one of the weights is - // essentially zero then skip - // this part. the ConstraintMatrix - // can also deal with cases like - // x0 = 0 - // if necessary + // we get constraints of the form x0 = a_1*x1 + a2*x2 + ... if one of + // the weights is essentially zero then skip this part. the + // AffineConstraints can also deal with cases like x0 = 0 if + // necessary // - // there is a problem if we have a - // normal vector of the form - // (a,a,small) or (a,a,a). Depending on - // round-off we may choose the first or - // second component (or third, in the - // latter case) as the largest one, and - // depending on our choice one or - // another degree of freedom will be - // constrained. On a single processor - // this is not much of a problem, but - // it's a nightmare when we run in - // parallel and two processors disagree - // on which DoF should be - // constrained. This led to an - // incredibly difficult to find bug in - // step-32 when running in parallel - // with 9 or more processors. + // there is a problem if we have a normal vector of the form + // (a,a,small) or (a,a,a). Depending on round-off we may choose the + // first or second component (or third, in the latter case) as the + // largest one, and depending on our choice one or another degree of + // freedom will be constrained. On a single processor this is not + // much of a problem, but it's a nightmare when we run in parallel + // and two processors disagree on which DoF should be constrained. + // This led to an incredibly difficult to find bug in step-32 when + // running in parallel with 9 or more processors. // - // in practice, such normal vectors of - // the form (a,a,small) or (a,a,a) - // happen not infrequently since they - // lie on the diagonals where vertices - // frequently happen to land upon mesh - // refinement if one starts from a - // symmetric and regular body. we work - // around this problem in the following - // way: if we have a normal vector of - // the form (a,b) (similarly algorithm - // in 3d), we choose 'a' as the largest - // coefficient not if a>b but if - // a>b+1e-10. this shifts the problem - // away from the frequently visited - // diagonal to a line that is off the - // diagonal. there will of course be - // problems where the exact values of a - // and b differ by exactly 1e-10 and we - // get into the same instability, but - // from a practical viewpoint such - // problems should be much rarer. in - // particular, meshes have to be very - // fine for a vertex to land on - // this line if the original body had a + // in practice, such normal vectors of the form (a,a,small) or + // (a,a,a) happen not infrequently since they lie on the diagonals + // where vertices frequently happen to land upon mesh refinement if + // one starts from a symmetric and regular body. we work around this + // problem in the following way: if we have a normal vector of the + // form (a,b) (similarly algorithm in 3d), we choose 'a' as the + // largest coefficient not if a>b but if a>b+1e-10. this shifts the + // problem away from the frequently visited diagonal to a line that + // is off the diagonal. there will of course be problems where the + // exact values of a and b differ by exactly 1e-10 and we get into + // the same instability, but from a practical viewpoint such problems + // should be much rarer. in particular, meshes have to be very fine + // for a vertex to land on this line if the original body had a // vertex on the diagonal as well switch (dim) { @@ -4008,7 +3979,7 @@ namespace VectorTools add_tangentiality_constraints( const VectorDoFTuple &dof_indices, const Tensor<1, dim> & tangent_vector, - ConstraintMatrix & constraints, + AffineConstraints &constraints, const Vector & b_values = Vector(dim)) { // choose the DoF that has the @@ -4808,12 +4779,12 @@ namespace VectorTools void project_boundary_values_curl_conforming( - const DoFHandler & dof_handler, - const unsigned int first_vector_component, - const Function & boundary_function, - const types::boundary_id boundary_component, - ConstraintMatrix & constraints, - const Mapping & mapping) + const DoFHandler & dof_handler, + const unsigned int first_vector_component, + const Function & boundary_function, + const types::boundary_id boundary_component, + AffineConstraints &constraints, + const Mapping & mapping) { // Projection-based interpolation is performed in two (in 2D) respectively // three (in 3D) steps. First the tangential component of the function is @@ -5061,7 +5032,7 @@ namespace VectorTools const unsigned int first_vector_component, const Function & boundary_function, const types::boundary_id boundary_component, - ConstraintMatrix & constraints, + AffineConstraints & constraints, const hp::MappingCollection &mapping_collection) { const hp::FECollection &fe_collection(dof_handler.get_fe_collection()); @@ -5986,7 +5957,7 @@ namespace VectorTools const unsigned int first_vector_component, const Function & boundary_function, const types::boundary_id boundary_component, - ConstraintMatrix & constraints, + AffineConstraints & constraints, const hp::MappingCollection &mapping_collection) { // L2-projection based interpolation formed in one (in 2D) or two (in 3D) @@ -6295,12 +6266,12 @@ namespace VectorTools template void project_boundary_values_curl_conforming_l2( - const DoFHandler & dof_handler, - const unsigned int first_vector_component, - const Function & boundary_function, - const types::boundary_id boundary_component, - ConstraintMatrix & constraints, - const Mapping & mapping) + const DoFHandler & dof_handler, + const unsigned int first_vector_component, + const Function & boundary_function, + const types::boundary_id boundary_component, + AffineConstraints &constraints, + const Mapping & mapping) { // non-hp version - calls the internal // compute_project_boundary_values_curl_conforming_l2() function @@ -6323,7 +6294,7 @@ namespace VectorTools const unsigned int first_vector_component, const Function & boundary_function, const types::boundary_id boundary_component, - ConstraintMatrix & constraints, + AffineConstraints & constraints, const hp::MappingCollection &mapping_collection) { // hp version - calls the internal @@ -6352,7 +6323,7 @@ namespace VectorTools const unsigned int first_vector_component, const Function<2> & boundary_function, const std::vector> &jacobians, - ConstraintMatrix & constraints) + AffineConstraints & constraints) { // Compute the integral over the product of the normal components of // the boundary function times the normal components of the shape @@ -6403,7 +6374,7 @@ namespace VectorTools cell->face(face)->get_dof_indices(face_dof_indices, cell->active_fe_index()); - // Copy the computed values in the ConstraintMatrix only, if the degree + // Copy the computed values in the AffineConstraints only, if the degree // of freedom is not already constrained. for (unsigned int i = 0; i < fe.dofs_per_face; ++i) if (!(constraints.is_constrained(face_dof_indices[i]))) @@ -6425,7 +6396,7 @@ namespace VectorTools const unsigned int, const Function &, const std::vector> &, - ConstraintMatrix &) + AffineConstraints &) { Assert(false, ExcNotImplemented()); } @@ -6532,12 +6503,12 @@ namespace VectorTools template void project_boundary_values_div_conforming( - const DoFHandler & dof_handler, - const unsigned int first_vector_component, - const Function & boundary_function, - const types::boundary_id boundary_component, - ConstraintMatrix & constraints, - const Mapping & mapping) + const DoFHandler & dof_handler, + const unsigned int first_vector_component, + const Function & boundary_function, + const types::boundary_id boundary_component, + AffineConstraints &constraints, + const Mapping & mapping) { const unsigned int spacedim = dim; // Interpolate the normal components @@ -6636,23 +6607,13 @@ namespace VectorTools case 3: { - // In three dimensions the - // edges between two faces - // are treated twice. - // Therefore we store the - // computed values in a - // vector and copy them over - // in the ConstraintMatrix - // after all values have been - // computed. - // If we have two values for - // one edge, we choose the one, - // which was computed with the - // higher order element. - // If both elements are of the - // same order, we just keep the - // first value and do not - // compute a second one. + // In three dimensions the edges between two faces are treated + // twice. Therefore we store the computed values in a vector + // and copy them over in the AffineConstraints after all values + // have been computed. If we have two values for one edge, we + // choose the one, which was computed with the higher order + // element. If both elements are of the same order, we just + // keep the first value and do not compute a second one. const unsigned int & n_dofs = dof_handler.n_dofs(); std::vector dof_values(n_dofs); std::vector projected_dofs(n_dofs); @@ -6670,12 +6631,9 @@ namespace VectorTools ++face) if (cell->face(face)->boundary_id() == boundary_component) { - // This is only - // implemented, if the - // FE is a Raviart-Thomas - // element. If the FE is - // a FESystem we cannot - // check this. + // This is only implemented, if the FE is a + // Raviart-Thomas element. If the FE is a FESystem we + // cannot check this. if (dynamic_cast *>( &cell->get_fe()) == nullptr) { @@ -6733,7 +6691,7 @@ namespace VectorTools const unsigned int first_vector_component, const Function & boundary_function, const types::boundary_id boundary_component, - ConstraintMatrix & constraints, + AffineConstraints & constraints, const hp::MappingCollection &mapping_collection) { const unsigned int spacedim = dim; @@ -6900,7 +6858,7 @@ namespace VectorTools const DoFHandlerType &dof_handler, const unsigned int first_vector_component, const std::set & boundary_ids, - ConstraintMatrix & constraints, + AffineConstraints & constraints, const Mapping & mapping) { ZeroFunction zero_function(dim); @@ -6923,7 +6881,7 @@ namespace VectorTools const unsigned int first_vector_component, const std::set & boundary_ids, typename FunctionMap::type &function_map, - ConstraintMatrix & constraints, + AffineConstraints & constraints, const Mapping & mapping) { Assert(dim > 1, @@ -7506,7 +7464,7 @@ namespace VectorTools const DoFHandlerType &dof_handler, const unsigned int first_vector_component, const std::set & boundary_ids, - ConstraintMatrix & constraints, + AffineConstraints & constraints, const Mapping & mapping) { ZeroFunction zero_function(dim); @@ -7529,10 +7487,11 @@ namespace VectorTools const unsigned int first_vector_component, const std::set & boundary_ids, typename FunctionMap::type &function_map, - ConstraintMatrix & constraints, + AffineConstraints & constraints, const Mapping & mapping) { - ConstraintMatrix no_normal_flux_constraints(constraints.get_local_lines()); + AffineConstraints no_normal_flux_constraints( + constraints.get_local_lines()); compute_nonzero_normal_flux_constraints(dof_handler, first_vector_component, boundary_ids, diff --git a/source/numerics/vector_tools_boundary.inst.in b/source/numerics/vector_tools_boundary.inst.in index d93bce9546..0ede9e500f 100644 --- a/source/numerics/vector_tools_boundary.inst.in +++ b/source/numerics/vector_tools_boundary.inst.in @@ -93,7 +93,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS; const DH &, const std::map *> &, - ConstraintMatrix &, + AffineConstraints &, const ComponentMask &); template void @@ -102,7 +102,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS; const DH &, const types::boundary_id, const Function &, - ConstraintMatrix &, + AffineConstraints &, const ComponentMask &); template void @@ -110,7 +110,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS; const DH &, const types::boundary_id, const Function &, - ConstraintMatrix &, + AffineConstraints &, const ComponentMask &); template void @@ -118,7 +118,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS; const DH &, const std::map *> &, - ConstraintMatrix &, + AffineConstraints &, const ComponentMask &); \} #endif @@ -177,7 +177,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) const std::map *> &, const Quadrature &, - ConstraintMatrix &, + AffineConstraints &, std::vector); template void @@ -186,7 +186,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) const std::map *> &, const Quadrature &, - ConstraintMatrix &, + AffineConstraints &, std::vector); # if deal_II_dimension != 1 @@ -196,7 +196,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) const unsigned int, const Function &, const types::boundary_id, - ConstraintMatrix &, + AffineConstraints &, const Mapping &); template void project_boundary_values_curl_conforming( @@ -204,7 +204,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) const unsigned int, const Function &, const types::boundary_id, - ConstraintMatrix &, + AffineConstraints &, const hp::MappingCollection &); template void project_boundary_values_curl_conforming_l2( @@ -212,7 +212,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) const unsigned int, const Function &, const types::boundary_id, - ConstraintMatrix &, + AffineConstraints &, const Mapping &); template void project_boundary_values_curl_conforming_l2( @@ -220,7 +220,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) const unsigned int, const Function &, const types::boundary_id, - ConstraintMatrix &, + AffineConstraints &, const hp::MappingCollection &); template void project_boundary_values_div_conforming( @@ -228,7 +228,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) const unsigned int, const Function &, const types::boundary_id, - ConstraintMatrix &, + AffineConstraints &, const Mapping &); template void project_boundary_values_div_conforming( @@ -236,7 +236,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) const unsigned int, const Function &, const types::boundary_id, - ConstraintMatrix &, + AffineConstraints &, const hp::MappingCollection &); # endif #endif diff --git a/source/numerics/vector_tools_constraints.inst.in b/source/numerics/vector_tools_constraints.inst.in index eb6833d4fc..2d6b205c1e 100644 --- a/source/numerics/vector_tools_constraints.inst.in +++ b/source/numerics/vector_tools_constraints.inst.in @@ -30,7 +30,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) const unsigned int first_vector_component, const std::set & boundary_ids, FunctionMap::type &function_map, - ConstraintMatrix & constraints, + AffineConstraints & constraints, const Mapping & mapping); template void @@ -39,7 +39,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) const unsigned int first_vector_component, const std::set & boundary_ids, FunctionMap::type & function_map, - ConstraintMatrix & constraints, + AffineConstraints & constraints, const Mapping & mapping); template void @@ -48,7 +48,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) const unsigned int first_vector_component, const std::set & boundary_ids, FunctionMap::type &function_map, - ConstraintMatrix & constraints, + AffineConstraints & constraints, const Mapping & mapping); template void compute_nonzero_tangential_flux_constraints( @@ -56,7 +56,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) const unsigned int first_vector_component, const std::set & boundary_ids, FunctionMap::type & function_map, - ConstraintMatrix & constraints, + AffineConstraints & constraints, const Mapping & mapping); template void @@ -64,7 +64,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) const DoFHandler &dof_handler, const unsigned int first_vector_component, const std::set & boundary_ids, - ConstraintMatrix & constraints, + AffineConstraints & constraints, const Mapping & mapping); template void @@ -72,7 +72,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) const hp::DoFHandler &dof_handler, const unsigned int first_vector_component, const std::set & boundary_ids, - ConstraintMatrix & constraints, + AffineConstraints & constraints, const Mapping & mapping); template void @@ -80,7 +80,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) const DoFHandler &dof_handler, const unsigned int first_vector_component, const std::set & boundary_ids, - ConstraintMatrix & constraints, + AffineConstraints & constraints, const Mapping & mapping); template void @@ -88,7 +88,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) const hp::DoFHandler &dof_handler, const unsigned int first_vector_component, const std::set & boundary_ids, - ConstraintMatrix & constraints, + AffineConstraints & constraints, const Mapping & mapping); # endif #endif diff --git a/source/numerics/vector_tools_interpolate.inst.in b/source/numerics/vector_tools_interpolate.inst.in index 59d5ab794c..53d8fef3dd 100644 --- a/source/numerics/vector_tools_interpolate.inst.in +++ b/source/numerics/vector_tools_interpolate.inst.in @@ -122,7 +122,7 @@ for (VEC : VECTOR_TYPES; deal_II_dimension : DIMENSIONS; const DoFHandler &, const VEC &, const DoFHandler &, - const ConstraintMatrix &, + const AffineConstraints &, VEC &); template void @@ -130,7 +130,7 @@ for (VEC : VECTOR_TYPES; deal_II_dimension : DIMENSIONS; const InterGridMap< DoFHandler> &, const VEC &, - const ConstraintMatrix &, + const AffineConstraints &, VEC &); \} #endif diff --git a/source/numerics/vector_tools_project.inst.in b/source/numerics/vector_tools_project.inst.in index f71c34a37e..cf73984404 100644 --- a/source/numerics/vector_tools_project.inst.in +++ b/source/numerics/vector_tools_project.inst.in @@ -25,7 +25,7 @@ for (VEC : REAL_VECTOR_TYPES; deal_II_dimension : DIMENSIONS; project( const Mapping &, const DoFHandler &, - const ConstraintMatrix &, + const AffineConstraints &, const Quadrature &, const Function &, VEC &, @@ -36,7 +36,7 @@ for (VEC : REAL_VECTOR_TYPES; deal_II_dimension : DIMENSIONS; template void project( const DoFHandler &, - const ConstraintMatrix &, + const AffineConstraints &, const Quadrature &, const Function &, VEC &, diff --git a/source/numerics/vector_tools_project_codim.inst.in b/source/numerics/vector_tools_project_codim.inst.in index fc51d36cc9..a3c87c3475 100644 --- a/source/numerics/vector_tools_project_codim.inst.in +++ b/source/numerics/vector_tools_project_codim.inst.in @@ -25,7 +25,7 @@ for (VEC : REAL_VECTOR_TYPES; deal_II_dimension : DIMENSIONS; project( const Mapping &, const DoFHandler &, - const ConstraintMatrix &, + const AffineConstraints &, const Quadrature &, const Function &, VEC &, @@ -36,7 +36,7 @@ for (VEC : REAL_VECTOR_TYPES; deal_II_dimension : DIMENSIONS; template void project( const DoFHandler &, - const ConstraintMatrix &, + const AffineConstraints &, const Quadrature &, const Function &, VEC &, diff --git a/source/numerics/vector_tools_project_hp.inst.in b/source/numerics/vector_tools_project_hp.inst.in index 9421dd4b1b..96f2259bf4 100644 --- a/source/numerics/vector_tools_project_hp.inst.in +++ b/source/numerics/vector_tools_project_hp.inst.in @@ -26,7 +26,7 @@ for (VEC : REAL_VECTOR_TYPES; deal_II_dimension : DIMENSIONS; const hp::MappingCollection &, const hp::DoFHandler &, - const ConstraintMatrix &, + const AffineConstraints &, const hp::QCollection &, const Function &, VEC &, @@ -37,7 +37,7 @@ for (VEC : REAL_VECTOR_TYPES; deal_II_dimension : DIMENSIONS; template void project( const hp::DoFHandler &, - const ConstraintMatrix &, + const AffineConstraints &, const hp::QCollection &, const Function &, VEC &, diff --git a/source/numerics/vector_tools_project_qp.inst.in b/source/numerics/vector_tools_project_qp.inst.in index 767255686c..0cc4547371 100644 --- a/source/numerics/vector_tools_project_qp.inst.in +++ b/source/numerics/vector_tools_project_qp.inst.in @@ -22,7 +22,7 @@ for (VEC : REAL_NONBLOCK_VECTORS; deal_II_dimension : DIMENSIONS) project( const Mapping &, const DoFHandler &, - const ConstraintMatrix &, + const AffineConstraints &, const Quadrature &, const std::function( std::shared_ptr> - matrix_free, - const ConstraintMatrix &constraints, + matrix_free, + const AffineConstraints &constraints, const std::function(const unsigned int, const unsigned int)> &, VEC &, @@ -33,8 +33,8 @@ for (VEC : REAL_NONBLOCK_VECTORS; deal_II_dimension : DIMENSIONS) template void project( std::shared_ptr> - matrix_free, - const ConstraintMatrix &constraints, + matrix_free, + const AffineConstraints &constraints, const unsigned int, const std::function(const unsigned int, const unsigned int)> &, diff --git a/source/numerics/vector_tools_rhs.inst.in b/source/numerics/vector_tools_rhs.inst.in index 54d58792ba..1288f5e5b1 100644 --- a/source/numerics/vector_tools_rhs.inst.in +++ b/source/numerics/vector_tools_rhs.inst.in @@ -27,7 +27,7 @@ for (VEC : REAL_VECTOR_TYPES; deal_II_dimension : DIMENSIONS; const Quadrature &, const Function &, VEC &, - const ConstraintMatrix &); + const AffineConstraints &); template void create_right_hand_side( @@ -35,7 +35,7 @@ for (VEC : REAL_VECTOR_TYPES; deal_II_dimension : DIMENSIONS; const Quadrature &, const Function &, VEC &, - const ConstraintMatrix &); + const AffineConstraints &); \} #endif } @@ -58,7 +58,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) const hp::QCollection &, const Function::value_type> &, Vector &, - const ConstraintMatrix &); + const AffineConstraints &); template void create_right_hand_side &, const Function::value_type> &, Vector &, - const ConstraintMatrix &); + const AffineConstraints &); # if deal_II_dimension > 1 template void