From: Wolfgang Bangerth Date: Fri, 8 Mar 2024 06:18:12 +0000 (-0700) Subject: Implement VectorTools::interpolate_to_finer/coarser_mesh(). X-Git-Tag: v9.6.0-rc1~489^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b874a87ca73a63dbc6d14660c017d59637f10ecb;p=dealii.git Implement VectorTools::interpolate_to_finer/coarser_mesh(). --- diff --git a/include/deal.II/numerics/vector_tools_interpolate.h b/include/deal.II/numerics/vector_tools_interpolate.h index 9880f4ae4e..8407108a24 100644 --- a/include/deal.II/numerics/vector_tools_interpolate.h +++ b/include/deal.II/numerics/vector_tools_interpolate.h @@ -227,7 +227,11 @@ namespace VectorTools * partitioned in their own ways, will not typically have corresponding * cells owned by the same process, and implementing the interpolation * procedure would require transfering data between processes in ways - * that are difficult to implement efficiently. + * that are difficult to implement efficiently. However, some special + * cases can more easily be implemented, namely the case where one + * of the meshes is strictly coarser or finer than the other. For these + * cases, see the interpolate_to_coarser_mesh() and + * interpolate_to_finer_mesh(). * * @dealiiConceptRequires{concepts::is_writable_dealii_vector_type} */ @@ -278,6 +282,71 @@ namespace VectorTools const AffineConstraints &constraints, VectorType &u2); + /** + * This function addresses one of the limitations of the + * interpolate_to_different_mesh() function, namely that the latter only + * works on parallel triangulations in very specific cases where both + * triangulations involved happen to be partitioned in such a way that + * if a process owns a cell of one mesh, it also needs to own the + * corresponding parent of child cells of the other mesh. In practice, this is + * rarely the case. + * + * This function does not have this restriction, and consequently also works + * in parallel, as long as the mesh we interpolate onto can be obtained from + * the mesh we interpolate off by *coarsening*. In other words, the target + * mesh is coarser than the source mesh. + * + * The function takes an additional constraints argument that is used after + * interpolation to ensure that the output vector is conforming (that is, + * that all entries of the output vector conform to their constraints). + * @p constraints_coarse therefore needs to correspond to + * @p dof_handler_coarse . + * + * The opposite operation, interpolation from a coarser to a finer mesh, + * is implemented in the interpolate_to_finer_mesh() function. + */ + template + DEAL_II_CXX20_REQUIRES(concepts::is_writable_dealii_vector_type) + void interpolate_to_coarser_mesh( + const DoFHandler &dof_handler_fine, + const VectorType &u_fine, + const DoFHandler &dof_handler_coarse, + const AffineConstraints + &constraints_coarse, + VectorType &u_coarse); + + /** + * This function addresses one of the limitations of the + * interpolate_to_different_mesh() function, namely that the latter only + * works on parallel triangulations in very specific cases where both + * triangulations involved happen to be partitioned in such a way that + * if a process owns a cell of one mesh, it also needs to own the + * corresponding parent of child cells of the other mesh. In practice, this is + * rarely the case. + * + * This function does not have this restriction, and consequently also works + * in parallel, as long as the mesh we interpolate onto can be obtained from + * the mesh we interpolate off by *refinement*. In other words, the target + * mesh is finer than the source mesh. + * + * The function takes an additional constraints argument that is used after + * interpolation to ensure that the output vector is conforming (that is, + * that all entries of the output vector conform to their constraints). + * @p constraints_finee therefore needs to correspond to + * @p dof_handler_fine . + * + * The opposite operation, interpolation from a finer to a coarser mesh, + * is implemented in the interpolate_to_coarser_mesh() function. + */ + template + DEAL_II_CXX20_REQUIRES(concepts::is_writable_dealii_vector_type) + void interpolate_to_finer_mesh( + const DoFHandler &dof_handler_coarse, + const VectorType &u_coarse, + const DoFHandler &dof_handler_fine, + const AffineConstraints &constraints_fine, + VectorType &u_fine); + /** @} */ /** diff --git a/include/deal.II/numerics/vector_tools_interpolate.templates.h b/include/deal.II/numerics/vector_tools_interpolate.templates.h index 63c2285941..4902dff684 100644 --- a/include/deal.II/numerics/vector_tools_interpolate.templates.h +++ b/include/deal.II/numerics/vector_tools_interpolate.templates.h @@ -16,6 +16,8 @@ #define dealii_vector_tools_interpolate_templates_h +#include + #include #include @@ -33,10 +35,21 @@ #include #include +#include + #include DEAL_II_NAMESPACE_OPEN +namespace LinearAlgebra +{ + namespace EpetraWrappers + { + class Vector; + } +} // namespace LinearAlgebra + + namespace VectorTools { // This namespace contains the actual implementation called @@ -989,6 +1002,248 @@ namespace VectorTools // Apply hanging node constraints. constraints.distribute(u2); } + + + // Some helper functions. Much of the functions here is taken from a similar + // implementation in data_out_dof_data.templates.h. + namespace InterpolateBetweenMeshes + { + template + void + create_vector(const DoFHandler &dof_handler, + LinearAlgebra::distributed::Vector &u) + { + const IndexSet &locally_owned_dofs = dof_handler.locally_owned_dofs(); + + const IndexSet locally_relevant_dofs = + DoFTools::extract_locally_relevant_dofs(dof_handler); + + Assert(locally_owned_dofs.is_contiguous(), + ExcMessage("You are trying to use vectors with non-contiguous " + "locally-owned index sets. This is not possible.")); + + u.reinit(locally_owned_dofs, + locally_relevant_dofs, + dof_handler.get_communicator()); + } + + + /** + * Copy the data from an arbitrary non-block vector to a + * LinearAlgebra::distributed::Vector. + */ + template + void + copy_locally_owned_data_from( + const VectorType &src, + LinearAlgebra::distributed::Vector &dst) + { + if constexpr (!IsBlockVector::value) + { + // If source and destination vector have the same underlying scalar, + // we can directly import elements by using only one temporary vector: + if constexpr (std::is_same_v) + { + LinearAlgebra::ReadWriteVector + temp; + temp.reinit(src.locally_owned_elements()); + temp.import_elements(src, VectorOperation::insert); + + dst.import_elements(temp, VectorOperation::insert); + } + else + // The source and destination vector have different scalar types. We + // need to split the parallel import and local copy operations into + // two phases + { + LinearAlgebra::ReadWriteVector + temp; + temp.reinit(src.locally_owned_elements()); + temp.import_elements(src, VectorOperation::insert); + + LinearAlgebra::ReadWriteVector temp2; + temp2.reinit(temp, true); + temp2 = temp; + + dst.import_elements(temp2, VectorOperation::insert); + } + } + else + DEAL_II_NOT_IMPLEMENTED(); + } + +#ifdef DEAL_II_WITH_TRILINOS + template + void + copy_locally_owned_data_from( + const TrilinosWrappers::MPI::Vector &src, + LinearAlgebra::distributed::Vector &dst) + { + // ReadWriteVector does not work for ghosted + // TrilinosWrappers::MPI::Vector objects. Fall back to copy the + // entries manually. + for (const auto i : dst.locally_owned_elements()) + dst[i] = src[i]; + } +#endif + +#ifdef DEAL_II_TRILINOS_WITH_TPETRA + template + void + copy_locally_owned_data_from( + const LinearAlgebra::TpetraWrappers::Vector &src, + LinearAlgebra::distributed::Vector &dst) + { + // ReadWriteVector does not work for ghosted + // TrilinosWrappers::MPI::Vector objects. Fall back to copy the + // entries manually. + for (const auto i : dst.locally_owned_elements()) + dst[i] = src[i]; + } +#endif + } // namespace InterpolateBetweenMeshes + + + template + DEAL_II_CXX20_REQUIRES(concepts::is_writable_dealii_vector_type) + void interpolate_to_coarser_mesh( + const DoFHandler &dof_handler_fine, + const VectorType &u_fine, + const DoFHandler &dof_handler_coarse, + const AffineConstraints + &constraints_coarse, + VectorType &u_coarse) + { + Assert(GridTools::have_same_coarse_mesh(dof_handler_coarse, + dof_handler_fine), + ExcMessage("The two DoF handlers must represent triangulations that " + "have the same coarse meshes")); + AssertDimension(dof_handler_fine.n_dofs(), u_fine.size()); + AssertDimension(dof_handler_coarse.n_dofs(), u_coarse.size()); + + using LAVector = + LinearAlgebra::distributed::Vector; + + // Create and initialize a version of the coarse vector using the + // p::d::Vector type: + LAVector my_u_coarse; + InterpolateBetweenMeshes::create_vector(dof_handler_coarse, my_u_coarse); + + // Then do the same for the fine vector, and initialize it with a copy + // of the source vector + LAVector my_u_fine; + InterpolateBetweenMeshes::create_vector(dof_handler_fine, my_u_fine); + InterpolateBetweenMeshes::copy_locally_owned_data_from(u_fine, my_u_fine); + my_u_fine.update_ghost_values(); + + + // Set up transfer operator. The transfer object also wants a constraints + // object for the fine level, because it implements both the up- and + // down-transfers. But we here only ever need the down-transfer, which + // never touches the fine constraints. So we can use a dummy object for + // that. + AffineConstraints constraints_fine( + dof_handler_fine.locally_owned_dofs(), + DoFTools::extract_locally_relevant_dofs(dof_handler_fine)); + constraints_fine.close(); + + MGTwoLevelTransfer transfer; + transfer.reinit(dof_handler_fine, + dof_handler_coarse, + constraints_fine, + constraints_coarse); + + // Then perform the interpolation from fine to coarse mesh: + transfer.interpolate(my_u_coarse, my_u_fine); + + // Finally, apply the constraints to make sure that the resulting + // object is conforming. Then copy it back into the user vector. + // (At the time of writing, LinearAlgebra::EpetraWrappers::Vector + // does not allow individual element access via operator() or + // operator[]. So disallow this vector type for the moment here.) + constraints_coarse.distribute(my_u_coarse); + if constexpr (!std::is_same_v) + { + for (const auto i : u_coarse.locally_owned_elements()) + u_coarse(i) = my_u_coarse(i); + } + else + DEAL_II_NOT_IMPLEMENTED(); + } + + + + template + DEAL_II_CXX20_REQUIRES(concepts::is_writable_dealii_vector_type) + void interpolate_to_finer_mesh( + const DoFHandler &dof_handler_coarse, + const VectorType &u_coarse, + const DoFHandler &dof_handler_fine, + const AffineConstraints &constraints_fine, + VectorType &u_fine) + { + Assert(GridTools::have_same_coarse_mesh(dof_handler_fine, + dof_handler_coarse), + ExcMessage("The two DoF handlers must represent triangulations that " + "have the same coarse meshes")); + AssertDimension(dof_handler_fine.n_dofs(), u_fine.size()); + AssertDimension(dof_handler_coarse.n_dofs(), u_coarse.size()); + + using LAVector = + LinearAlgebra::distributed::Vector; + + // Create and initialize a version of the coarse vector using the + // p::d::Vector type. Copy the source vector + LAVector my_u_coarse; + InterpolateBetweenMeshes::create_vector(dof_handler_coarse, my_u_coarse); + InterpolateBetweenMeshes::copy_locally_owned_data_from(u_coarse, + my_u_coarse); + my_u_coarse.update_ghost_values(); + + // Then do the same for the fine vector, and initialize it with a copy + // of the source vector + LAVector my_u_fine; + InterpolateBetweenMeshes::create_vector(dof_handler_fine, my_u_fine); + + // Set up transfer operator. The transfer object also wants a constraints + // object for the coarse level, because it implements both the up- and + // down-transfers. But we here only ever need the up-transfer, which + // never touches the coarse constraints. So we can use a dummy object for + // that. + AffineConstraints constraints_coarse( + dof_handler_coarse.locally_owned_dofs(), + DoFTools::extract_locally_relevant_dofs(dof_handler_coarse)); + constraints_coarse.close(); + + MGTwoLevelTransfer transfer; + transfer.reinit(dof_handler_fine, + dof_handler_coarse, + constraints_fine, + constraints_coarse); + + // Then perform the interpolation from coarse to fine mesh. (Note that + // we add to my_u_fine, but that vector is initially a zero vector, + // so we are really just writing into it.) + transfer.prolongate_and_add(my_u_fine, my_u_coarse); + + // Finally, apply the constraints to make sure that the resulting + // object is conforming. Then copy it back into the user vector. + // (At the time of writing, LinearAlgebra::EpetraWrappers::Vector + // does not allow individual element access via operator() or + // operator[]. So disallow this vector type for the moment here.) + constraints_fine.distribute(my_u_fine); + if constexpr (!std::is_same_v) + { + for (const auto i : u_fine.locally_owned_elements()) + u_fine(i) = my_u_fine(i); + } + else + DEAL_II_NOT_IMPLEMENTED(); + } + + } // namespace VectorTools DEAL_II_NAMESPACE_CLOSE diff --git a/source/numerics/vector_tools_interpolate.inst.in b/source/numerics/vector_tools_interpolate.inst.in index 592964c707..a2500e3592 100644 --- a/source/numerics/vector_tools_interpolate.inst.in +++ b/source/numerics/vector_tools_interpolate.inst.in @@ -127,3 +127,28 @@ for (VEC : VECTOR_TYPES; deal_II_dimension : DIMENSIONS; \} #endif } + +for (VEC : REAL_VECTOR_TYPES; deal_II_dimension : DIMENSIONS; + deal_II_space_dimension : SPACE_DIMENSIONS) + { +#if deal_II_dimension == deal_II_space_dimension + namespace VectorTools + \{ + template void + interpolate_to_coarser_mesh( + const DoFHandler &, + const VEC &, + const DoFHandler &, + const AffineConstraints &, + VEC &); + + template void + interpolate_to_finer_mesh( + const DoFHandler &, + const VEC &, + const DoFHandler &, + const AffineConstraints &, + VEC &); + \} +#endif + }