From f5dd5329f6e63b40c87622e19d2b90c0ba1b7d2b Mon Sep 17 00:00:00 2001 From: Bruno Turcksin Date: Tue, 14 Sep 2021 18:17:50 +0000 Subject: [PATCH] Add FieldTransfer for the transfer of field with FE_Nothing --- include/deal.II/distributed/field_transfer.h | 113 +++++++++ source/distributed/CMakeLists.txt | 2 + source/distributed/field_transfer.cc | 242 +++++++++++++++++++ source/distributed/field_transfer.inst.in | 39 +++ 4 files changed, 396 insertions(+) create mode 100644 include/deal.II/distributed/field_transfer.h create mode 100644 source/distributed/field_transfer.cc create mode 100644 source/distributed/field_transfer.inst.in diff --git a/include/deal.II/distributed/field_transfer.h b/include/deal.II/distributed/field_transfer.h new file mode 100644 index 0000000000..6b5ef36508 --- /dev/null +++ b/include/deal.II/distributed/field_transfer.h @@ -0,0 +1,113 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + +#ifndef dealii_distributed_field_transfer_h +#define dealii_distributed_field_transfer_h + +#include + +#include + +#include + +#include + +DEAL_II_NAMESPACE_OPEN + +namespace parallel +{ + namespace distributed + { + /** + * Classes and functions in the experimental namespace can be modified or + * removed without being deprecated first. + */ + namespace experimental + { + /** + * This class is similar to SolutionTransfer but it supports the case + * where elements have been activated during refinement, i.e., FE_Nothing + * elements have been associated with a finite elements during refinement. + */ + template + class FieldTransfer + { + private: + using Number = typename VectorType::value_type; + + public: + /** + * 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. + */ + 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 + * @param[in] fe_nothing_index The finite element index associated with + * FE_Nothing + */ + void + prepare_for_coarsening_and_refinement( + const VectorType & in, + const unsigned int fe_nothing_index); + + /** + * Interpolate the data previously stored in this object before the mesh + * was refined or coarsenend 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. + */ + void + interpolate(const Number & new_value, + const AffineConstraints &affine_constraints, + VectorType & out); + + private: + /** + * DoFHandler associated with the object. + */ + const DoFHandler &dof_handler; + + /** + * Data transfered by cell_data_transfer. + */ + std::vector> data_to_transfer; + + /** + * CellDataTransfer used to perform the field transfer. + */ + std::unique_ptr< + CellDataTransfer>>> + cell_data_transfer; + }; + } // namespace experimental + } // namespace distributed +} // namespace parallel + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/source/distributed/CMakeLists.txt b/source/distributed/CMakeLists.txt index 836602ef8d..7754496a0e 100644 --- a/source/distributed/CMakeLists.txt +++ b/source/distributed/CMakeLists.txt @@ -26,6 +26,7 @@ SET(_unity_include_src tria_base.cc shared_tria.cc p4est_wrappers.cc + field_transfer.cc ) SET(_separate_src @@ -44,6 +45,7 @@ SET(_inst grid_refinement.inst.in cell_weights.inst.in cell_data_transfer.inst.in + field_transfer.inst.in fully_distributed_tria.inst.in repartitioning_policy_tools.inst.in solution_transfer.inst.in diff --git a/source/distributed/field_transfer.cc b/source/distributed/field_transfer.cc new file mode 100644 index 0000000000..75b2d32319 --- /dev/null +++ b/source/distributed/field_transfer.cc @@ -0,0 +1,242 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + +#include + +#ifdef DEAL_II_WITH_P4EST + +# include + + +DEAL_II_NAMESPACE_OPEN + +namespace parallel +{ + namespace distributed + { + namespace experimental + { + template + FieldTransfer::FieldTransfer( + const DoFHandler &dof) + : dof_handler(dof) + { + cell_data_transfer = std::make_unique< + CellDataTransfer>>>( + dynamic_cast< + dealii::parallel::distributed::Triangulation &>( + const_cast &>( + dof_handler.get_triangulation()))); + } + + + + template + void + FieldTransfer:: + prepare_for_coarsening_and_refinement( + const VectorType & in, + const unsigned int fe_nothing_index) + { + const unsigned int dofs_per_cell = + dof_handler.get_fe_collection().max_dofs_per_cell(); + + Vector cell_solution(dofs_per_cell); + Vector dummy_cell_solution(dofs_per_cell); + + for (auto &val : dummy_cell_solution) + { + val = std::numeric_limits::infinity(); + } + + in.update_ghost_values(); + + std::vector dof_indices(dofs_per_cell); + for (const auto &cell : dof_handler.active_cell_iterators()) + { + if ((cell->is_locally_owned()) && + (cell->active_fe_index() != fe_nothing_index)) + { + cell->get_dof_values(in, cell_solution); + data_to_transfer.push_back(cell_solution); + } + else + { + data_to_transfer.push_back(dummy_cell_solution); + } + } + + cell_data_transfer->prepare_for_coarsening_and_refinement( + data_to_transfer); + } + + + + template + void + FieldTransfer::interpolate( + const Number & new_value, + const AffineConstraints &affine_constraints, + VectorType & out) + { + const unsigned int dofs_per_cell = + dof_handler.get_fe_collection().max_dofs_per_cell(); + std::vector> transferred_data( + dof_handler.get_triangulation().n_active_cells(), + Vector(dofs_per_cell)); + + cell_data_transfer->unpack(transferred_data); + + // Free the memory allocated by data_to_transfer + data_to_transfer.clear(); + + for (unsigned int i = 0; i < out.locally_owned_size(); ++i) + out.local_element(i) = std::numeric_limits::infinity(); + + unsigned int cell_i = 0; + for (auto const &cell : dof_handler.active_cell_iterators()) + { + if ((cell->is_locally_owned()) && + (transferred_data[cell_i][0] != + std::numeric_limits::infinity())) + { + cell->set_dof_values(transferred_data[cell_i], out); + } + ++cell_i; + } + + + // Communicate the results. + out.compress(dealii::VectorOperation::min); + + // Treat hanging nodes + std::vector dof_indices; + std::vector dofs_map; + std::vector>> + constraint_lines; + std::vector constraint_values; + IndexSet locally_owned_dofs = dof_handler.locally_owned_dofs(); + for (auto constrained_dof : locally_owned_dofs) + if (affine_constraints.is_constrained(constrained_dof)) + { + auto *constraint = + affine_constraints.get_constraint_entries(constrained_dof); + const unsigned int line_size = constraint->size(); + bool add_line = false; + for (unsigned int i = 0; i < line_size; ++i) + { + types::global_dof_index constraining_dof = + (*constraint)[i].first; + // If one of the constraining value is infinity, we need to + // reverse the relationship + if (out[constraining_dof] == + std::numeric_limits::infinity()) + { + add_line = true; + break; + } + } + + if (add_line) + { + std::vector> line; + Number val = out[constrained_dof]; + + for (unsigned int i = 0; i < line_size; ++i) + { + types::global_dof_index constraining_dof = + (*constraint)[i].first; + + if (out[constraining_dof] == + std::numeric_limits::infinity()) + { + auto constraining_dof_map_it = + std::find(dofs_map.begin(), + dofs_map.end(), + constraining_dof); + if (constraining_dof_map_it == dofs_map.end()) + { + dofs_map.push_back(constraining_dof); + } + line.push_back((*constraint)[i]); + } + else + { + val -= + out[constraining_dof] * (*constraint)[i].second; + } + } + constraint_lines.push_back(line); + constraint_values.push_back(val); + } + } + + // Build a constraint matrix that we invert + const unsigned int n_rows = constraint_lines.size(); + if (n_rows > 0) + { + const unsigned int n_cols = dofs_map.size(); + dealii::LAPACKFullMatrix constraints_matrix(n_rows, n_cols); + dealii::Vector constraint_values_vec(n_rows); + + for (unsigned int i = 0; i < n_rows; ++i) + { + for (unsigned int j = 0; j < n_cols; ++j) + { + if (j < constraint_lines[i].size()) + { + auto constraint_it = + std::find(dofs_map.begin(), + dofs_map.end(), + constraint_lines[i][j].first); + constraints_matrix(i, + constraint_it - dofs_map.begin()) = + constraint_lines[i][j].second; + } + } + + constraint_values_vec[i] = constraint_values[i]; + } + + constraints_matrix.compute_inverse_svd(); + + dealii::Vector new_constrained_values(n_cols); + constraints_matrix.vmult(new_constrained_values, + constraint_values_vec); + + for (unsigned int i = 0; i < n_cols; ++i) + { + out[dofs_map[i]] = new_constrained_values[i]; + } + } + + // Set the value to the newly create DoFs. + std::for_each(out.begin(), out.end(), [&](Number &val) { + if (val == std::numeric_limits::infinity()) + { + val = new_value; + } + }); + } + } // namespace experimental + } // namespace distributed +} // namespace parallel + +// explicit instantiations +# include "field_transfer.inst" + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/source/distributed/field_transfer.inst.in b/source/distributed/field_transfer.inst.in new file mode 100644 index 0000000000..3dcbc6219b --- /dev/null +++ b/source/distributed/field_transfer.inst.in @@ -0,0 +1,39 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2022 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. +// +// --------------------------------------------------------------------- + + + +for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) + { + namespace parallel + \{ + namespace distributed + \{ + namespace experimental + \{ +#if deal_II_dimension <= deal_II_space_dimension + template class FieldTransfer< + deal_II_dimension, + ::dealii::LinearAlgebra::distributed::Vector, + deal_II_space_dimension>; + template class FieldTransfer< + deal_II_dimension, + ::dealii::LinearAlgebra::distributed::Vector, + deal_II_space_dimension>; +#endif + \} + \} + \} + } -- 2.39.5