// 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<number> tmp(this->get_fe().dofs_per_cell);
- this->get_dof_values(values, tmp);
-
- FullMatrix<double> 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<number> tmp(dofs_per_cell);
+ this->get_dof_values(values, tmp);
+
+ FullMatrix<double> 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
internal::extract_interpolation_matrices(*dof_handler, interpolation_hp);
Vector<typename VectorType::value_type> 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()));
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)
{
}
}
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());
// 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<typename VectorType::value_type> *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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/hp/fe_collection.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/solution_transfer.h>
+
+#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<double> 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<double>, 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<double> 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);
+}
--- /dev/null
+
+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
projects / dealii.git / commitdiff