From d42cd5f0318d5fb04f04a6861483cdff581ab945 Mon Sep 17 00:00:00 2001 From: Sebastian Proell Date: Tue, 14 Feb 2023 19:10:22 +0100 Subject: [PATCH] FieldTransfer: enable coarsening --- include/deal.II/distributed/field_transfer.h | 14 +- source/distributed/field_transfer.cc | 86 ++++++++- tests/hp/field_transfer_05.cc | 170 ++++++++++++++++++ ...ransfer_05.with_p4est=true.mpirun=2.output | 27 +++ 4 files changed, 286 insertions(+), 11 deletions(-) create mode 100644 tests/hp/field_transfer_05.cc create mode 100644 tests/hp/field_transfer_05.with_p4est=true.mpirun=2.output diff --git a/include/deal.II/distributed/field_transfer.h b/include/deal.II/distributed/field_transfer.h index 0793e4fcce..95c0ec8f8f 100644 --- a/include/deal.II/distributed/field_transfer.h +++ b/include/deal.II/distributed/field_transfer.h @@ -52,17 +52,11 @@ namespace parallel * Constructor. * * @param[in] dof_handler The DoFHandler on which all the operations - * will happen. At the time when this constructor is call, the - * DoFHandler still points to the Triangulation before the refinement in - * question happens. + * will happen. This constructor must be called before the underlying + * Triangulation is coarsened/refined. */ FieldTransfer(const DoFHandler &dof_handler); - /** - * Destructor. - */ - ~FieldTransfer() = default; - /** * Prepare the current object for coarsening and refinement. * @param[in] in The vector that will be interpolated @@ -76,10 +70,10 @@ namespace parallel /** * Interpolate the data previously stored in this object before the mesh - * was refined or coarsenend onto the current set of cells. @p new_value + * was refined or coarsened onto the current set of cells. @p new_value * is the value associated to the new degrees of freedom that where * created during the element activation. @p affine_constraints is the - * AffinieConstraints after refinement. + * AffineConstraints after refinement. */ void interpolate(const Number & new_value, diff --git a/source/distributed/field_transfer.cc b/source/distributed/field_transfer.cc index 2c8366d38a..2f633debc2 100644 --- a/source/distributed/field_transfer.cc +++ b/source/distributed/field_transfer.cc @@ -34,12 +34,96 @@ namespace parallel const DoFHandler &dof) : dof_handler(dof) { + // When coarsening, we want to mimic the behavior of SolutionTransfer + // and interpolate from child cells to parent. Define this strategy here + // since it is not readily available + const auto coarsening_strategy = + [this]( + const typename dealii::Triangulation::cell_iterator + & parent, + const std::vector> &children_values) { + // get the equivalent DoFCellAccessor + typename DoFHandler::cell_iterator dof_cell_iterator( + &dof_handler.get_triangulation(), + parent->level(), + parent->index(), + &dof_handler); + + int fe_index = 0; + if (dof_handler.has_hp_capabilities()) + fe_index = dealii::internal::hp::DoFHandlerImplementation:: + dominated_future_fe_on_children( + dof_cell_iterator); + + const auto &fe = dof_handler.get_fe(fe_index); + Assert(fe.n_dofs_per_cell() > 0, + ExcMessage( + "Cannot coarsen onto a FiniteElement with no DoFs.")); + AssertDimension(dof_cell_iterator->n_children(), + children_values.size()); + + const auto child_iterators = dof_cell_iterator->child_iterators(); + const unsigned int n_children_with_fe_nothing = + std::count_if(child_iterators.begin(), + child_iterators.end(), + [](const auto &child_cell) { + return child_cell->get_fe().n_dofs_per_cell() == + 0; + }); + + Assert( + n_children_with_fe_nothing == 0 || + n_children_with_fe_nothing == dof_cell_iterator->n_children(), + ExcMessage( + "Coarsening is only supported for parent cells where either all" + " or none of the child cells are FE_Nothing.")); + + // in case all children are FE_Nothing there is nothing to + // interpolate and we just return the first entry from the children + // values (containing invalid entries) + if (n_children_with_fe_nothing == dof_cell_iterator->n_children()) + { + return children_values[0]; + } + + const unsigned int dofs_per_cell = fe.n_dofs_per_cell(); + Vector tmp(dofs_per_cell); + Vector interpolated_values(dofs_per_cell); + + // Otherwise, perform the actual interpolation here. Due to the + // assert above, we know that all child cells have data to + // interpolate. + for (unsigned int child = 0; + child < dof_cell_iterator->n_children(); + ++child) + { + // interpolate the previously stored values on a child to the + // mother cell + fe.get_restriction_matrix(child, + dof_cell_iterator->refinement_case()) + .vmult(tmp, children_values[child]); + + // and add up or set them in the output vector + for (unsigned int i = 0; i < dofs_per_cell; ++i) + if (fe.restriction_is_additive(i)) + interpolated_values(i) += tmp(i); + else if (tmp(i) != Number()) + interpolated_values(i) = tmp(i); + } + + return interpolated_values; + }; + cell_data_transfer = std::make_unique< CellDataTransfer>>>( dynamic_cast< dealii::parallel::distributed::Triangulation &>( const_cast &>( - dof_handler.get_triangulation()))); + dof_handler.get_triangulation())), + false, + &dealii::AdaptationStrategies::Refinement:: + preserve>, + coarsening_strategy); } diff --git a/tests/hp/field_transfer_05.cc b/tests/hp/field_transfer_05.cc new file mode 100644 index 0000000000..c48f5f97f9 --- /dev/null +++ b/tests/hp/field_transfer_05.cc @@ -0,0 +1,170 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2021 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +// Check that we can use FieldTransfer when coarsening cells with non-uniform +// data + +#include + +#include +#include + +#include +#include + +#include +#include +#include + +#include + +#include +#include + +#include "../tests.h" + + +class TestFunction : public Function<2> +{ +public: + double + value(const dealii::Point<2> &p, const unsigned int component) const override + { + return std::exp(10 - 10 * p[1]); + } +}; + +void +output(const DoFHandler<2> & dof_handler, + const LinearAlgebra::distributed::Vector &vector, + const std::string & file_name) +{ + DataOut<2> data_out; + data_out.attach_dof_handler(dof_handler); + data_out.add_data_vector(vector, "solution"); + data_out.build_patches(); + data_out.write_vtu_in_parallel(file_name, MPI_COMM_WORLD); +} + +int +main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + + MPILogInitAll mpi_log; + + parallel::distributed::Triangulation<2> triangulation(MPI_COMM_WORLD); + GridGenerator::hyper_cube(triangulation); + triangulation.refine_global(3); + + hp::FECollection<2> fe_collection; + fe_collection.push_back(FE_Q<2>(1)); + fe_collection.push_back(FE_Nothing<2>()); + + DoFHandler<2> dof_handler(triangulation); + + // Assign FE_Nothing to half of the domain + for (auto cell : dof_handler.active_cell_iterators()) + { + if (cell->is_locally_owned()) + { + if (cell->center()[0] < 0.5) + { + cell->set_active_fe_index(0); + } + else + { + cell->set_active_fe_index(1); + } + } + } + + dof_handler.distribute_dofs(fe_collection); + + // Initialize solution + auto locally_relevant_dofs = + DoFTools::extract_locally_relevant_dofs(dof_handler); + LinearAlgebra::distributed::Vector solution( + dof_handler.locally_owned_dofs(), locally_relevant_dofs, MPI_COMM_WORLD); + + // Initialize with a non-constant function + VectorTools::interpolate(dof_handler, TestFunction(), solution); + + { + if (false) + output(dof_handler, solution, "result_0.vtu"); + + deallog << "solution before: " << std::endl; + std::stringstream ss; + solution.print(ss); + deallog << ss.str() << std::endl; + } + + parallel::distributed::experimental:: + FieldTransfer<2, LinearAlgebra::distributed::Vector> + field_transfer(dof_handler); + // Assign FE_Q to all the cells + for (auto cell : dof_handler.active_cell_iterators()) + { + if (cell->is_locally_owned()) + { + cell->set_future_fe_index(0); + } + } + + // coarsen some cells + for (auto cell : dof_handler.active_cell_iterators()) + { + if (cell->is_locally_owned()) + { + if (cell->center()[0] < 0.25 && cell->center()[1] < 0.25) + { + cell->set_coarsen_flag(); + } + } + } + + triangulation.prepare_coarsening_and_refinement(); + + field_transfer.prepare_for_coarsening_and_refinement(solution, 1); + + triangulation.execute_coarsening_and_refinement(); + + dof_handler.distribute_dofs(fe_collection); + + locally_relevant_dofs = DoFTools::extract_locally_relevant_dofs(dof_handler); + LinearAlgebra::distributed::Vector new_solution( + dof_handler.locally_owned_dofs(), locally_relevant_dofs, MPI_COMM_WORLD); + + AffineConstraints affine_constraints; + DoFTools::make_hanging_node_constraints(dof_handler, affine_constraints); + affine_constraints.close(); + + const double new_value = 11.; + + field_transfer.interpolate(new_value, affine_constraints, new_solution); + + affine_constraints.distribute(new_solution); + + { + if (false) + output(dof_handler, new_solution, "result_1.vtu"); + + deallog << "solution after coarsening and activation: " << std::endl; + std::stringstream ss; + new_solution.print(ss); + deallog << ss.str() << std::endl; + } +} diff --git a/tests/hp/field_transfer_05.with_p4est=true.mpirun=2.output b/tests/hp/field_transfer_05.with_p4est=true.mpirun=2.output new file mode 100644 index 0000000000..5bde4cb74b --- /dev/null +++ b/tests/hp/field_transfer_05.with_p4est=true.mpirun=2.output @@ -0,0 +1,27 @@ + +DEAL:0::solution before: +DEAL:0::Process #0 +Local range: [0, 25), global size: 45 +Vector data: +2.203e+04 2.203e+04 6.311e+03 6.311e+03 2.203e+04 6.311e+03 1.808e+03 1.808e+03 1.808e+03 2.203e+04 6.311e+03 2.203e+04 6.311e+03 1.808e+03 1.808e+03 5.180e+02 5.180e+02 5.180e+02 1.484e+02 1.484e+02 1.484e+02 5.180e+02 5.180e+02 1.484e+02 1.484e+02 + +DEAL:0::solution after coarsening and activation: +DEAL:0::Process #0 +Local range: [0, 42), global size: 78 +Vector data: +2.203e+04 2.203e+04 1.808e+03 1.808e+03 2.203e+04 1.192e+04 6.311e+03 2.203e+04 6.311e+03 1.808e+03 1.808e+03 1.808e+03 5.180e+02 5.180e+02 5.180e+02 1.484e+02 1.484e+02 1.484e+02 5.180e+02 5.180e+02 1.484e+02 1.484e+02 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 + + +DEAL:1::solution before: +DEAL:1::Process #1 +Local range: [25, 45), global size: 45 +Vector data: +4.252e+01 4.252e+01 4.252e+01 1.218e+01 1.218e+01 1.218e+01 4.252e+01 4.252e+01 1.218e+01 1.218e+01 3.490e+00 3.490e+00 3.490e+00 1.000e+00 1.000e+00 1.000e+00 3.490e+00 3.490e+00 1.000e+00 1.000e+00 + +DEAL:1::solution after coarsening and activation: +DEAL:1::Process #1 +Local range: [42, 78), global size: 78 +Vector data: +4.252e+01 4.252e+01 4.252e+01 1.218e+01 1.218e+01 1.218e+01 4.252e+01 4.252e+01 1.218e+01 1.218e+01 3.490e+00 3.490e+00 3.490e+00 1.000e+00 1.000e+00 1.000e+00 3.490e+00 3.490e+00 1.000e+00 1.000e+00 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 + + -- 2.39.5