From 476da6f30178d8eebfffeda5fc3c2e3f9c4ee1c0 Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Thu, 26 Sep 2019 00:01:44 -0400 Subject: [PATCH] Fix SolutionTransfer when corsening from FE_Nothing --- source/dofs/dof_accessor_get.cc | 26 +++++--- source/numerics/solution_transfer.cc | 59 ++++++++--------- tests/hp/solution_transfer_16.cc | 94 ++++++++++++++++++++++++++++ tests/hp/solution_transfer_16.output | 7 +++ 4 files changed, 145 insertions(+), 41 deletions(-) create mode 100644 tests/hp/solution_transfer_16.cc create mode 100644 tests/hp/solution_transfer_16.output diff --git a/source/dofs/dof_accessor_get.cc b/source/dofs/dof_accessor_get.cc index c557fb5fec..8329d7373c 100644 --- a/source/dofs/dof_accessor_get.cc +++ b/source/dofs/dof_accessor_get.cc @@ -69,15 +69,23 @@ DoFCellAccessor::get_interpolated_dof_values( // well, here we need to first get the values from the current // cell and then interpolate it to the element requested. this // can clearly only happen for hp::DoFHandler objects - Vector tmp(this->get_fe().dofs_per_cell); - this->get_dof_values(values, tmp); - - FullMatrix interpolation( - this->dof_handler->get_fe(fe_index).dofs_per_cell, - this->get_fe().dofs_per_cell); - this->dof_handler->get_fe(fe_index).get_interpolation_matrix( - this->get_fe(), interpolation); - interpolation.vmult(interpolated_values, tmp); + const unsigned int dofs_per_cell = this->get_fe().dofs_per_cell; + if (dofs_per_cell == 0) + { + interpolated_values = 0; + } + else + { + Vector tmp(dofs_per_cell); + this->get_dof_values(values, tmp); + + FullMatrix interpolation( + this->dof_handler->get_fe(fe_index).dofs_per_cell, + this->get_fe().dofs_per_cell); + this->dof_handler->get_fe(fe_index).get_interpolation_matrix( + this->get_fe(), interpolation); + interpolation.vmult(interpolated_values, tmp); + } } } else diff --git a/source/numerics/solution_transfer.cc b/source/numerics/solution_transfer.cc index 3456035e60..ede47089dc 100644 --- a/source/numerics/solution_transfer.cc +++ b/source/numerics/solution_transfer.cc @@ -456,9 +456,7 @@ SolutionTransfer::interpolate( internal::extract_interpolation_matrices(*dof_handler, interpolation_hp); Vector tmp, tmp2; - typename DoFHandlerType::cell_iterator cell = dof_handler->begin(), - endc = dof_handler->end(); - for (; cell != endc; ++cell) + for (const auto &cell : dof_handler->cell_iterators()) { pointerstruct = cell_map.find(std::make_pair(cell->level(), cell->index())); @@ -479,11 +477,8 @@ SolutionTransfer::interpolate( 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 + // 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) { @@ -499,8 +494,7 @@ SolutionTransfer::interpolate( } } else if (valuesptr) - // the children of this cell were - // deleted + // the children of this cell were deleted { Assert(!cell->has_children(), ExcInternalError()); Assert(indexptr == nullptr, ExcInternalError()); @@ -511,36 +505,37 @@ SolutionTransfer::interpolate( // indices cell->get_dof_indices(dofs); - // distribute the - // stored data to the - // new vectors + // 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 + // 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; - tmp.reinit(dofs_per_cell, true); - AssertDimension( - (*valuesptr)[j].size(), - interpolation_hp(active_fe_index, old_index).n()); - AssertDimension( - tmp.size(), - interpolation_hp(active_fe_index, old_index).m()); - interpolation_hp(active_fe_index, old_index) - .vmult(tmp, (*valuesptr)[j]); + // The interpolation matrix might be empty when using + // FE_Nothing. + if (interpolation_hp(active_fe_index, old_index).empty()) + tmp.reinit(dofs_per_cell, false); + else + { + tmp.reinit(dofs_per_cell, true); + const unsigned int old_index = + pointerstruct->second.active_fe_index; + AssertDimension( + (*valuesptr)[j].size(), + interpolation_hp(active_fe_index, old_index).n()); + AssertDimension( + tmp.size(), + interpolation_hp(active_fe_index, old_index).m()); + interpolation_hp(active_fe_index, old_index) + .vmult(tmp, (*valuesptr)[j]); + } data = &tmp; } else diff --git a/tests/hp/solution_transfer_16.cc b/tests/hp/solution_transfer_16.cc new file mode 100644 index 0000000000..874ed573c2 --- /dev/null +++ b/tests/hp/solution_transfer_16.cc @@ -0,0 +1,94 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 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. +// +// --------------------------------------------------------------------- + + +// Verify that we can run SolutionTransfer when coarsening a cell that has +// FE_Nothing assigned. + + +#include + +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include + +#include "../tests.h" + +using namespace dealii; + +int +main() +{ + initlog(); + + Triangulation<2> triangulation(Triangulation<2>::none); + GridGenerator::hyper_cube(triangulation); + triangulation.refine_global(1); + + hp::FECollection<2> fe_collection; + fe_collection.push_back(FE_Q<2>(1)); + fe_collection.push_back(FE_Nothing<2>()); + + hp::DoFHandler<2> dof_handler(triangulation); + + // Assign FE_Nothing to the first cell + dof_handler.begin_active()->set_active_fe_index(1); + + dof_handler.distribute_dofs(fe_collection); + deallog << "Initial number of dofs: " << dof_handler.n_dofs() << std::endl; + + // Initialize solution + Vector solution(dof_handler.n_dofs()); + solution = 10; + + deallog << "Initial Vector:" << std::endl; + solution.print(deallog); + + // Coarsen all cells + for (const auto &cell : dof_handler.active_cell_iterators()) + cell->set_coarsen_flag(); + + triangulation.prepare_coarsening_and_refinement(); + + // Interpolate solution + SolutionTransfer<2, Vector, hp::DoFHandler<2>> solution_trans( + dof_handler); + solution_trans.prepare_for_coarsening_and_refinement(solution); + + triangulation.execute_coarsening_and_refinement(); + + // Assign FE_Q(1) to all cells + for (const auto &cell : dof_handler.active_cell_iterators()) + cell->set_active_fe_index(0); + + dof_handler.distribute_dofs(fe_collection); + deallog << "Final number of dofs: " << dof_handler.n_dofs() << std::endl; + + Vector new_solution(dof_handler.n_dofs()); + new_solution = 1.; + solution_trans.interpolate(solution, new_solution); + + deallog << "Vector after solution transfer:" << std::endl; + new_solution.print(deallog); +} diff --git a/tests/hp/solution_transfer_16.output b/tests/hp/solution_transfer_16.output new file mode 100644 index 0000000000..f29e61fb18 --- /dev/null +++ b/tests/hp/solution_transfer_16.output @@ -0,0 +1,7 @@ + +DEAL::Initial number of dofs: 8 +DEAL::Initial Vector: +DEAL::10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 +DEAL::Final number of dofs: 4 +DEAL::Vector after solution transfer: +DEAL::0.00000 10.0000 10.0000 10.0000 -- 2.39.5