Merge pull request #8860 from masterleinad/fix_fe_nothing_coarsening

Fix SolutionTransfer when coarsening from FE_Nothing

source/dofs/dof_accessor_get.cc

@@ -69,15 +69,23 @@ DoFCellAccessor<DoFHandlerType, lda>::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<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

source/numerics/solution_transfer.cc

@@ -456,9 +456,7 @@ SolutionTransfer<dim, VectorType, DoFHandlerType>::interpolate(
   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()));
@@ -479,11 +477,8 @@ SolutionTransfer<dim, VectorType, DoFHandlerType>::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<dim, VectorType, DoFHandlerType>::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,33 @@ SolutionTransfer<dim, VectorType, DoFHandlerType>::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<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]);
+                      const FullMatrix<double> &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

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 <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);
+}

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