From 4120eecf1bbfd42dc724e17af0d3e53368649dc1 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 26 May 2025 17:20:50 -0600 Subject: [PATCH] Move some code to a .templates.h file. --- .../numerics/solution_transfer.templates.h | 1135 +++++++++++++++++ source/numerics/solution_transfer.cc | 1112 +--------------- 2 files changed, 1136 insertions(+), 1111 deletions(-) create mode 100644 include/deal.II/numerics/solution_transfer.templates.h diff --git a/include/deal.II/numerics/solution_transfer.templates.h b/include/deal.II/numerics/solution_transfer.templates.h new file mode 100644 index 0000000000..b01392efc2 --- /dev/null +++ b/include/deal.II/numerics/solution_transfer.templates.h @@ -0,0 +1,1135 @@ +// ------------------------------------------------------------------------ +// +// SPDX-License-Identifier: LGPL-2.1-or-later +// Copyright (C) 1999 - 2023 by the deal.II authors +// +// This file is part of the deal.II library. +// +// Part of the source code is dual licensed under Apache-2.0 WITH +// LLVM-exception OR LGPL-2.1-or-later. Detailed license information +// governing the source code and code contributions can be found in +// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II. +// +// ------------------------------------------------------------------------ + +#ifndef dealii_solution_transfer_templates_h +#define dealii_solution_transfer_templates_h + +#include + +#include + +#include +#include + +#include +#include +#include + +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include + +#include +#include + + +DEAL_II_NAMESPACE_OPEN + + +namespace +{ + /** + * Optimized pack function for values assigned on degrees of freedom. + * + * Given that the elements of @p dof_values are stored in consecutive + * locations, we can just memcpy them. Since floating point values don't + * compress well, we also waive the compression that the default + * Utilities::pack() and Utilities::unpack() functions offer. + */ + template + std::vector + pack_dof_values(std::vector> &dof_values, + const unsigned int dofs_per_cell) + { + for (const auto &values : dof_values) + { + AssertDimension(values.size(), dofs_per_cell); + (void)values; + } + + const std::size_t bytes_per_entry = sizeof(value_type) * dofs_per_cell; + + std::vector buffer(dof_values.size() * bytes_per_entry); + for (unsigned int i = 0; i < dof_values.size(); ++i) + std::memcpy(&buffer[i * bytes_per_entry], + &dof_values[i](0), + bytes_per_entry); + + return buffer; + } + + + + /** + * Optimized unpack function for values assigned on degrees of freedom. + */ + template + std::vector> + unpack_dof_values( + const boost::iterator_range::const_iterator> &data_range, + const unsigned int dofs_per_cell) + { + const std::size_t bytes_per_entry = sizeof(value_type) * dofs_per_cell; + const unsigned int n_elements = data_range.size() / bytes_per_entry; + + Assert((data_range.size() % bytes_per_entry == 0), ExcInternalError()); + + std::vector> unpacked_data; + unpacked_data.reserve(n_elements); + for (unsigned int i = 0; i < n_elements; ++i) + { + Vector dof_values(dofs_per_cell); + std::memcpy(&dof_values(0), + &(*std::next(data_range.begin(), i * bytes_per_entry)), + bytes_per_entry); + unpacked_data.emplace_back(std::move(dof_values)); + } + + return unpacked_data; + } +} // namespace + + + +template +SolutionTransfer::SolutionTransfer( + const DoFHandler &dof, + const bool average_values) + : dof_handler(&dof, typeid(*this).name()) + , average_values(average_values) + , handle(numbers::invalid_unsigned_int) +{} + + + +template +void +SolutionTransfer:: + prepare_for_coarsening_and_refinement( + const std::vector &all_in) +{ + const dealii::internal::parallel::shared:: + TemporarilyRestoreSubdomainIds + subdomain_modifier(dof_handler->get_triangulation()); + + for (unsigned int i = 0; i < all_in.size(); ++i) + Assert(all_in[i]->size() == dof_handler->n_dofs(), + ExcDimensionMismatch(all_in[i]->size(), dof_handler->n_dofs())); + + input_vectors = all_in; + register_data_attach(); +} + + + +template +void +SolutionTransfer:: + prepare_for_coarsening_and_refinement(const std::vector &all_in) +{ + std::vector temp(all_in.size()); + + for (std::size_t i = 0; i < temp.size(); ++i) + temp[i] = &(all_in[i]); + + this->prepare_for_coarsening_and_refinement(temp); +} + + + +template +void +SolutionTransfer::register_data_attach() +{ + // TODO: casting away constness is bad + auto tria = const_cast *>( + &dof_handler->get_triangulation()); + Assert(tria != nullptr, ExcInternalError()); + + Assert(handle == numbers::invalid_unsigned_int, + ExcMessage("You can only add one solution per " + "SolutionTransfer object.")); + + handle = tria->register_data_attach( + [this](const typename Triangulation::cell_iterator &cell_, + const CellStatus status) { + return this->pack_callback(cell_, status); + }, + /*returns_variable_size_data=*/dof_handler->has_hp_capabilities()); +} + + + +template +void +SolutionTransfer:: + prepare_for_coarsening_and_refinement(const VectorType &in) +{ + std::vector all_in(1, &in); + prepare_for_coarsening_and_refinement(all_in); +} + + + +template +void +SolutionTransfer::prepare_for_serialization( + const VectorType &in) +{ + std::vector all_in(1, &in); + prepare_for_serialization(all_in); +} + + + +template +void +SolutionTransfer::prepare_for_serialization( + const std::vector &all_in) +{ + prepare_for_coarsening_and_refinement(all_in); +} + + + +template +void +SolutionTransfer::deserialize(VectorType &in) +{ + std::vector all_in(1, &in); + deserialize(all_in); +} + + + +template +void +SolutionTransfer::deserialize( + std::vector &all_in) +{ + register_data_attach(); + + // this makes interpolate() happy + input_vectors.resize(all_in.size()); + + interpolate(all_in); +} + + +template +void +SolutionTransfer::interpolate( + std::vector &all_out) +{ + const dealii::internal::parallel::shared:: + TemporarilyRestoreSubdomainIds + subdomain_modifier(dof_handler->get_triangulation()); + + Assert(input_vectors.size() == all_out.size(), + ExcDimensionMismatch(input_vectors.size(), all_out.size())); + for (unsigned int i = 0; i < all_out.size(); ++i) + Assert(all_out[i]->size() == dof_handler->n_dofs(), + ExcDimensionMismatch(all_out[i]->size(), dof_handler->n_dofs())); + + // TODO: casting away constness is bad + auto tria = const_cast *>( + &dof_handler->get_triangulation()); + Assert(tria != nullptr, ExcInternalError()); + Assert( + handle != numbers::invalid_unsigned_int, + ExcMessage( + "You can only call interpolate() once per SolutionTransfer object.")); + + using Number = typename VectorType::value_type; + + if (average_values) + for (auto *const vec : all_out) + *vec = Number(); + + VectorType valence; + + // initialize valence vector only if we need to average + if (average_values) + valence.reinit(*all_out[0]); + + tria->notify_ready_to_unpack( + handle, + [this, &all_out, &valence]( + const typename Triangulation::cell_iterator &cell_, + const CellStatus status, + const boost::iterator_range::const_iterator> + &data_range) { + this->unpack_callback(cell_, status, data_range, all_out, valence); + }); + + if (average_values) + { + // finalize valence: compress and invert + valence.compress(VectorOperation::add); + for (const auto i : valence.locally_owned_elements()) + valence[i] = (static_cast(valence[i]) == Number() ? + Number() : + (Number(1.0) / static_cast(valence[i]))); + valence.compress(VectorOperation::insert); + + for (auto *const vec : all_out) + { + // compress and weight with valence + vec->compress(VectorOperation::add); + vec->scale(valence); + } + } + else + { + for (auto *const vec : all_out) + vec->compress(VectorOperation::insert); + } + + input_vectors.clear(); + handle = numbers::invalid_unsigned_int; +} + + +template +void +SolutionTransfer::interpolate( + std::vector &all_out) +{ + std::vector temp(all_out.size()); + + for (std::size_t i = 0; i < temp.size(); ++i) + temp[i] = &(all_out[i]); + + this->interpolate(temp); +} + + + +template +void +SolutionTransfer::interpolate(VectorType &out) +{ + std::vector all_out(1, &out); + interpolate(all_out); +} + + + +template +std::vector +SolutionTransfer::pack_callback( + const typename Triangulation::cell_iterator &cell_, + const CellStatus status) +{ + typename DoFHandler::cell_iterator cell(*cell_, dof_handler); + + // create buffer for each individual object + std::vector<::dealii::Vector> dof_values( + input_vectors.size()); + + unsigned int fe_index = 0; + if (dof_handler->has_hp_capabilities()) + { + switch (status) + { + case CellStatus::cell_will_persist: + case CellStatus::cell_will_be_refined: + { + fe_index = cell->future_fe_index(); + break; + } + + case CellStatus::children_will_be_coarsened: + { + // In case of coarsening, we need to find a suitable FE index + // for the parent cell. We choose the 'least dominant fe' + // on all children from the associated FECollection. + if constexpr (running_in_debug_mode()) + { + for (const auto &child : cell->child_iterators()) + Assert(child->is_active() && child->coarsen_flag_set(), + typename dealii::Triangulation< + dim>::ExcInconsistentCoarseningFlags()); + } + + fe_index = dealii::internal::hp::DoFHandlerImplementation:: + dominated_future_fe_on_children(cell); + break; + } + + default: + Assert(false, ExcInternalError()); + break; + } + } + + const unsigned int dofs_per_cell = + dof_handler->get_fe(fe_index).n_dofs_per_cell(); + + if (dofs_per_cell == 0) + return std::vector(); // nothing to do for FE_Nothing + + auto it_input = input_vectors.cbegin(); + auto it_output = dof_values.begin(); + for (; it_input != input_vectors.cend(); ++it_input, ++it_output) + { + it_output->reinit(dofs_per_cell); + cell->get_interpolated_dof_values(*(*it_input), *it_output, fe_index); + } + + return pack_dof_values(dof_values, + dofs_per_cell); +} + + + +template +void +SolutionTransfer::unpack_callback( + const typename Triangulation::cell_iterator &cell_, + const CellStatus status, + const boost::iterator_range::const_iterator> &data_range, + std::vector &all_out, + VectorType &valence) +{ + typename DoFHandler::cell_iterator cell(*cell_, dof_handler); + + unsigned int fe_index = 0; + if (dof_handler->has_hp_capabilities()) + { + switch (status) + { + case CellStatus::cell_will_persist: + case CellStatus::children_will_be_coarsened: + { + fe_index = cell->active_fe_index(); + break; + } + + case CellStatus::cell_will_be_refined: + { + // After refinement, this particular cell is no longer active, + // and its children have inherited its FE index. However, to + // unpack the data on the old cell, we need to recover its FE + // index from one of the children. Just to be sure, we also + // check if all children have the same FE index. + fe_index = cell->child(0)->active_fe_index(); + for (unsigned int child_index = 1; + child_index < cell->n_children(); + ++child_index) + Assert(cell->child(child_index)->active_fe_index() == fe_index, + ExcInternalError()); + break; + } + + default: + Assert(false, ExcInternalError()); + break; + } + } + + const unsigned int dofs_per_cell = + dof_handler->get_fe(fe_index).n_dofs_per_cell(); + + if (dofs_per_cell == 0) + return; // nothing to do for FE_Nothing + + const std::vector<::dealii::Vector> + dof_values = + unpack_dof_values(data_range, + dofs_per_cell); + + // check if sizes match + AssertDimension(dof_values.size(), all_out.size()); + + // check if we have enough dofs provided by the FE object + // to interpolate the transferred data correctly + for (auto it_dof_values = dof_values.begin(); + it_dof_values != dof_values.end(); + ++it_dof_values) + Assert( + dofs_per_cell == it_dof_values->size(), + ExcMessage( + "The transferred data was packed with a different number of dofs than the " + "currently registered FE object assigned to the DoFHandler has.")); + + // distribute data for each registered vector on mesh + auto it_input = dof_values.cbegin(); + auto it_output = all_out.begin(); + for (; it_input != dof_values.cend(); ++it_input, ++it_output) + if (average_values) + cell->distribute_local_to_global_by_interpolation(*it_input, + *(*it_output), + fe_index); + else + cell->set_dof_values_by_interpolation(*it_input, + *(*it_output), + fe_index, + true); + + if (average_values) + { + // compute valence vector if averaging should be performed + Vector ones(dofs_per_cell); + ones = 1.0; + cell->distribute_local_to_global_by_interpolation(ones, + valence, + fe_index); + } +} + + + +template +void +SolutionTransfer::clear() +{ + // nothing to do +} + + + +namespace Legacy +{ + + template + SolutionTransfer::SolutionTransfer( + const DoFHandler &dof) + : dof_handler(&dof, typeid(*this).name()) + , n_dofs_old(0) + , prepared_for(none) + { + Assert((dynamic_cast< + const parallel::distributed::Triangulation *>( + &dof_handler->get_triangulation()) == nullptr), + ExcMessage( + "You are calling the dealii::SolutionTransfer class " + "with a DoFHandler that is built on a " + "parallel::distributed::Triangulation. This will not " + "work for parallel computations. You probably want to " + "use the parallel::distributed::SolutionTransfer class.")); + } + + + + template + SolutionTransfer::~SolutionTransfer() + { + clear(); + } + + + + template + void + SolutionTransfer::clear() + { + indices_on_cell.clear(); + dof_values_on_cell.clear(); + cell_map.clear(); + + prepared_for = none; + } + + + + template + void + SolutionTransfer::prepare_for_pure_refinement() + { + Assert(prepared_for != pure_refinement, ExcAlreadyPrepForRef()); + Assert(prepared_for != coarsening_and_refinement, + ExcAlreadyPrepForCoarseAndRef()); + + clear(); + + // We need to access dof indices on the entire domain. For + // parallel::shared::Triangulations, ownership of cells might change. If + // they allow artificial cells, we need to restore the "true" cell owners + // temporarily. + // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save + // the current set of subdomain ids, set subdomain ids to the "true" owner + // of each cell upon construction of the TemporarilyRestoreSubdomainIds + // object, and later restore these flags when it is destroyed. + const dealii::internal::parallel::shared:: + TemporarilyRestoreSubdomainIds + subdomain_modifier(dof_handler->get_triangulation()); + + const unsigned int n_active_cells = + dof_handler->get_triangulation().n_active_cells(); + n_dofs_old = dof_handler->n_dofs(); + + // efficient reallocation of indices_on_cell + std::vector>(n_active_cells) + .swap(indices_on_cell); + + for (const auto &cell : dof_handler->active_cell_iterators()) + { + const unsigned int i = cell->active_cell_index(); + indices_on_cell[i].resize(cell->get_fe().n_dofs_per_cell()); + // on each cell store the indices of the + // dofs. after refining we get the values + // on the children by taking these + // indices, getting the respective values + // out of the data vectors and prolonging + // them to the children + cell->get_dof_indices(indices_on_cell[i]); + cell_map[std::make_pair(cell->level(), cell->index())] = + Pointerstruct(&indices_on_cell[i], cell->active_fe_index()); + } + prepared_for = pure_refinement; + } + + + + template + void + SolutionTransfer::refine_interpolate( + const VectorType &in, + VectorType &out) const + { + Assert(prepared_for == pure_refinement, ExcNotPrepared()); + Assert(in.size() == n_dofs_old, + ExcDimensionMismatch(in.size(), n_dofs_old)); + Assert(out.size() == dof_handler->n_dofs(), + ExcDimensionMismatch(out.size(), dof_handler->n_dofs())); + Assert(&in != &out, + ExcMessage("Vectors cannot be used as input and output" + " at the same time!")); + + // We need to access dof indices on the entire domain. For + // parallel::shared::Triangulations, ownership of cells might change. If + // they allow artificial cells, we need to restore the "true" cell owners + // temporarily. + // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save + // the current set of subdomain ids, set subdomain ids to the "true" owner + // of each cell upon construction of the TemporarilyRestoreSubdomainIds + // object, and later restore these flags when it is destroyed. + const dealii::internal::parallel::shared:: + TemporarilyRestoreSubdomainIds + subdomain_modifier(dof_handler->get_triangulation()); + + Vector local_values(0); + + typename std::map, + Pointerstruct>::const_iterator pointerstruct, + cell_map_end = cell_map.end(); + + for (const auto &cell : dof_handler->cell_iterators()) + { + pointerstruct = + cell_map.find(std::make_pair(cell->level(), cell->index())); + + if (pointerstruct != cell_map_end) + // this cell was refined or not + // touched at all, so we can get + // the new values by just setting + // or interpolating to the children, + // which is both done by one + // function + { + const unsigned int this_fe_index = + pointerstruct->second.active_fe_index; + const unsigned int dofs_per_cell = + cell->get_dof_handler().get_fe(this_fe_index).n_dofs_per_cell(); + local_values.reinit(dofs_per_cell, true); + + // make sure that the size of the stored indices is the same as + // dofs_per_cell. since we store the desired fe_index, we know + // what this size should be + Assert(dofs_per_cell == (*pointerstruct->second.indices_ptr).size(), + ExcInternalError()); + for (unsigned int i = 0; i < dofs_per_cell; ++i) + local_values(i) = + dealii::internal::ElementAccess::get( + in, (*pointerstruct->second.indices_ptr)[i]); + cell->set_dof_values_by_interpolation(local_values, + out, + this_fe_index, + true); + } + } + } + + + + namespace internal + { + /** + * Generate a table that contains + * interpolation matrices between + * each combination of finite + * elements used in a DoFHandler of + * some kind. Since not all + * elements can be interpolated + * onto each other, the table may + * contain empty matrices for those + * combinations of elements for + * which no such interpolation is + * implemented. + */ + template + void + extract_interpolation_matrices( + const DoFHandler &dof, + dealii::Table<2, FullMatrix> &matrices) + { + if (dof.has_hp_capabilities() == false) + return; + + const dealii::hp::FECollection &fe = + dof.get_fe_collection(); + matrices.reinit(fe.size(), fe.size()); + for (unsigned int i = 0; i < fe.size(); ++i) + for (unsigned int j = 0; j < fe.size(); ++j) + if (i != j) + { + matrices(i, j).reinit(fe[i].n_dofs_per_cell(), + fe[j].n_dofs_per_cell()); + + // see if we can get the interpolation matrices for this + // combination of elements. if not, reset the matrix sizes to zero + // to indicate that this particular combination isn't + // supported. this isn't an outright error right away since we may + // never need to actually interpolate between these two elements + // on actual cells; we simply have to trigger an error if someone + // actually tries + try + { + fe[i].get_interpolation_matrix(fe[j], matrices(i, j)); + } + catch (const typename FiniteElement:: + ExcInterpolationNotImplemented &) + { + matrices(i, j).reinit(0, 0); + } + } + } + + + template + void + restriction_additive(const FiniteElement &, + std::vector> &) + {} + + template + void + restriction_additive( + const dealii::hp::FECollection &fe, + std::vector> &restriction_is_additive) + { + restriction_is_additive.resize(fe.size()); + for (unsigned int f = 0; f < fe.size(); ++f) + { + restriction_is_additive[f].resize(fe[f].n_dofs_per_cell()); + for (unsigned int i = 0; i < fe[f].n_dofs_per_cell(); ++i) + restriction_is_additive[f][i] = fe[f].restriction_is_additive(i); + } + } + } // namespace internal + + + + template + void + SolutionTransfer:: + prepare_for_coarsening_and_refinement(const std::vector &all_in) + { + Assert(prepared_for != pure_refinement, ExcAlreadyPrepForRef()); + Assert(prepared_for != coarsening_and_refinement, + ExcAlreadyPrepForCoarseAndRef()); + + clear(); + n_dofs_old = dof_handler->n_dofs(); + const unsigned int in_size = all_in.size(); + + if constexpr (running_in_debug_mode()) + { + Assert(in_size != 0, + ExcMessage("The array of input vectors you pass to this " + "function has no elements. This is not useful.")); + for (unsigned int i = 0; i < in_size; ++i) + { + Assert(all_in[i].size() == n_dofs_old, + ExcDimensionMismatch(all_in[i].size(), n_dofs_old)); + } + } + + // We need to access dof indices on the entire domain. For + // parallel::shared::Triangulations, ownership of cells might change. If + // they allow artificial cells, we need to restore the "true" cell owners + // temporarily. + // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save + // the current set of subdomain ids, set subdomain ids to the "true" owner + // of each cell upon construction of the TemporarilyRestoreSubdomainIds + // object, and later restore these flags when it is destroyed. + const dealii::internal::parallel::shared:: + TemporarilyRestoreSubdomainIds + subdomain_modifier(dof_handler->get_triangulation()); + + // first count the number + // of cells that will be coarsened + // and that'll stay or be refined + unsigned int n_cells_to_coarsen = 0; + unsigned int n_cells_to_stay_or_refine = 0; + for (const auto &act_cell : dof_handler->active_cell_iterators()) + { + if (act_cell->coarsen_flag_set()) + ++n_cells_to_coarsen; + else + ++n_cells_to_stay_or_refine; + } + Assert((n_cells_to_coarsen + n_cells_to_stay_or_refine) == + dof_handler->get_triangulation().n_active_cells(), + ExcInternalError()); + + unsigned int n_coarsen_fathers = 0; + for (const auto &cell : dof_handler->cell_iterators()) + if (!cell->is_active() && cell->child(0)->coarsen_flag_set()) + ++n_coarsen_fathers; + Assert(n_cells_to_coarsen >= 2 * n_coarsen_fathers, ExcInternalError()); + (void)n_cells_to_coarsen; + + // allocate the needed memory. initialize + // the following arrays in an efficient + // way, without copying much + std::vector>(n_cells_to_stay_or_refine) + .swap(indices_on_cell); + + std::vector>>( + n_coarsen_fathers, + std::vector>(in_size)) + .swap(dof_values_on_cell); + + Table<2, FullMatrix> interpolation_hp; + std::vector> restriction_is_additive; + + internal::extract_interpolation_matrices(*dof_handler, interpolation_hp); + internal::restriction_additive(dof_handler->get_fe_collection(), + restriction_is_additive); + + // we need counters for + // the 'to_stay_or_refine' cells 'n_sr' and + // the 'coarsen_fathers' cells 'n_cf', + unsigned int n_sr = 0, n_cf = 0; + for (const auto &cell : dof_handler->cell_iterators()) + { + // CASE 1: active cell that remains as it is + if (cell->is_active() && !cell->coarsen_flag_set()) + { + const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell(); + indices_on_cell[n_sr].resize(dofs_per_cell); + // cell will not be coarsened, + // so we get away by storing the + // dof indices and later + // interpolating to the children + cell->get_dof_indices(indices_on_cell[n_sr]); + cell_map[std::make_pair(cell->level(), cell->index())] = + Pointerstruct(&indices_on_cell[n_sr], cell->active_fe_index()); + ++n_sr; + } + + // CASE 2: cell is inactive but will become active + else if (cell->has_children() && cell->child(0)->coarsen_flag_set()) + { + // we will need to interpolate from the children of this cell + // to the current one. in the hp-context, this also means + // we need to figure out which finite element space to interpolate + // to since that is not implied by the global FE as in the non-hp- + // case. we choose the 'least dominant fe' on all children from + // the associated FECollection. + std::set fe_indices_children; + for (const auto &child : cell->child_iterators()) + { + Assert(child->is_active() && child->coarsen_flag_set(), + typename dealii::Triangulation< + dim>::ExcInconsistentCoarseningFlags()); + + fe_indices_children.insert(child->active_fe_index()); + } + Assert(!fe_indices_children.empty(), ExcInternalError()); + + const unsigned int target_fe_index = + dof_handler->get_fe_collection().find_dominated_fe_extended( + fe_indices_children, /*codim=*/0); + + Assert(target_fe_index != numbers::invalid_unsigned_int, + dealii::internal::hp::DoFHandlerImplementation:: + ExcNoDominatedFiniteElementOnChildren()); + + const unsigned int dofs_per_cell = + dof_handler->get_fe(target_fe_index).n_dofs_per_cell(); + + std::vector>( + in_size, Vector(dofs_per_cell)) + .swap(dof_values_on_cell[n_cf]); + + + // store the data of each of the input vectors. get this data + // as interpolated onto a finite element space that encompasses + // that of all the children. note that + // cell->get_interpolated_dof_values already does all of the + // interpolations between spaces + for (unsigned int j = 0; j < in_size; ++j) + cell->get_interpolated_dof_values(all_in[j], + dof_values_on_cell[n_cf][j], + target_fe_index); + cell_map[std::make_pair(cell->level(), cell->index())] = + Pointerstruct(&dof_values_on_cell[n_cf], target_fe_index); + ++n_cf; + } + } + Assert(n_sr == n_cells_to_stay_or_refine, ExcInternalError()); + Assert(n_cf == n_coarsen_fathers, ExcInternalError()); + + prepared_for = coarsening_and_refinement; + } + + + + template + void + SolutionTransfer:: + prepare_for_coarsening_and_refinement(const VectorType &in) + { + std::vector all_in(1, in); + prepare_for_coarsening_and_refinement(all_in); + } + + + + template + void + SolutionTransfer::interpolate( + const std::vector &all_in, + std::vector &all_out) const + { + const unsigned int size = all_in.size(); + if constexpr (running_in_debug_mode()) + { + Assert(prepared_for == coarsening_and_refinement, ExcNotPrepared()); + Assert(all_out.size() == size, + ExcDimensionMismatch(all_out.size(), size)); + for (unsigned int i = 0; i < size; ++i) + Assert(all_in[i].size() == n_dofs_old, + ExcDimensionMismatch(all_in[i].size(), n_dofs_old)); + for (unsigned int i = 0; i < all_out.size(); ++i) + Assert(all_out[i].size() == dof_handler->n_dofs(), + ExcDimensionMismatch(all_out[i].size(), + dof_handler->n_dofs())); + for (unsigned int i = 0; i < size; ++i) + for (unsigned int j = 0; j < size; ++j) + Assert(&all_in[i] != &all_out[j], + ExcMessage("Vectors cannot be used as input and output" + " at the same time!")); + } + + // We need to access dof indices on the entire domain. For + // parallel::shared::Triangulations, ownership of cells might change. If + // they allow artificial cells, we need to restore the "true" cell owners + // temporarily. + // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save + // the current set of subdomain ids, set subdomain ids to the "true" owner + // of each cell upon construction of the TemporarilyRestoreSubdomainIds + // object, and later restore these flags when it is destroyed. + const dealii::internal::parallel::shared:: + TemporarilyRestoreSubdomainIds + subdomain_modifier(dof_handler->get_triangulation()); + + Vector local_values; + std::vector dofs; + + typename std::map, + Pointerstruct>::const_iterator pointerstruct, + cell_map_end = cell_map.end(); + + Table<2, FullMatrix> interpolation_hp; + internal::extract_interpolation_matrices(*dof_handler, interpolation_hp); + Vector tmp, tmp2; + + for (const auto &cell : dof_handler->cell_iterators()) + { + pointerstruct = + cell_map.find(std::make_pair(cell->level(), cell->index())); + + if (pointerstruct != cell_map_end) + { + const std::vector *const indexptr = + pointerstruct->second.indices_ptr; + + const std::vector> + *const valuesptr = pointerstruct->second.dof_values_ptr; + + // cell stayed as it was or was refined + if (indexptr != nullptr) + { + Assert(valuesptr == nullptr, ExcInternalError()); + + const unsigned int old_fe_index = + pointerstruct->second.active_fe_index; + + // get the values of each of the input data vectors on this cell + // and prolong it to its children + unsigned int in_size = indexptr->size(); + for (unsigned int j = 0; j < size; ++j) + { + tmp.reinit(in_size, true); + for (unsigned int i = 0; i < in_size; ++i) + tmp(i) = dealii::internal::ElementAccess::get( + all_in[j], (*indexptr)[i]); + + cell->set_dof_values_by_interpolation(tmp, + all_out[j], + old_fe_index, + true); + } + } + else if (valuesptr) + // the children of this cell were deleted + { + Assert(!cell->has_children(), ExcInternalError()); + Assert(indexptr == nullptr, ExcInternalError()); + + const unsigned int dofs_per_cell = + cell->get_fe().n_dofs_per_cell(); + dofs.resize(dofs_per_cell); + // get the local + // indices + cell->get_dof_indices(dofs); + + // distribute the stored data to the new vectors + for (unsigned int j = 0; j < size; ++j) + { + // make sure that the size of the stored indices is the same + // as dofs_per_cell. this is kind of a test if we use the + // same FE in the hp-case. to really do that test we would + // have to store the fe_index of all cells + const Vector *data = + nullptr; + const unsigned int active_fe_index = + cell->active_fe_index(); + if (active_fe_index != + pointerstruct->second.active_fe_index) + { + const unsigned int old_index = + pointerstruct->second.active_fe_index; + const FullMatrix &interpolation_matrix = + interpolation_hp(active_fe_index, old_index); + // The interpolation matrix might be empty when using + // FE_Nothing. + if (interpolation_matrix.empty()) + tmp.reinit(dofs_per_cell, false); + else + { + tmp.reinit(dofs_per_cell, true); + AssertDimension((*valuesptr)[j].size(), + interpolation_matrix.n()); + AssertDimension(tmp.size(), + interpolation_matrix.m()); + interpolation_matrix.vmult(tmp, (*valuesptr)[j]); + } + data = &tmp; + } + else + data = &(*valuesptr)[j]; + + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + dealii::internal::ElementAccess::set( + (*data)(i), dofs[i], all_out[j]); + } + } + // undefined status + else + Assert(false, ExcInternalError()); + } + } + + // We have written into the output vectors. If this was a PETSc vector, for + // example, then we need to compress these to make future operations safe: + for (auto &vec : all_out) + vec.compress(VectorOperation::insert); + } + + + + template + void + SolutionTransfer::interpolate( + const VectorType &in, + VectorType &out) const + { + Assert(in.size() == n_dofs_old, + ExcDimensionMismatch(in.size(), n_dofs_old)); + Assert(out.size() == dof_handler->n_dofs(), + ExcDimensionMismatch(out.size(), dof_handler->n_dofs())); + + std::vector all_in = {in}; + std::vector all_out = {out}; + + interpolate(all_in, all_out); + + out = all_out[0]; + } + + + + template + std::size_t + SolutionTransfer::memory_consumption() const + { + // at the moment we do not include the memory + // consumption of the cell_map as we have no + // real idea about memory consumption of a + // std::map + return (MemoryConsumption::memory_consumption(dof_handler) + + MemoryConsumption::memory_consumption(n_dofs_old) + + sizeof(prepared_for) + + MemoryConsumption::memory_consumption(indices_on_cell) + + MemoryConsumption::memory_consumption(dof_values_on_cell)); + } + + + + template + std::size_t + SolutionTransfer::Pointerstruct:: + memory_consumption() const + { + return sizeof(*this); + } + +} // namespace Legacy + +DEAL_II_NAMESPACE_CLOSE + + +#endif diff --git a/source/numerics/solution_transfer.cc b/source/numerics/solution_transfer.cc index d8c0970e0d..7bb3620d1e 100644 --- a/source/numerics/solution_transfer.cc +++ b/source/numerics/solution_transfer.cc @@ -12,1122 +12,12 @@ // // ------------------------------------------------------------------------ -#include -#include - -#include -#include - -#include -#include -#include - -#include - -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include - -#include - -#include -#include +#include DEAL_II_NAMESPACE_OPEN - -namespace -{ - /** - * Optimized pack function for values assigned on degrees of freedom. - * - * Given that the elements of @p dof_values are stored in consecutive - * locations, we can just memcpy them. Since floating point values don't - * compress well, we also waive the compression that the default - * Utilities::pack() and Utilities::unpack() functions offer. - */ - template - std::vector - pack_dof_values(std::vector> &dof_values, - const unsigned int dofs_per_cell) - { - for (const auto &values : dof_values) - { - AssertDimension(values.size(), dofs_per_cell); - (void)values; - } - - const std::size_t bytes_per_entry = sizeof(value_type) * dofs_per_cell; - - std::vector buffer(dof_values.size() * bytes_per_entry); - for (unsigned int i = 0; i < dof_values.size(); ++i) - std::memcpy(&buffer[i * bytes_per_entry], - &dof_values[i](0), - bytes_per_entry); - - return buffer; - } - - - - /** - * Optimized unpack function for values assigned on degrees of freedom. - */ - template - std::vector> - unpack_dof_values( - const boost::iterator_range::const_iterator> &data_range, - const unsigned int dofs_per_cell) - { - const std::size_t bytes_per_entry = sizeof(value_type) * dofs_per_cell; - const unsigned int n_elements = data_range.size() / bytes_per_entry; - - Assert((data_range.size() % bytes_per_entry == 0), ExcInternalError()); - - std::vector> unpacked_data; - unpacked_data.reserve(n_elements); - for (unsigned int i = 0; i < n_elements; ++i) - { - Vector dof_values(dofs_per_cell); - std::memcpy(&dof_values(0), - &(*std::next(data_range.begin(), i * bytes_per_entry)), - bytes_per_entry); - unpacked_data.emplace_back(std::move(dof_values)); - } - - return unpacked_data; - } -} // namespace - - - -template -SolutionTransfer::SolutionTransfer( - const DoFHandler &dof, - const bool average_values) - : dof_handler(&dof, typeid(*this).name()) - , average_values(average_values) - , handle(numbers::invalid_unsigned_int) -{} - - - -template -void -SolutionTransfer:: - prepare_for_coarsening_and_refinement( - const std::vector &all_in) -{ - const dealii::internal::parallel::shared:: - TemporarilyRestoreSubdomainIds - subdomain_modifier(dof_handler->get_triangulation()); - - for (unsigned int i = 0; i < all_in.size(); ++i) - Assert(all_in[i]->size() == dof_handler->n_dofs(), - ExcDimensionMismatch(all_in[i]->size(), dof_handler->n_dofs())); - - input_vectors = all_in; - register_data_attach(); -} - - - -template -void -SolutionTransfer:: - prepare_for_coarsening_and_refinement(const std::vector &all_in) -{ - std::vector temp(all_in.size()); - - for (std::size_t i = 0; i < temp.size(); ++i) - temp[i] = &(all_in[i]); - - this->prepare_for_coarsening_and_refinement(temp); -} - - - -template -void -SolutionTransfer::register_data_attach() -{ - // TODO: casting away constness is bad - auto tria = const_cast *>( - &dof_handler->get_triangulation()); - Assert(tria != nullptr, ExcInternalError()); - - Assert(handle == numbers::invalid_unsigned_int, - ExcMessage("You can only add one solution per " - "SolutionTransfer object.")); - - handle = tria->register_data_attach( - [this](const typename Triangulation::cell_iterator &cell_, - const CellStatus status) { - return this->pack_callback(cell_, status); - }, - /*returns_variable_size_data=*/dof_handler->has_hp_capabilities()); -} - - - -template -void -SolutionTransfer:: - prepare_for_coarsening_and_refinement(const VectorType &in) -{ - std::vector all_in(1, &in); - prepare_for_coarsening_and_refinement(all_in); -} - - - -template -void -SolutionTransfer::prepare_for_serialization( - const VectorType &in) -{ - std::vector all_in(1, &in); - prepare_for_serialization(all_in); -} - - - -template -void -SolutionTransfer::prepare_for_serialization( - const std::vector &all_in) -{ - prepare_for_coarsening_and_refinement(all_in); -} - - - -template -void -SolutionTransfer::deserialize(VectorType &in) -{ - std::vector all_in(1, &in); - deserialize(all_in); -} - - - -template -void -SolutionTransfer::deserialize( - std::vector &all_in) -{ - register_data_attach(); - - // this makes interpolate() happy - input_vectors.resize(all_in.size()); - - interpolate(all_in); -} - - -template -void -SolutionTransfer::interpolate( - std::vector &all_out) -{ - const dealii::internal::parallel::shared:: - TemporarilyRestoreSubdomainIds - subdomain_modifier(dof_handler->get_triangulation()); - - Assert(input_vectors.size() == all_out.size(), - ExcDimensionMismatch(input_vectors.size(), all_out.size())); - for (unsigned int i = 0; i < all_out.size(); ++i) - Assert(all_out[i]->size() == dof_handler->n_dofs(), - ExcDimensionMismatch(all_out[i]->size(), dof_handler->n_dofs())); - - // TODO: casting away constness is bad - auto tria = const_cast *>( - &dof_handler->get_triangulation()); - Assert(tria != nullptr, ExcInternalError()); - Assert( - handle != numbers::invalid_unsigned_int, - ExcMessage( - "You can only call interpolate() once per SolutionTransfer object.")); - - using Number = typename VectorType::value_type; - - if (average_values) - for (auto *const vec : all_out) - *vec = Number(); - - VectorType valence; - - // initialize valence vector only if we need to average - if (average_values) - valence.reinit(*all_out[0]); - - tria->notify_ready_to_unpack( - handle, - [this, &all_out, &valence]( - const typename Triangulation::cell_iterator &cell_, - const CellStatus status, - const boost::iterator_range::const_iterator> - &data_range) { - this->unpack_callback(cell_, status, data_range, all_out, valence); - }); - - if (average_values) - { - // finalize valence: compress and invert - valence.compress(VectorOperation::add); - for (const auto i : valence.locally_owned_elements()) - valence[i] = (static_cast(valence[i]) == Number() ? - Number() : - (Number(1.0) / static_cast(valence[i]))); - valence.compress(VectorOperation::insert); - - for (auto *const vec : all_out) - { - // compress and weight with valence - vec->compress(VectorOperation::add); - vec->scale(valence); - } - } - else - { - for (auto *const vec : all_out) - vec->compress(VectorOperation::insert); - } - - input_vectors.clear(); - handle = numbers::invalid_unsigned_int; -} - - -template -void -SolutionTransfer::interpolate( - std::vector &all_out) -{ - std::vector temp(all_out.size()); - - for (std::size_t i = 0; i < temp.size(); ++i) - temp[i] = &(all_out[i]); - - this->interpolate(temp); -} - - - -template -void -SolutionTransfer::interpolate(VectorType &out) -{ - std::vector all_out(1, &out); - interpolate(all_out); -} - - - -template -std::vector -SolutionTransfer::pack_callback( - const typename Triangulation::cell_iterator &cell_, - const CellStatus status) -{ - typename DoFHandler::cell_iterator cell(*cell_, dof_handler); - - // create buffer for each individual object - std::vector<::dealii::Vector> dof_values( - input_vectors.size()); - - unsigned int fe_index = 0; - if (dof_handler->has_hp_capabilities()) - { - switch (status) - { - case CellStatus::cell_will_persist: - case CellStatus::cell_will_be_refined: - { - fe_index = cell->future_fe_index(); - break; - } - - case CellStatus::children_will_be_coarsened: - { - // In case of coarsening, we need to find a suitable FE index - // for the parent cell. We choose the 'least dominant fe' - // on all children from the associated FECollection. - if constexpr (running_in_debug_mode()) - { - for (const auto &child : cell->child_iterators()) - Assert(child->is_active() && child->coarsen_flag_set(), - typename dealii::Triangulation< - dim>::ExcInconsistentCoarseningFlags()); - } - - fe_index = dealii::internal::hp::DoFHandlerImplementation:: - dominated_future_fe_on_children(cell); - break; - } - - default: - Assert(false, ExcInternalError()); - break; - } - } - - const unsigned int dofs_per_cell = - dof_handler->get_fe(fe_index).n_dofs_per_cell(); - - if (dofs_per_cell == 0) - return std::vector(); // nothing to do for FE_Nothing - - auto it_input = input_vectors.cbegin(); - auto it_output = dof_values.begin(); - for (; it_input != input_vectors.cend(); ++it_input, ++it_output) - { - it_output->reinit(dofs_per_cell); - cell->get_interpolated_dof_values(*(*it_input), *it_output, fe_index); - } - - return pack_dof_values(dof_values, - dofs_per_cell); -} - - - -template -void -SolutionTransfer::unpack_callback( - const typename Triangulation::cell_iterator &cell_, - const CellStatus status, - const boost::iterator_range::const_iterator> &data_range, - std::vector &all_out, - VectorType &valence) -{ - typename DoFHandler::cell_iterator cell(*cell_, dof_handler); - - unsigned int fe_index = 0; - if (dof_handler->has_hp_capabilities()) - { - switch (status) - { - case CellStatus::cell_will_persist: - case CellStatus::children_will_be_coarsened: - { - fe_index = cell->active_fe_index(); - break; - } - - case CellStatus::cell_will_be_refined: - { - // After refinement, this particular cell is no longer active, - // and its children have inherited its FE index. However, to - // unpack the data on the old cell, we need to recover its FE - // index from one of the children. Just to be sure, we also - // check if all children have the same FE index. - fe_index = cell->child(0)->active_fe_index(); - for (unsigned int child_index = 1; - child_index < cell->n_children(); - ++child_index) - Assert(cell->child(child_index)->active_fe_index() == fe_index, - ExcInternalError()); - break; - } - - default: - Assert(false, ExcInternalError()); - break; - } - } - - const unsigned int dofs_per_cell = - dof_handler->get_fe(fe_index).n_dofs_per_cell(); - - if (dofs_per_cell == 0) - return; // nothing to do for FE_Nothing - - const std::vector<::dealii::Vector> - dof_values = - unpack_dof_values(data_range, - dofs_per_cell); - - // check if sizes match - AssertDimension(dof_values.size(), all_out.size()); - - // check if we have enough dofs provided by the FE object - // to interpolate the transferred data correctly - for (auto it_dof_values = dof_values.begin(); - it_dof_values != dof_values.end(); - ++it_dof_values) - Assert( - dofs_per_cell == it_dof_values->size(), - ExcMessage( - "The transferred data was packed with a different number of dofs than the " - "currently registered FE object assigned to the DoFHandler has.")); - - // distribute data for each registered vector on mesh - auto it_input = dof_values.cbegin(); - auto it_output = all_out.begin(); - for (; it_input != dof_values.cend(); ++it_input, ++it_output) - if (average_values) - cell->distribute_local_to_global_by_interpolation(*it_input, - *(*it_output), - fe_index); - else - cell->set_dof_values_by_interpolation(*it_input, - *(*it_output), - fe_index, - true); - - if (average_values) - { - // compute valence vector if averaging should be performed - Vector ones(dofs_per_cell); - ones = 1.0; - cell->distribute_local_to_global_by_interpolation(ones, - valence, - fe_index); - } -} - - - -template -void -SolutionTransfer::clear() -{ - // nothing to do -} - - - -namespace Legacy -{ - - template - SolutionTransfer::SolutionTransfer( - const DoFHandler &dof) - : dof_handler(&dof, typeid(*this).name()) - , n_dofs_old(0) - , prepared_for(none) - { - Assert((dynamic_cast< - const parallel::distributed::Triangulation *>( - &dof_handler->get_triangulation()) == nullptr), - ExcMessage( - "You are calling the dealii::SolutionTransfer class " - "with a DoFHandler that is built on a " - "parallel::distributed::Triangulation. This will not " - "work for parallel computations. You probably want to " - "use the parallel::distributed::SolutionTransfer class.")); - } - - - - template - SolutionTransfer::~SolutionTransfer() - { - clear(); - } - - - - template - void - SolutionTransfer::clear() - { - indices_on_cell.clear(); - dof_values_on_cell.clear(); - cell_map.clear(); - - prepared_for = none; - } - - - - template - void - SolutionTransfer::prepare_for_pure_refinement() - { - Assert(prepared_for != pure_refinement, ExcAlreadyPrepForRef()); - Assert(prepared_for != coarsening_and_refinement, - ExcAlreadyPrepForCoarseAndRef()); - - clear(); - - // We need to access dof indices on the entire domain. For - // parallel::shared::Triangulations, ownership of cells might change. If - // they allow artificial cells, we need to restore the "true" cell owners - // temporarily. - // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save - // the current set of subdomain ids, set subdomain ids to the "true" owner - // of each cell upon construction of the TemporarilyRestoreSubdomainIds - // object, and later restore these flags when it is destroyed. - const dealii::internal::parallel::shared:: - TemporarilyRestoreSubdomainIds - subdomain_modifier(dof_handler->get_triangulation()); - - const unsigned int n_active_cells = - dof_handler->get_triangulation().n_active_cells(); - n_dofs_old = dof_handler->n_dofs(); - - // efficient reallocation of indices_on_cell - std::vector>(n_active_cells) - .swap(indices_on_cell); - - for (const auto &cell : dof_handler->active_cell_iterators()) - { - const unsigned int i = cell->active_cell_index(); - indices_on_cell[i].resize(cell->get_fe().n_dofs_per_cell()); - // on each cell store the indices of the - // dofs. after refining we get the values - // on the children by taking these - // indices, getting the respective values - // out of the data vectors and prolonging - // them to the children - cell->get_dof_indices(indices_on_cell[i]); - cell_map[std::make_pair(cell->level(), cell->index())] = - Pointerstruct(&indices_on_cell[i], cell->active_fe_index()); - } - prepared_for = pure_refinement; - } - - - - template - void - SolutionTransfer::refine_interpolate( - const VectorType &in, - VectorType &out) const - { - Assert(prepared_for == pure_refinement, ExcNotPrepared()); - Assert(in.size() == n_dofs_old, - ExcDimensionMismatch(in.size(), n_dofs_old)); - Assert(out.size() == dof_handler->n_dofs(), - ExcDimensionMismatch(out.size(), dof_handler->n_dofs())); - Assert(&in != &out, - ExcMessage("Vectors cannot be used as input and output" - " at the same time!")); - - // We need to access dof indices on the entire domain. For - // parallel::shared::Triangulations, ownership of cells might change. If - // they allow artificial cells, we need to restore the "true" cell owners - // temporarily. - // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save - // the current set of subdomain ids, set subdomain ids to the "true" owner - // of each cell upon construction of the TemporarilyRestoreSubdomainIds - // object, and later restore these flags when it is destroyed. - const dealii::internal::parallel::shared:: - TemporarilyRestoreSubdomainIds - subdomain_modifier(dof_handler->get_triangulation()); - - Vector local_values(0); - - typename std::map, - Pointerstruct>::const_iterator pointerstruct, - cell_map_end = cell_map.end(); - - for (const auto &cell : dof_handler->cell_iterators()) - { - pointerstruct = - cell_map.find(std::make_pair(cell->level(), cell->index())); - - if (pointerstruct != cell_map_end) - // this cell was refined or not - // touched at all, so we can get - // the new values by just setting - // or interpolating to the children, - // which is both done by one - // function - { - const unsigned int this_fe_index = - pointerstruct->second.active_fe_index; - const unsigned int dofs_per_cell = - cell->get_dof_handler().get_fe(this_fe_index).n_dofs_per_cell(); - local_values.reinit(dofs_per_cell, true); - - // make sure that the size of the stored indices is the same as - // dofs_per_cell. since we store the desired fe_index, we know - // what this size should be - Assert(dofs_per_cell == (*pointerstruct->second.indices_ptr).size(), - ExcInternalError()); - for (unsigned int i = 0; i < dofs_per_cell; ++i) - local_values(i) = - dealii::internal::ElementAccess::get( - in, (*pointerstruct->second.indices_ptr)[i]); - cell->set_dof_values_by_interpolation(local_values, - out, - this_fe_index, - true); - } - } - } - - - - namespace internal - { - /** - * Generate a table that contains - * interpolation matrices between - * each combination of finite - * elements used in a DoFHandler of - * some kind. Since not all - * elements can be interpolated - * onto each other, the table may - * contain empty matrices for those - * combinations of elements for - * which no such interpolation is - * implemented. - */ - template - void - extract_interpolation_matrices( - const DoFHandler &dof, - dealii::Table<2, FullMatrix> &matrices) - { - if (dof.has_hp_capabilities() == false) - return; - - const dealii::hp::FECollection &fe = - dof.get_fe_collection(); - matrices.reinit(fe.size(), fe.size()); - for (unsigned int i = 0; i < fe.size(); ++i) - for (unsigned int j = 0; j < fe.size(); ++j) - if (i != j) - { - matrices(i, j).reinit(fe[i].n_dofs_per_cell(), - fe[j].n_dofs_per_cell()); - - // see if we can get the interpolation matrices for this - // combination of elements. if not, reset the matrix sizes to zero - // to indicate that this particular combination isn't - // supported. this isn't an outright error right away since we may - // never need to actually interpolate between these two elements - // on actual cells; we simply have to trigger an error if someone - // actually tries - try - { - fe[i].get_interpolation_matrix(fe[j], matrices(i, j)); - } - catch (const typename FiniteElement:: - ExcInterpolationNotImplemented &) - { - matrices(i, j).reinit(0, 0); - } - } - } - - - template - void - restriction_additive(const FiniteElement &, - std::vector> &) - {} - - template - void - restriction_additive( - const dealii::hp::FECollection &fe, - std::vector> &restriction_is_additive) - { - restriction_is_additive.resize(fe.size()); - for (unsigned int f = 0; f < fe.size(); ++f) - { - restriction_is_additive[f].resize(fe[f].n_dofs_per_cell()); - for (unsigned int i = 0; i < fe[f].n_dofs_per_cell(); ++i) - restriction_is_additive[f][i] = fe[f].restriction_is_additive(i); - } - } - } // namespace internal - - - - template - void - SolutionTransfer:: - prepare_for_coarsening_and_refinement(const std::vector &all_in) - { - Assert(prepared_for != pure_refinement, ExcAlreadyPrepForRef()); - Assert(prepared_for != coarsening_and_refinement, - ExcAlreadyPrepForCoarseAndRef()); - - clear(); - n_dofs_old = dof_handler->n_dofs(); - const unsigned int in_size = all_in.size(); - - if constexpr (running_in_debug_mode()) - { - Assert(in_size != 0, - ExcMessage("The array of input vectors you pass to this " - "function has no elements. This is not useful.")); - for (unsigned int i = 0; i < in_size; ++i) - { - Assert(all_in[i].size() == n_dofs_old, - ExcDimensionMismatch(all_in[i].size(), n_dofs_old)); - } - } - - // We need to access dof indices on the entire domain. For - // parallel::shared::Triangulations, ownership of cells might change. If - // they allow artificial cells, we need to restore the "true" cell owners - // temporarily. - // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save - // the current set of subdomain ids, set subdomain ids to the "true" owner - // of each cell upon construction of the TemporarilyRestoreSubdomainIds - // object, and later restore these flags when it is destroyed. - const dealii::internal::parallel::shared:: - TemporarilyRestoreSubdomainIds - subdomain_modifier(dof_handler->get_triangulation()); - - // first count the number - // of cells that will be coarsened - // and that'll stay or be refined - unsigned int n_cells_to_coarsen = 0; - unsigned int n_cells_to_stay_or_refine = 0; - for (const auto &act_cell : dof_handler->active_cell_iterators()) - { - if (act_cell->coarsen_flag_set()) - ++n_cells_to_coarsen; - else - ++n_cells_to_stay_or_refine; - } - Assert((n_cells_to_coarsen + n_cells_to_stay_or_refine) == - dof_handler->get_triangulation().n_active_cells(), - ExcInternalError()); - - unsigned int n_coarsen_fathers = 0; - for (const auto &cell : dof_handler->cell_iterators()) - if (!cell->is_active() && cell->child(0)->coarsen_flag_set()) - ++n_coarsen_fathers; - Assert(n_cells_to_coarsen >= 2 * n_coarsen_fathers, ExcInternalError()); - (void)n_cells_to_coarsen; - - // allocate the needed memory. initialize - // the following arrays in an efficient - // way, without copying much - std::vector>(n_cells_to_stay_or_refine) - .swap(indices_on_cell); - - std::vector>>( - n_coarsen_fathers, - std::vector>(in_size)) - .swap(dof_values_on_cell); - - Table<2, FullMatrix> interpolation_hp; - std::vector> restriction_is_additive; - - internal::extract_interpolation_matrices(*dof_handler, interpolation_hp); - internal::restriction_additive(dof_handler->get_fe_collection(), - restriction_is_additive); - - // we need counters for - // the 'to_stay_or_refine' cells 'n_sr' and - // the 'coarsen_fathers' cells 'n_cf', - unsigned int n_sr = 0, n_cf = 0; - for (const auto &cell : dof_handler->cell_iterators()) - { - // CASE 1: active cell that remains as it is - if (cell->is_active() && !cell->coarsen_flag_set()) - { - const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell(); - indices_on_cell[n_sr].resize(dofs_per_cell); - // cell will not be coarsened, - // so we get away by storing the - // dof indices and later - // interpolating to the children - cell->get_dof_indices(indices_on_cell[n_sr]); - cell_map[std::make_pair(cell->level(), cell->index())] = - Pointerstruct(&indices_on_cell[n_sr], cell->active_fe_index()); - ++n_sr; - } - - // CASE 2: cell is inactive but will become active - else if (cell->has_children() && cell->child(0)->coarsen_flag_set()) - { - // we will need to interpolate from the children of this cell - // to the current one. in the hp-context, this also means - // we need to figure out which finite element space to interpolate - // to since that is not implied by the global FE as in the non-hp- - // case. we choose the 'least dominant fe' on all children from - // the associated FECollection. - std::set fe_indices_children; - for (const auto &child : cell->child_iterators()) - { - Assert(child->is_active() && child->coarsen_flag_set(), - typename dealii::Triangulation< - dim>::ExcInconsistentCoarseningFlags()); - - fe_indices_children.insert(child->active_fe_index()); - } - Assert(!fe_indices_children.empty(), ExcInternalError()); - - const unsigned int target_fe_index = - dof_handler->get_fe_collection().find_dominated_fe_extended( - fe_indices_children, /*codim=*/0); - - Assert(target_fe_index != numbers::invalid_unsigned_int, - dealii::internal::hp::DoFHandlerImplementation:: - ExcNoDominatedFiniteElementOnChildren()); - - const unsigned int dofs_per_cell = - dof_handler->get_fe(target_fe_index).n_dofs_per_cell(); - - std::vector>( - in_size, Vector(dofs_per_cell)) - .swap(dof_values_on_cell[n_cf]); - - - // store the data of each of the input vectors. get this data - // as interpolated onto a finite element space that encompasses - // that of all the children. note that - // cell->get_interpolated_dof_values already does all of the - // interpolations between spaces - for (unsigned int j = 0; j < in_size; ++j) - cell->get_interpolated_dof_values(all_in[j], - dof_values_on_cell[n_cf][j], - target_fe_index); - cell_map[std::make_pair(cell->level(), cell->index())] = - Pointerstruct(&dof_values_on_cell[n_cf], target_fe_index); - ++n_cf; - } - } - Assert(n_sr == n_cells_to_stay_or_refine, ExcInternalError()); - Assert(n_cf == n_coarsen_fathers, ExcInternalError()); - - prepared_for = coarsening_and_refinement; - } - - - - template - void - SolutionTransfer:: - prepare_for_coarsening_and_refinement(const VectorType &in) - { - std::vector all_in(1, in); - prepare_for_coarsening_and_refinement(all_in); - } - - - - template - void - SolutionTransfer::interpolate( - const std::vector &all_in, - std::vector &all_out) const - { - const unsigned int size = all_in.size(); - if constexpr (running_in_debug_mode()) - { - Assert(prepared_for == coarsening_and_refinement, ExcNotPrepared()); - Assert(all_out.size() == size, - ExcDimensionMismatch(all_out.size(), size)); - for (unsigned int i = 0; i < size; ++i) - Assert(all_in[i].size() == n_dofs_old, - ExcDimensionMismatch(all_in[i].size(), n_dofs_old)); - for (unsigned int i = 0; i < all_out.size(); ++i) - Assert(all_out[i].size() == dof_handler->n_dofs(), - ExcDimensionMismatch(all_out[i].size(), - dof_handler->n_dofs())); - for (unsigned int i = 0; i < size; ++i) - for (unsigned int j = 0; j < size; ++j) - Assert(&all_in[i] != &all_out[j], - ExcMessage("Vectors cannot be used as input and output" - " at the same time!")); - } - - // We need to access dof indices on the entire domain. For - // parallel::shared::Triangulations, ownership of cells might change. If - // they allow artificial cells, we need to restore the "true" cell owners - // temporarily. - // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save - // the current set of subdomain ids, set subdomain ids to the "true" owner - // of each cell upon construction of the TemporarilyRestoreSubdomainIds - // object, and later restore these flags when it is destroyed. - const dealii::internal::parallel::shared:: - TemporarilyRestoreSubdomainIds - subdomain_modifier(dof_handler->get_triangulation()); - - Vector local_values; - std::vector dofs; - - typename std::map, - Pointerstruct>::const_iterator pointerstruct, - cell_map_end = cell_map.end(); - - Table<2, FullMatrix> interpolation_hp; - internal::extract_interpolation_matrices(*dof_handler, interpolation_hp); - Vector tmp, tmp2; - - for (const auto &cell : dof_handler->cell_iterators()) - { - pointerstruct = - cell_map.find(std::make_pair(cell->level(), cell->index())); - - if (pointerstruct != cell_map_end) - { - const std::vector *const indexptr = - pointerstruct->second.indices_ptr; - - const std::vector> - *const valuesptr = pointerstruct->second.dof_values_ptr; - - // cell stayed as it was or was refined - if (indexptr != nullptr) - { - Assert(valuesptr == nullptr, ExcInternalError()); - - const unsigned int old_fe_index = - pointerstruct->second.active_fe_index; - - // get the values of each of the input data vectors on this cell - // and prolong it to its children - unsigned int in_size = indexptr->size(); - for (unsigned int j = 0; j < size; ++j) - { - tmp.reinit(in_size, true); - for (unsigned int i = 0; i < in_size; ++i) - tmp(i) = dealii::internal::ElementAccess::get( - all_in[j], (*indexptr)[i]); - - cell->set_dof_values_by_interpolation(tmp, - all_out[j], - old_fe_index, - true); - } - } - else if (valuesptr) - // the children of this cell were deleted - { - Assert(!cell->has_children(), ExcInternalError()); - Assert(indexptr == nullptr, ExcInternalError()); - - const unsigned int dofs_per_cell = - cell->get_fe().n_dofs_per_cell(); - dofs.resize(dofs_per_cell); - // get the local - // indices - cell->get_dof_indices(dofs); - - // distribute the stored data to the new vectors - for (unsigned int j = 0; j < size; ++j) - { - // make sure that the size of the stored indices is the same - // as dofs_per_cell. this is kind of a test if we use the - // same FE in the hp-case. to really do that test we would - // have to store the fe_index of all cells - const Vector *data = - nullptr; - const unsigned int active_fe_index = - cell->active_fe_index(); - if (active_fe_index != - pointerstruct->second.active_fe_index) - { - const unsigned int old_index = - pointerstruct->second.active_fe_index; - const FullMatrix &interpolation_matrix = - interpolation_hp(active_fe_index, old_index); - // The interpolation matrix might be empty when using - // FE_Nothing. - if (interpolation_matrix.empty()) - tmp.reinit(dofs_per_cell, false); - else - { - tmp.reinit(dofs_per_cell, true); - AssertDimension((*valuesptr)[j].size(), - interpolation_matrix.n()); - AssertDimension(tmp.size(), - interpolation_matrix.m()); - interpolation_matrix.vmult(tmp, (*valuesptr)[j]); - } - data = &tmp; - } - else - data = &(*valuesptr)[j]; - - - for (unsigned int i = 0; i < dofs_per_cell; ++i) - dealii::internal::ElementAccess::set( - (*data)(i), dofs[i], all_out[j]); - } - } - // undefined status - else - Assert(false, ExcInternalError()); - } - } - - // We have written into the output vectors. If this was a PETSc vector, for - // example, then we need to compress these to make future operations safe: - for (auto &vec : all_out) - vec.compress(VectorOperation::insert); - } - - - - template - void - SolutionTransfer::interpolate( - const VectorType &in, - VectorType &out) const - { - Assert(in.size() == n_dofs_old, - ExcDimensionMismatch(in.size(), n_dofs_old)); - Assert(out.size() == dof_handler->n_dofs(), - ExcDimensionMismatch(out.size(), dof_handler->n_dofs())); - - std::vector all_in = {in}; - std::vector all_out = {out}; - - interpolate(all_in, all_out); - - out = all_out[0]; - } - - - - template - std::size_t - SolutionTransfer::memory_consumption() const - { - // at the moment we do not include the memory - // consumption of the cell_map as we have no - // real idea about memory consumption of a - // std::map - return (MemoryConsumption::memory_consumption(dof_handler) + - MemoryConsumption::memory_consumption(n_dofs_old) + - sizeof(prepared_for) + - MemoryConsumption::memory_consumption(indices_on_cell) + - MemoryConsumption::memory_consumption(dof_values_on_cell)); - } - - - - template - std::size_t - SolutionTransfer::Pointerstruct:: - memory_consumption() const - { - return sizeof(*this); - } - -} // namespace Legacy - - -/*-------------- Explicit Instantiations -------------------------------*/ #define SPLIT_INSTANTIATIONS_COUNT 4 #ifndef SPLIT_INSTANTIATIONS_INDEX # define SPLIT_INSTANTIATIONS_INDEX 0 -- 2.39.5