From: Daniel Arndt Date: Fri, 4 Nov 2016 10:30:59 +0000 (+0100) Subject: Allow a template VectorType in VectorTools::create_right_hand_side and VectorTools... X-Git-Tag: v8.5.0-rc1~384^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=795dcdf42807ccc05f7dd47c3b6c6c48135793cf;p=dealii.git Allow a template VectorType in VectorTools::create_right_hand_side and VectorTools::create_boundary_right_hand_side --- diff --git a/doc/news/changes.h b/doc/news/changes.h index 1a2aee72cb..c2fc5e7416 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -38,6 +38,13 @@ inconvenience this causes.

    +
  1. Changed: VectorTools::create_right_hand_side and + VectorTools::create_boundary_right_hand_side now take an additional template + parameter VectorType. +
    + (Daniel Arndt, 2016/10/25) +
    2. +
  2. Deprecated: ParameterHandler::read_input, ParameterHandler::read_input_from_xml, and ParameterHandler::read_input_from_string are now deprecated in favor of @@ -495,7 +502,7 @@ inconvenience this causes. (Rajat Arora, 2016/10/29)
    3. -
  3. New: Add MatrixFreeOperators::MassOperator representing a mass matrix. +
  4. New: Add MatrixFreeOperators::MassOperator representing a mass matrix.
    (Daniel Arndt, 2016/10/27)
    5. @@ -515,6 +522,12 @@ inconvenience this causes. (Martin Kronbichler, 2016/10/27) +
  5. >New: VectorTools::create_right_hand_side can be used for parallel + Triangulations and homogeneous constraints using a ConstraintMatrix. +
    + (Daniel Arndt, 2016/10/25) +
    6. +
  6. Fixed: GridIn::read_vtk() had off-by-one errors in reading face boundary indicators from VTK files. Consequently, not all boundary indicators were correctly set. diff --git a/include/deal.II/numerics/vector_tools.h b/include/deal.II/numerics/vector_tools.h index 7071137f61..8509ed417b 100644 --- a/include/deal.II/numerics/vector_tools.h +++ b/include/deal.II/numerics/vector_tools.h @@ -739,7 +739,6 @@ namespace VectorTools const ConstraintMatrix &constraints, VectorType &u2); - /** * The same function as above, but takes an InterGridMap object directly as * a parameter. Useful for interpolating several vectors at the same time. @@ -1925,41 +1924,45 @@ namespace VectorTools * * See the general documentation of this namespace for further information. */ - template + template void create_right_hand_side (const Mapping &mapping, const DoFHandler &dof, const Quadrature &q, - const Function &rhs, - Vector &rhs_vector); + const Function &rhs, + VectorType &rhs_vector, + const ConstraintMatrix &constraints=ConstraintMatrix()); /** * Call the create_right_hand_side() function, see above, with * mapping=MappingQGeneric@(1). */ - template + template void create_right_hand_side (const DoFHandler &dof, const Quadrature &q, - const Function &rhs, - Vector &rhs_vector); + const Function &rhs, + VectorType &rhs_vector, + const ConstraintMatrix &constraints=ConstraintMatrix()); /** * Like the previous set of functions, but for hp objects. */ - template + template void create_right_hand_side (const hp::MappingCollection &mapping, const hp::DoFHandler &dof, const hp::QCollection &q, - const Function &rhs, - Vector &rhs_vector); + const Function &rhs, + VectorType &rhs_vector, + const ConstraintMatrix &constraints=ConstraintMatrix()); /** * Like the previous set of functions, but for hp objects. */ - template + template void create_right_hand_side (const hp::DoFHandler &dof, const hp::QCollection &q, - const Function &rhs, - Vector &rhs_vector); + const Function &rhs, + VectorType &rhs_vector, + const ConstraintMatrix &constraints=ConstraintMatrix()); /** * Create a right hand side vector for a point source at point @p p. In @@ -2069,12 +2072,12 @@ namespace VectorTools * @see * @ref GlossBoundaryIndicator "Glossary entry on boundary indicators" */ - template - void create_boundary_right_hand_side (const Mapping &mapping, - const DoFHandler &dof, - const Quadrature &q, - const Function &rhs, - Vector &rhs_vector, + template + void create_boundary_right_hand_side (const Mapping &mapping, + const DoFHandler &dof, + const Quadrature &q, + const Function &rhs, + VectorType &rhs_vector, const std::set &boundary_ids = std::set()); /** @@ -2084,12 +2087,12 @@ namespace VectorTools * @see * @ref GlossBoundaryIndicator "Glossary entry on boundary indicators" */ - template + template void create_boundary_right_hand_side (const DoFHandler &dof, const Quadrature &q, - const Function &rhs, - Vector &rhs_vector, + const Function &rhs, + VectorType &rhs_vector, const std::set &boundary_ids = std::set()); /** @@ -2098,13 +2101,13 @@ namespace VectorTools * @see * @ref GlossBoundaryIndicator "Glossary entry on boundary indicators" */ - template + template void create_boundary_right_hand_side (const hp::MappingCollection &mapping, const hp::DoFHandler &dof, const hp::QCollection &q, - const Function &rhs, - Vector &rhs_vector, + const Function &rhs, + VectorType &rhs_vector, const std::set &boundary_ids = std::set()); /** @@ -2115,12 +2118,12 @@ namespace VectorTools * @see * @ref GlossBoundaryIndicator "Glossary entry on boundary indicators" */ - template + template void create_boundary_right_hand_side (const hp::DoFHandler &dof, const hp::QCollection &q, - const Function &rhs, - Vector &rhs_vector, + const Function &rhs, + VectorType &rhs_vector, const std::set &boundary_ids = std::set()); //@} diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h index 6f04aa2ace..b6bfa29e04 100644 --- a/include/deal.II/numerics/vector_tools.templates.h +++ b/include/deal.II/numerics/vector_tools.templates.h @@ -777,6 +777,8 @@ namespace VectorTools return true; } + + template void invert_mass_matrix(const SparseMatrix &mass_matrix, const Vector &rhs, @@ -804,6 +806,7 @@ namespace VectorTools } + /** * Generic implementation of the project() function */ @@ -1178,6 +1181,7 @@ namespace VectorTools } + template void project (const DoFHandler &dof, const ConstraintMatrix &constraints, @@ -1228,13 +1232,17 @@ namespace VectorTools } - template + + template void create_right_hand_side (const Mapping &mapping, const DoFHandler &dof_handler, const Quadrature &quadrature, - const Function &rhs_function, - Vector &rhs_vector) + const Function &rhs_function, + VectorType &rhs_vector, + const ConstraintMatrix &constraints) { + typedef typename VectorType::value_type Number; + const FiniteElement &fe = dof_handler.get_fe(); Assert (fe.n_components() == rhs_function.n_components, ExcDimensionMismatch(fe.n_components(), rhs_function.n_components)); @@ -1252,7 +1260,7 @@ namespace VectorTools n_components = fe.n_components(); std::vector dofs (dofs_per_cell); - Vector cell_vector (dofs_per_cell); + Vector cell_vector (dofs_per_cell); typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), @@ -1260,105 +1268,107 @@ namespace VectorTools if (n_components==1) { - std::vector rhs_values(n_q_points); + std::vector rhs_values(n_q_points); for (; cell!=endc; ++cell) - { - fe_values.reinit(cell); - - const std::vector &weights = fe_values.get_JxW_values (); - rhs_function.value_list (fe_values.get_quadrature_points(), - rhs_values); + if (cell->is_locally_owned()) + { + fe_values.reinit(cell); - cell_vector = 0; - for (unsigned int point=0; point &weights = fe_values.get_JxW_values (); + rhs_function.value_list (fe_values.get_quadrature_points(), + rhs_values); - cell->get_dof_indices (dofs); + cell_vector = 0; + for (unsigned int point=0; pointget_dof_indices(dofs); + constraints.distribute_local_to_global(cell_vector, dofs, rhs_vector); + } } else { - std::vector > rhs_values(n_q_points, - Vector(n_components)); + std::vector > rhs_values(n_q_points, + Vector(n_components)); for (; cell!=endc; ++cell) - { - fe_values.reinit(cell); - - const std::vector &weights = fe_values.get_JxW_values (); - rhs_function.vector_value_list (fe_values.get_quadrature_points(), - rhs_values); + if (cell->is_locally_owned()) + { + fe_values.reinit(cell); - cell_vector = 0; - // Use the faster code if the - // FiniteElement is primitive - if (fe.is_primitive ()) - { - for (unsigned int point=0; point &weights = fe_values.get_JxW_values (); + rhs_function.vector_value_list (fe_values.get_quadrature_points(), + rhs_values); - cell_vector(i) += rhs_values[point](component) * - fe_values.shape_value(i,point) * - weights[point]; - } - } - else - { - // Otherwise do it the way - // proposed for vector valued - // elements - for (unsigned int point=0; pointget_dof_indices (dofs); + cell_vector(i) += rhs_values[point](component) * + fe_values.shape_value(i,point) * + weights[point]; + } + } + else + { + // Otherwise do it the way + // proposed for vector valued + // elements + for (unsigned int point=0; pointget_dof_indices (dofs); - for (unsigned int i=0; i + template void create_right_hand_side (const DoFHandler &dof_handler, const Quadrature &quadrature, - const Function &rhs_function, - Vector &rhs_vector) + const Function &rhs_function, + VectorType &rhs_vector, + const ConstraintMatrix &constraints) { create_right_hand_side(StaticMappingQ1::mapping, dof_handler, quadrature, - rhs_function, rhs_vector); + rhs_function, rhs_vector, constraints); } - template + template void create_right_hand_side (const hp::MappingCollection &mapping, const hp::DoFHandler &dof_handler, const hp::QCollection &quadrature, - const Function &rhs_function, - Vector &rhs_vector) + const Function &rhs_function, + VectorType &rhs_vector, + const ConstraintMatrix &constraints) { + typedef typename VectorType::value_type Number; + const hp::FECollection &fe = dof_handler.get_fe(); Assert (fe.n_components() == rhs_function.n_components, ExcDimensionMismatch(fe.n_components(), rhs_function.n_components)); @@ -1374,7 +1384,7 @@ namespace VectorTools const unsigned int n_components = fe.n_components(); std::vector dofs (fe.max_dofs_per_cell()); - Vector cell_vector (fe.max_dofs_per_cell()); + Vector cell_vector (fe.max_dofs_per_cell()); typename hp::DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), @@ -1382,115 +1392,115 @@ namespace VectorTools if (n_components==1) { - std::vector rhs_values; + std::vector rhs_values; for (; cell!=endc; ++cell) - { - x_fe_values.reinit(cell); + if (cell->is_locally_owned()) + { + x_fe_values.reinit(cell); - const FEValues &fe_values = x_fe_values.get_present_fe_values(); + const FEValues &fe_values = x_fe_values.get_present_fe_values(); - const unsigned int dofs_per_cell = fe_values.dofs_per_cell, - n_q_points = fe_values.n_quadrature_points; - rhs_values.resize (n_q_points); - dofs.resize (dofs_per_cell); - cell_vector.reinit (dofs_per_cell); + const unsigned int dofs_per_cell = fe_values.dofs_per_cell, + n_q_points = fe_values.n_quadrature_points; + rhs_values.resize (n_q_points); + dofs.resize (dofs_per_cell); + cell_vector.reinit (dofs_per_cell); - const std::vector &weights = fe_values.get_JxW_values (); - rhs_function.value_list (fe_values.get_quadrature_points(), - rhs_values); + const std::vector &weights = fe_values.get_JxW_values (); + rhs_function.value_list (fe_values.get_quadrature_points(), + rhs_values); - cell_vector = 0; - for (unsigned int point=0; pointget_dof_indices (dofs); + cell->get_dof_indices (dofs); - for (unsigned int i=0; i > rhs_values; + std::vector > rhs_values; for (; cell!=endc; ++cell) - { - x_fe_values.reinit(cell); + if (cell->is_locally_owned()) + { + x_fe_values.reinit(cell); - const FEValues &fe_values = x_fe_values.get_present_fe_values(); + const FEValues &fe_values = x_fe_values.get_present_fe_values(); - const unsigned int dofs_per_cell = fe_values.dofs_per_cell, - n_q_points = fe_values.n_quadrature_points; - rhs_values.resize (n_q_points, - Vector(n_components)); - dofs.resize (dofs_per_cell); - cell_vector.reinit (dofs_per_cell); + const unsigned int dofs_per_cell = fe_values.dofs_per_cell, + n_q_points = fe_values.n_quadrature_points; + rhs_values.resize (n_q_points, + Vector(n_components)); + dofs.resize (dofs_per_cell); + cell_vector.reinit (dofs_per_cell); - const std::vector &weights = fe_values.get_JxW_values (); - rhs_function.vector_value_list (fe_values.get_quadrature_points(), - rhs_values); + const std::vector &weights = fe_values.get_JxW_values (); + rhs_function.vector_value_list (fe_values.get_quadrature_points(), + rhs_values); - cell_vector = 0; + cell_vector = 0; - // Use the faster code if the - // FiniteElement is primitive - if (cell->get_fe().is_primitive ()) - { - for (unsigned int point=0; pointget_fe().system_to_component_index(i).first; + // Use the faster code if the + // FiniteElement is primitive + if (cell->get_fe().is_primitive ()) + { + for (unsigned int point=0; pointget_fe().system_to_component_index(i).first; - cell_vector(i) += rhs_values[point](component) * - fe_values.shape_value(i,point) * - weights[point]; - } - } - else - { - // Otherwise do it the way proposed - // for vector valued elements - for (unsigned int point=0; pointget_fe().get_nonzero_components(i)[comp_i]) - { - cell_vector(i) += rhs_values[point](comp_i) * - fe_values.shape_value_component(i,point,comp_i) * - weights[point]; - } - } + cell_vector(i) += rhs_values[point](component) * + fe_values.shape_value(i,point) * + weights[point]; + } + } + else + { + // Otherwise do it the way proposed + // for vector valued elements + for (unsigned int point=0; pointget_fe().get_nonzero_components(i)[comp_i]) + { + cell_vector(i) += rhs_values[point](comp_i) * + fe_values.shape_value_component(i,point,comp_i) * + weights[point]; + } + } - cell->get_dof_indices (dofs); + cell->get_dof_indices (dofs); - for (unsigned int i=0; i + template void create_right_hand_side (const hp::DoFHandler &dof_handler, const hp::QCollection &quadrature, - const Function &rhs_function, - Vector &rhs_vector) + const Function &rhs_function, + VectorType &rhs_vector, + const ConstraintMatrix &constraints) { create_right_hand_side(hp::StaticMappingQ1::mapping_collection, dof_handler, quadrature, - rhs_function, rhs_vector); + rhs_function, rhs_vector, constraints); } - template void create_point_source_vector (const Mapping &mapping, const DoFHandler &dof_handler, @@ -1678,13 +1688,13 @@ namespace VectorTools - template + template void create_boundary_right_hand_side (const Mapping &mapping, const DoFHandler &dof_handler, const Quadrature &quadrature, - const Function &rhs_function, - Vector &rhs_vector, + const Function &rhs_function, + VectorType &rhs_vector, const std::set &boundary_ids) { const FiniteElement &fe = dof_handler.get_fe(); @@ -1801,12 +1811,12 @@ namespace VectorTools - template + template void create_boundary_right_hand_side (const DoFHandler &dof_handler, const Quadrature &quadrature, - const Function &rhs_function, - Vector &rhs_vector, + const Function &rhs_function, + VectorType &rhs_vector, const std::set &boundary_ids) { create_boundary_right_hand_side(StaticMappingQ1::mapping, dof_handler, @@ -1817,13 +1827,13 @@ namespace VectorTools - template + template void create_boundary_right_hand_side (const hp::MappingCollection &mapping, const hp::DoFHandler &dof_handler, const hp::QCollection &quadrature, - const Function &rhs_function, - Vector &rhs_vector, + const Function &rhs_function, + VectorType &rhs_vector, const std::set &boundary_ids) { const hp::FECollection &fe = dof_handler.get_fe(); @@ -1952,12 +1962,12 @@ namespace VectorTools - template + template void create_boundary_right_hand_side (const hp::DoFHandler &dof_handler, const hp::QCollection &quadrature, - const Function &rhs_function, - Vector &rhs_vector, + const Function &rhs_function, + VectorType &rhs_vector, const std::set &boundary_ids) { create_boundary_right_hand_side(hp::StaticMappingQ1::mapping_collection, @@ -2495,7 +2505,7 @@ namespace VectorTools // 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 - ReductionControl control(5*rhs.size(), 0., 1.e-12, false, false); + ReductionControl control(5.*rhs.size(), 0., 1.e-12, false, false); GrowingVectorMemory > memory; SolverCG > cg(control,memory); diff --git a/source/numerics/vector_tools_rhs.inst.in b/source/numerics/vector_tools_rhs.inst.in index 2c38f5a6d2..2de583e733 100644 --- a/source/numerics/vector_tools_rhs.inst.in +++ b/source/numerics/vector_tools_rhs.inst.in @@ -14,25 +14,26 @@ // --------------------------------------------------------------------- -for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) +for (VEC : REAL_SERIAL_VECTORS; deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) { #if deal_II_dimension <= deal_II_space_dimension namespace VectorTools \{ template - void create_right_hand_side - (const Mapping &, - const DoFHandler &, - const Quadrature &, - const Function &, - Vector &); + void create_right_hand_side + (const Mapping &, + const DoFHandler &, + const Quadrature &, + const Function &, + VEC &, + const ConstraintMatrix &); template - void create_right_hand_side - (const DoFHandler &, - const Quadrature &, - const Function &, - Vector &); - + void create_right_hand_side + (const DoFHandler &, + const Quadrature &, + const Function &, + VEC &, + const ConstraintMatrix &); \} #endif } @@ -45,61 +46,63 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS #if deal_II_dimension == deal_II_space_dimension template - void create_right_hand_side - (const hp::MappingCollection &, - const hp::DoFHandler &, - const hp::QCollection &, - const Function &, - Vector &); + void create_right_hand_side > + (const hp::MappingCollection &, + const hp::DoFHandler &, + const hp::QCollection &, + const Function::value_type> &, + Vector &, + const ConstraintMatrix &); template - void create_right_hand_side - (const hp::DoFHandler &, + void create_right_hand_side > + (const hp::DoFHandler &, const hp::QCollection &, - const Function &, - Vector &); + const Function::value_type> &, + Vector &, + const ConstraintMatrix &); #if deal_II_dimension > 1 template void - create_boundary_right_hand_side - (const Mapping &, - const DoFHandler &, + create_boundary_right_hand_side > + (const Mapping &, + const DoFHandler &, const Quadrature &, - const Function &, - Vector &, - const std::set &); + const Function::value_type> &, + Vector &, + const std::set &); #endif template void - create_boundary_right_hand_side - (const DoFHandler &, + create_boundary_right_hand_side > + (const DoFHandler &, const Quadrature &, - const Function &, - Vector &, - const std::set &); + const Function::value_type> &, + Vector &, + const std::set &); #if deal_II_dimension > 1 template void - create_boundary_right_hand_side - (const hp::MappingCollection &, - const hp::DoFHandler &, - const hp::QCollection &, - const Function &, - Vector &, - const std::set &); + create_boundary_right_hand_side > + (const hp::MappingCollection &, + const hp::DoFHandler &, + const hp::QCollection &, + const Function::value_type> &, + Vector &, + const std::set &); #endif template void - create_boundary_right_hand_side - (const hp::DoFHandler &, + create_boundary_right_hand_side > + (const hp::DoFHandler &, const hp::QCollection &, - const Function &, - Vector &, - const std::set &); + const Function::value_type> &, + Vector &, + const std::set &); #endif