]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add dealii comment blob and split extrapolate into an own file
authorDaniel Arndt <d.arndt@math.uni-goettingen.de>
Tue, 27 Sep 2016 15:33:34 +0000 (17:33 +0200)
committerDaniel Arndt <d.arndt@math.uni-goettingen.de>
Mon, 3 Oct 2016 16:20:58 +0000 (18:20 +0200)
include/deal.II/distributed/p4est_wrappers.h
source/distributed/p4est_wrappers.inst.in
source/fe/CMakeLists.txt
source/fe/fe_tools_extrapolate.cc [new file with mode: 0755]
source/fe/fe_tools_extrapolate.inst.in [new file with mode: 0755]
source/fe/fe_tools_interpolate.cc [changed mode: 0644->0755]
source/fe/fe_tools_interpolate.inst.in [changed mode: 0644->0755]
tests/mpi/fe_tools_extrapolate_03.mpirun=10.output [new file with mode: 0644]
tests/mpi/fe_tools_extrapolate_common.h [changed mode: 0644->0755]

index 1bbad8f11427237d4cc5f48c7ba7b00deb7b430e..67181ec79e9f00e49d03f3ae9d175c09153dbbb7 100644 (file)
@@ -1,5 +1,20 @@
-#ifndef P4EST_WRAPPERS_H
-#define P4EST_WRAPPERS_H
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2016 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii__p4est_wrappers_h
+#define dealii__p4est_wrappers_h
 
 #include <deal.II/base/geometry_info.h>
 
@@ -531,4 +546,4 @@ DEAL_II_NAMESPACE_CLOSE
 
 #endif //DEAL_II_WITH_P4EST
 
-#endif //P4EST_WRAPPERS_H
+#endif //dealii__p4est_wrappers_h
index ec5547c9e49e65d62d2f79f6e7ede42c127c2f56..129cc271aa2503f74e044aa513bb975312785960 100644 (file)
@@ -1,6 +1,6 @@
 // ---------------------------------------------------------------------
 //
-// Copyright (C) 2010 - 2014 by the deal.II authors
+// Copyright (C) 2016 - 2016 by the deal.II authors
 //
 // This file is part of the deal.II library.
 //
index e71ac66a0e9caa57b633cd9275a8acbb31c3699f..1c744274b46cee4c117a28aac11ed9095ba3c45c 100644 (file)
@@ -48,6 +48,7 @@ SET(_src
   fe_enriched.cc
   fe_tools.cc
   fe_tools_interpolate.cc
+  fe_tools_extrapolate.cc
   fe_trace.cc
   fe_values.cc
   fe_values_inst2.cc
@@ -91,6 +92,7 @@ SET(_inst
   fe_enriched.inst.in
   fe_tools.inst.in
   fe_tools_interpolate.inst.in
+  fe_tools_extrapolate.inst.in
   fe_trace.inst.in
   fe_values.decl.1.inst.in
   fe_values.decl.2.inst.in
diff --git a/source/fe/fe_tools_extrapolate.cc b/source/fe/fe_tools_extrapolate.cc
new file mode 100755 (executable)
index 0000000..6c82be1
--- /dev/null
@@ -0,0 +1,1567 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2000 - 2016 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/distributed/p4est_wrappers.h>
+#include <deal.II/distributed/tria.h>
+#include <deal.II/fe/fe_tools.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/la_parallel_block_vector.h>
+#include <deal.II/lac/petsc_parallel_vector.h>
+#include <deal.II/lac/petsc_block_vector.h>
+#include <deal.II/lac/petsc_parallel_block_vector.h>
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/lac/trilinos_block_vector.h>
+
+#include <queue>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace FETools
+{
+
+  namespace internal
+  {
+#ifdef DEAL_II_WITH_P4EST
+    // Implementation of the @p extrapolate function
+    // on parallel distributed grids.
+    template <int dim,int spacedim>
+    class ExtrapolateImplementation
+    {
+    public:
+
+      ExtrapolateImplementation ();
+
+      ~ExtrapolateImplementation ();
+
+      template <class InVector,class OutVector>
+      void extrapolate_parallel (const InVector &u2_relevant,
+                                 const DoFHandler<dim,spacedim> &dof2,
+                                 OutVector &u2);
+
+    private:
+
+      // A structure holding all data to
+      // set dofs recursively on cells of arbitrary level
+      struct WorkPackage
+      {
+        const typename dealii::internal::p4est::types<dim>::forest    forest;
+        const typename dealii::internal::p4est::types<dim>::tree      tree;
+        const typename dealii::internal::p4est::types<dim>::locidx    tree_index;
+        const typename DoFHandler<dim,spacedim>::cell_iterator        dealii_cell;
+        const typename dealii::internal::p4est::types<dim>::quadrant  p4est_cell;
+
+        WorkPackage(const typename dealii::internal::p4est::types<dim>::forest    &forest_,
+                    const typename dealii::internal::p4est::types<dim>::tree      &tree_,
+                    const typename dealii::internal::p4est::types<dim>::locidx    &tree_index_,
+                    const typename DoFHandler<dim,spacedim>::cell_iterator        &dealii_cell_,
+                    const typename dealii::internal::p4est::types<dim>::quadrant  &p4est_cell_)
+          :
+          forest(forest_),
+          tree(tree_),
+          tree_index(tree_index_),
+          dealii_cell(dealii_cell_),
+          p4est_cell(p4est_cell_)
+        {}
+      };
+
+
+      // A structure holding all data
+      // of cells needed from other processes
+      // for the extrapolate algorithm.
+      struct CellData
+      {
+        CellData ();
+
+        CellData (const unsigned int dofs_per_cell);
+
+        Vector<double>  dof_values;
+
+        unsigned int tree_index;
+
+        typename dealii::internal::p4est::types<dim>::quadrant  quadrant;
+
+        int receiver;
+
+        bool operator < (const CellData &rhs) const
+        {
+          if (dealii::internal::p4est::functions<dim>::quadrant_compare (&quadrant, &rhs.quadrant) < 0)
+            return true;
+
+          return false;
+        }
+
+        unsigned int bytes_for_buffer () const
+        {
+          return (sizeof(unsigned int) +                                              // dofs_per_cell
+                  dof_values.size() * sizeof(double) +                                // dof_values
+                  sizeof(unsigned int) +                                              // tree_index
+                  sizeof(typename dealii::internal::p4est::types<dim>::quadrant));    // quadrant
+        }
+
+        void pack_data (std::vector<char> &buffer) const
+        {
+          buffer.resize(bytes_for_buffer());
+
+          char *ptr = &buffer[0];
+
+          unsigned int n_dofs = dof_values.size ();
+          std::memcpy(ptr, &n_dofs, sizeof(unsigned int));
+          ptr += sizeof(unsigned int);
+
+          std::memcpy(ptr,dof_values.begin(),n_dofs*sizeof(double));
+          ptr += n_dofs*sizeof(double);
+
+          std::memcpy(ptr,&tree_index,sizeof(unsigned int));
+          ptr += sizeof(unsigned int);
+
+          std::memcpy(ptr,&quadrant,sizeof(typename dealii::internal::p4est::types<dim>::quadrant));
+          ptr += sizeof(typename dealii::internal::p4est::types<dim>::quadrant);
+
+          Assert (ptr == &buffer[0]+buffer.size(),
+                  ExcInternalError());
+        }
+
+        void unpack_data (const std::vector<char> &buffer)
+        {
+          const char *ptr = &buffer[0];
+          unsigned int n_dofs;
+          memcpy(&n_dofs, ptr, sizeof(unsigned int));
+          ptr += sizeof(unsigned int);
+
+          dof_values.reinit(n_dofs);
+          std::memcpy(dof_values.begin(),ptr,n_dofs * sizeof(double));
+          ptr += n_dofs * sizeof(double);
+
+          std::memcpy(&tree_index,ptr,sizeof(unsigned int));
+          ptr += sizeof(unsigned int);
+
+          std::memcpy(&quadrant,ptr,sizeof(typename dealii::internal::p4est::types<dim>::quadrant));
+          ptr += sizeof(typename dealii::internal::p4est::types<dim>::quadrant);
+
+          Assert (ptr == &buffer[0]+buffer.size(),
+                  ExcInternalError());
+        }
+      };
+
+      // Problem: The function extrapolates a polynomial
+      // function from a finer mesh of size $h$ to a polynmial
+      // function of higher degree but on a coarser mesh of
+      // size $2h$. Therefor the mesh has to consist of patches
+      // of four (2d) or eight (3d) cells and the function requires
+      // that the mesh is refined globally at least once.
+      // The algorithm starts on the coarsest level of the grid,
+      // loops over all cells and if a cell has at least one active child,
+      // dof values are set via
+      //   cell->get_interpolated_dof_values(u_input, dof_values)
+      //   cell->set_dof_values_by_interpolation(dof_values, u_output)
+      // both *_interpolation_* functions traverse recursively
+      // over all children down to the active cells
+      // and get/set dof values by interpolation.
+      //
+      // On distributed parallel grids the problem arises, that
+      // if a cell has at least one active child, there is no guarantee
+      // that all children of this cell belong to this process.
+      // There might be children which are owned by other processes
+      // and the algorithm needs to find and has to get the
+      // dof values from these processes and so on.
+      //
+      // Algorithm:
+      // 1) Loop over all coarse cells
+      //      From each coarse cell traverse down the tree
+      //      of refined cells and search for active children
+      //      If there is an active child, check, if all
+      //      other children down to the finest level are part
+      //      of this process, if not, add the cell to the list
+      //      of needs.
+      // 2) Send the list of needs to other processes
+      //      Each process has a list of needs and after
+      //      the loop over all coarse cells (all trees)
+      //      is finished, this list can be send to other processes.
+      // 3) Compute the needs required by other processes
+      //      While computing the values required from other
+      //      processes there can arise new needs and they are
+      //      stored in a list again.
+      // 4) Send the computed values and the list of new needs around
+      //
+      // This procedure has to be repeated until no process needs any
+      // new need and all needs are computed, but there are at most the
+      // "number of grid levels" rounds of sending/receiving cell data.
+      //
+      // Then each process has all data needed for extrapolation.
+
+      // driver function sending/receiving all values from
+      // different processes
+      template <class InVector>
+      void compute_all_non_local_data (const DoFHandler<dim,spacedim> &dof2,
+                                       const InVector                 &u);
+
+      // traverse recursively over
+      // the cells of this tree and
+      // interpolate over patches which
+      // are part of this process
+      template <class InVector,class OutVector>
+      void interpolate_recursively (const typename dealii::internal::p4est::types<dim>::forest    &forest,
+                                    const typename dealii::internal::p4est::types<dim>::tree      &tree,
+                                    const typename dealii::internal::p4est::types<dim>::locidx    &tree_index,
+                                    const typename DoFHandler<dim,spacedim>::cell_iterator        &dealii_cell,
+                                    const typename dealii::internal::p4est::types<dim>::quadrant  &p4est_cell,
+                                    const InVector                                                &u1,
+                                    OutVector                                                     &u2);
+
+      // get dof values for this
+      // cell by interpolation
+      // if a child is reached which
+      // is not part of this process
+      // a new need is created
+      template <class InVector,typename number>
+      void get_interpolated_dof_values (const typename dealii::internal::p4est::types<dim>::forest    &forest,
+                                        const typename dealii::internal::p4est::types<dim>::tree      &tree,
+                                        const typename dealii::internal::p4est::types<dim>::locidx    &tree_index,
+                                        const typename DoFHandler<dim,spacedim>::cell_iterator        &dealii_cell,
+                                        const typename dealii::internal::p4est::types<dim>::quadrant  &p4est_cell,
+                                        const InVector                                                &u,
+                                        Vector<number>                                                &interpolated_values,
+                                        std::vector<CellData>                                         &new_needs);
+
+      // set dof values for this
+      // cell by interpolation
+      template <class OutVector,typename number>
+      void set_dof_values_by_interpolation (const typename DoFHandler<dim,spacedim>::cell_iterator        &dealii_cell,
+                                            const typename dealii::internal::p4est::types<dim>::quadrant  &p4est_cell,
+                                            const Vector<number>                                          &interpolated_values,
+                                            OutVector                                                     &u);
+
+      // compute all cell_data
+      // needed from other processes
+      // to interpolate the part of
+      // this process
+      void compute_needs (const DoFHandler<dim,spacedim> &dof2,
+                          std::vector<CellData>          &new_needs);
+
+      // traverse over the tree
+      // and look for patches this
+      // process has to work on
+      void traverse_tree_recursively (const typename dealii::internal::p4est::types<dim>::forest    &forest,
+                                      const typename dealii::internal::p4est::types<dim>::tree      &tree,
+                                      const typename dealii::internal::p4est::types<dim>::locidx    &tree_index,
+                                      const typename DoFHandler<dim,spacedim>::cell_iterator        &dealii_cell,
+                                      const typename dealii::internal::p4est::types<dim>::quadrant  &p4est_cell,
+                                      std::vector<CellData>                                         &new_needs);
+
+      // traverse recursively
+      // over a patch and look
+      // for cells needed from
+      // other processes for interpolation
+      void traverse_patch_recursively (const typename dealii::internal::p4est::types<dim>::forest    &forest,
+                                       const typename dealii::internal::p4est::types<dim>::tree      &tree,
+                                       const typename dealii::internal::p4est::types<dim>::locidx    &tree_index,
+                                       const typename DoFHandler<dim,spacedim>::cell_iterator        &dealii_cell,
+                                       const typename dealii::internal::p4est::types<dim>::quadrant  &p4est_cell,
+                                       std::vector<CellData>                                         &new_needs);
+
+      // compute dof values of all
+      // cells collected in cells_to_compute
+      // computed cells are deleted
+      // from cells_to_compute and
+      // stored in computed_cells
+      // during computation there can
+      // be new cells needed from
+      // other processes, these cells
+      // are stored in new_needs
+      template <class InVector>
+      void compute_cells (const DoFHandler<dim,spacedim> &dof2,
+                          const InVector                 &u,
+                          std::vector<CellData>          &cells_to_compute,
+                          std::vector<CellData>          &computed_cells,
+                          std::vector<CellData>          &new_needs);
+
+      // traverse over the tree
+      // and compute cells
+      template <class InVector>
+      void compute_cells_in_tree_recursively (const typename dealii::internal::p4est::types<dim>::forest    &forest,
+                                              const typename dealii::internal::p4est::types<dim>::tree      &tree,
+                                              const typename dealii::internal::p4est::types<dim>::locidx    &tree_index,
+                                              const typename DoFHandler<dim,spacedim>::cell_iterator        &dealii_cell,
+                                              const typename dealii::internal::p4est::types<dim>::quadrant  &p4est_cell,
+                                              const InVector                                                &u,
+                                              std::vector<CellData>                                         &cells_to_compute,
+                                              std::vector<CellData>                                         &computed_cells,
+                                              std::vector<CellData>                                         &new_needs);
+
+      // send all cell_data in the vector
+      // cells_to_send to their receivers
+      // and receives a vector of cell_data
+      void send_cells (const std::vector<CellData>  &cells_to_send,
+                       std::vector<CellData>  &received_cells) const;
+
+      // add new cell_data to
+      // the ordered list new_needs
+      // uses cell_data_insert
+      void add_new_need (const typename dealii::internal::p4est::types<dim>::forest    &forest,
+                         const typename dealii::internal::p4est::types<dim>::locidx    &tree_index,
+                         const typename DoFHandler<dim,spacedim>::cell_iterator        &dealii_cell,
+                         const typename dealii::internal::p4est::types<dim>::quadrant  &p4est_cell,
+                         std::vector<CellData>                                         &new_needs) const;
+
+      // binary search in cells_list
+      // assume that cells_list
+      // is ordered
+      static int cell_data_search (const CellData               &cell_data,
+                                   const std::vector<CellData>  &cells_list);
+
+      // insert cell_data into a sorted
+      // vector cells_list at the correct
+      // position if cell_data
+      // not exists already in cells_list
+      static void cell_data_insert (const CellData         &cell_data,
+                                    std::vector<CellData>  &cells_list);
+
+      MPI_Comm  communicator;
+
+      // a vector of all cells this process has
+      // computed or received data
+      std::vector<CellData>  available_cells;
+
+      // stores the indices of dofs on more refined ghosted cells along
+      // with the maximum level
+      std::map<types::global_dof_index, int> dofs_on_refined_neighbors;
+    };
+
+    template <>
+    class ExtrapolateImplementation<1,1>
+    {
+    public:
+
+      ExtrapolateImplementation ()
+      {
+        AssertThrow(false, ExcNotImplemented())
+      }
+
+      ~ExtrapolateImplementation ()
+      {}
+
+      template <class InVector,class OutVector>
+      void extrapolate_parallel (const InVector &/*u2_relevant*/,
+                                 const DoFHandler<1,1> &/*dof2*/,
+                                 OutVector &/*u2*/)
+      {}
+    };
+
+    template <>
+    class ExtrapolateImplementation<1,2>
+    {
+    public:
+
+      ExtrapolateImplementation ()
+      {
+        AssertThrow(false, ExcNotImplemented())
+      }
+
+      ~ExtrapolateImplementation ()
+      {}
+
+      template <class InVector,class OutVector>
+      void extrapolate_parallel (const InVector &/*u2_relevant*/,
+                                 const DoFHandler<1,2> &/*dof2*/,
+                                 OutVector &/*u2*/)
+      {}
+    };
+
+    template <>
+    class ExtrapolateImplementation<1,3>
+    {
+    public:
+
+      ExtrapolateImplementation ()
+      {
+        AssertThrow(false, ExcNotImplemented())
+      }
+
+      ~ExtrapolateImplementation ()
+      {}
+
+      template <class InVector,class OutVector>
+      void extrapolate_parallel (const InVector &/*u2_relevant*/,
+                                 const DoFHandler<1,3> &/*dof2*/,
+                                 OutVector &/*u2*/)
+      {}
+    };
+
+
+
+    template <int dim,int spacedim>
+    ExtrapolateImplementation<dim,spacedim>::
+    ExtrapolateImplementation ()
+    {}
+
+
+
+    template <int dim,int spacedim>
+    ExtrapolateImplementation<dim,spacedim>::
+    ~ExtrapolateImplementation ()
+    {}
+
+
+
+    template <int dim,int spacedim>
+    ExtrapolateImplementation<dim,spacedim>::
+    CellData::CellData ()
+      : tree_index (0),
+        receiver (0)
+    {}
+
+
+
+    template <int dim,int spacedim>
+    ExtrapolateImplementation<dim,spacedim>::
+    CellData::CellData (const unsigned int dofs_per_cell)
+      : tree_index (0),
+        receiver (0)
+    {
+      dof_values.reinit (dofs_per_cell);
+    }
+
+
+
+    template <int dim,int spacedim>
+    template <class InVector,class OutVector>
+    void
+    ExtrapolateImplementation<dim,spacedim>::
+    interpolate_recursively (const typename dealii::internal::p4est::types<dim>::forest    &forest,
+                             const typename dealii::internal::p4est::types<dim>::tree      &tree,
+                             const typename dealii::internal::p4est::types<dim>::locidx    &tree_index,
+                             const typename DoFHandler<dim,spacedim>::cell_iterator        &dealii_cell,
+                             const typename dealii::internal::p4est::types<dim>::quadrant  &p4est_cell,
+                             const InVector                                                &u1,
+                             OutVector                                                     &u2)
+    {
+      // check if this cell exists in the local p4est
+      int idx = sc_array_bsearch(const_cast<sc_array_t *>(&tree.quadrants),
+                                 &p4est_cell,
+                                 dealii::internal::p4est::functions<dim>::quadrant_compare);
+
+      // if neither this cell nor one of it's children belongs to us, don't do anything
+      if (idx == -1 && (dealii::internal::p4est::functions<dim>::
+                        quadrant_overlaps_tree (const_cast<typename dealii::internal::p4est::types<dim>::tree *>(&tree),
+                                                &p4est_cell)
+                        == false))
+        return;
+
+      bool p4est_has_children = (idx == -1);
+
+      bool locally_owned_children = false;
+      if (p4est_has_children)
+        {
+          Assert (dealii_cell->has_children (), ExcInternalError ());
+
+          // check if at least one child is locally owned on our process
+          for (unsigned int child_n=0; child_n<dealii_cell->n_children(); ++child_n)
+            if (dealii_cell->child(child_n)->active())
+              if (dealii_cell->child(child_n)->is_locally_owned())
+                {
+                  locally_owned_children=true;
+                  break;
+                }
+        }
+
+      if (locally_owned_children)
+        {
+          const FiniteElement<dim,spacedim> &fe            = dealii_cell->get_dof_handler().get_fe();
+          const unsigned int                 dofs_per_cell = fe.dofs_per_cell;
+
+          Vector<typename OutVector::value_type> interpolated_values(dofs_per_cell);
+
+          std::vector<CellData>  new_needs;
+          get_interpolated_dof_values (forest,
+                                       tree,
+                                       tree_index,
+                                       dealii_cell,
+                                       p4est_cell,
+                                       u1,
+                                       interpolated_values,
+                                       new_needs);
+
+          // at this point of
+          // the procedure no new
+          // needs should come up
+          Assert (new_needs.size () == 0,
+                  ExcInternalError ());
+
+          set_dof_values_by_interpolation (dealii_cell,
+                                           p4est_cell,
+                                           interpolated_values,
+                                           u2);
+        }
+
+
+
+    }
+
+
+
+    template <int dim,int spacedim>
+    template <class InVector,typename number>
+    void
+    ExtrapolateImplementation<dim,spacedim>::
+    get_interpolated_dof_values (const typename dealii::internal::p4est::types<dim>::forest    &forest,
+                                 const typename dealii::internal::p4est::types<dim>::tree      &tree,
+                                 const typename dealii::internal::p4est::types<dim>::locidx    &tree_index,
+                                 const typename DoFHandler<dim,spacedim>::cell_iterator        &dealii_cell,
+                                 const typename dealii::internal::p4est::types<dim>::quadrant  &p4est_cell,
+                                 const InVector                                                &u,
+                                 Vector<number>                                                &interpolated_values,
+                                 std::vector<CellData>                                         &new_needs)
+    {
+      if (!dealii_cell->has_children ())
+        {
+          if (dealii_cell->is_locally_owned ())
+            {
+              // if this is one of our cells,
+              // get dof values from input vector
+              dealii_cell->get_dof_values (u, interpolated_values);
+            }
+          else
+            {
+              add_new_need (forest,
+                            tree_index,
+                            dealii_cell,
+                            p4est_cell,
+                            new_needs);
+            }
+        }
+      else
+        {
+          const FiniteElement<dim,spacedim> &fe            = dealii_cell->get_dof_handler().get_fe();
+          const unsigned int                 dofs_per_cell = fe.dofs_per_cell;
+
+          Assert (&fe != 0,
+                  ExcNotInitialized());
+          Assert (interpolated_values.size() == dofs_per_cell,
+                  ExcDimensionMismatch(interpolated_values.size(), dofs_per_cell));
+          Assert (u.size() == dealii_cell->get_dof_handler().n_dofs(),
+                  ExcDimensionMismatch(u.size(), dealii_cell->get_dof_handler().n_dofs()));
+
+          Vector<number> tmp1(dofs_per_cell);
+          Vector<number> tmp2(dofs_per_cell);
+
+          interpolated_values = 0;
+          std::vector<bool> restriction_is_additive (dofs_per_cell);
+          for (unsigned int i=0; i<dofs_per_cell; ++i)
+            restriction_is_additive[i] = fe.restriction_is_additive(i);
+
+          typename dealii::internal::p4est::types<dim>::quadrant
+          p4est_child[GeometryInfo<dim>::max_children_per_cell];
+
+          dealii::internal::p4est::init_quadrant_children<dim> (p4est_cell,
+                                                                p4est_child);
+
+          bool found_child = true;
+          for (unsigned int c=0; c<GeometryInfo<dim>::max_children_per_cell; ++c)
+            {
+              if (dealii::internal::p4est::functions<dim>::
+                  quadrant_overlaps_tree (const_cast<typename dealii::internal::p4est::types<dim>::tree *>(&tree),
+                                          &p4est_child[c])
+                  == false)
+                {
+                  // this is a cell this process needs
+                  // data from another process
+
+                  // check if this cell is
+                  // available from other
+                  // computed patches
+                  CellData  cell_data;
+                  cell_data.quadrant = p4est_child[c];
+                  int pos = cell_data_search (cell_data, available_cells);
+
+                  if (pos == -1)
+                    {
+                      // data is not available
+                      // create a new need
+                      found_child = false;
+
+                      add_new_need (forest,
+                                    tree_index,
+                                    dealii_cell->child(c),
+                                    p4est_child[c],
+                                    new_needs);
+                    }
+                  else
+                    {
+                      Assert (available_cells[pos].dof_values.size() == dofs_per_cell,
+                              ExcDimensionMismatch(available_cells[pos].dof_values.size(), dofs_per_cell));
+
+                      tmp1 = available_cells[pos].dof_values;
+                    }
+                }
+              else
+                {
+                  // get the values from the present child, if necessary by
+                  // interpolation either from its own children or
+                  // by interpolating from the finite element on an active
+                  // child to the finite element space requested here
+                  get_interpolated_dof_values (forest,
+                                               tree,
+                                               tree_index,
+                                               dealii_cell->child(c),
+                                               p4est_child[c],
+                                               u,
+                                               tmp1,
+                                               new_needs);
+                }
+
+              if (found_child)
+                {
+                  // interpolate these to the mother cell
+                  fe.get_restriction_matrix(c, dealii_cell->refinement_case()).vmult (tmp2, tmp1);
+
+                  // and add up or set them in the output vector
+                  for (unsigned int i=0; i<dofs_per_cell; ++i)
+                    if (tmp2(i) != number())
+                      {
+                        if (restriction_is_additive[i])
+                          interpolated_values(i) += tmp2(i);
+                        else
+                          interpolated_values(i) = tmp2(i);
+                      }
+                }
+            }
+
+          if (found_child == false)
+            interpolated_values = 0;
+        }
+    }
+
+
+
+    template <int dim,int spacedim>
+    template <class OutVector,typename number>
+    void
+    ExtrapolateImplementation<dim,spacedim>::
+    set_dof_values_by_interpolation (const typename DoFHandler<dim,spacedim>::cell_iterator        &dealii_cell,
+                                     const typename dealii::internal::p4est::types<dim>::quadrant  &p4est_cell,
+                                     const Vector<number>                                          &local_values,
+                                     OutVector                                                     &u)
+    {
+      const FiniteElement<dim,spacedim> &fe            = dealii_cell->get_dof_handler().get_fe();
+      const unsigned int                 dofs_per_cell = fe.dofs_per_cell;
+
+      if (!dealii_cell->has_children ())
+        {
+          if (dealii_cell->is_locally_owned ())
+            {
+              std::vector<types::global_dof_index> indices(dofs_per_cell);
+              dealii_cell->get_dof_indices(indices);
+              for (unsigned int j=0; j<dofs_per_cell; ++j)
+                {
+                  // don't set this dof if there is a more refined ghosted neighbor setting this dof.
+                  const bool on_refined_neighbor
+                    = (dofs_on_refined_neighbors.find(indices[j])!=dofs_on_refined_neighbors.end());
+                  if (!(on_refined_neighbor && dofs_on_refined_neighbors[indices[j]]>dealii_cell->level()))
+                    u(indices[j]) = local_values(j);
+                }
+            }
+        }
+      else
+        {
+          Assert (&fe != 0,
+                  ExcNotInitialized());
+          Assert (local_values.size() == dofs_per_cell,
+                  ExcDimensionMismatch(local_values.size(), dofs_per_cell));
+          Assert (u.size() == dealii_cell->get_dof_handler().n_dofs(),
+                  ExcDimensionMismatch(u.size(), dealii_cell->get_dof_handler().n_dofs()));
+
+          Vector<number> tmp(dofs_per_cell);
+
+          typename dealii::internal::p4est::types<dim>::quadrant
+          p4est_child[GeometryInfo<dim>::max_children_per_cell];
+
+          dealii::internal::p4est::init_quadrant_children<dim> (p4est_cell,
+                                                                p4est_child);
+
+          for (unsigned int c=0; c<GeometryInfo<dim>::max_children_per_cell; ++c)
+            {
+              if (tmp.size() > 0)
+                fe.get_prolongation_matrix(c, dealii_cell->refinement_case())
+                .vmult (tmp, local_values);
+
+              set_dof_values_by_interpolation (dealii_cell->child(c),
+                                               p4est_child[c],
+                                               tmp,
+                                               u);
+            }
+        }
+    }
+
+
+
+    template <int dim,int spacedim>
+    void
+    ExtrapolateImplementation<dim,spacedim>::
+    compute_needs (const DoFHandler<dim,spacedim> &dof2,
+                   std::vector<CellData>          &new_needs)
+    {
+      const parallel::distributed::Triangulation< dim, spacedim > *tr
+        = (dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>
+           (&dof2.get_tria()));
+
+      Assert (tr != 0, ExcInternalError());
+
+      typename DoFHandler<dim,spacedim>::cell_iterator
+      cell=dof2.begin(0),
+      endc=dof2.end(0);
+      for (; cell!=endc; ++cell)
+        {
+          if (dealii::internal::p4est::tree_exists_locally<dim> (tr->parallel_forest,
+                                                                 tr->coarse_cell_to_p4est_tree_permutation[cell->index()])
+              == false)
+            continue;
+
+          typename dealii::internal::p4est::types<dim>::quadrant p4est_coarse_cell;
+          const unsigned int tree_index = tr->coarse_cell_to_p4est_tree_permutation[cell->index()];
+          typename dealii::internal::p4est::types<dim>::tree *tree = tr->init_tree(cell->index());
+
+          dealii::internal::p4est::init_coarse_quadrant<dim>(p4est_coarse_cell);
+
+          // make sure that each cell on the
+          // coarsest level is at least once
+          // refined, otherwise, these cells
+          // can't be treated and would
+          // generate a bogus result
+          {
+            int idx = sc_array_bsearch(const_cast<sc_array_t *>(&tree->quadrants),
+                                       &p4est_coarse_cell,
+                                       dealii::internal::p4est::functions<dim>::quadrant_compare);
+
+            AssertThrow (idx == -1, ExcGridNotRefinedAtLeastOnce ());
+          }
+
+          traverse_tree_recursively (*tr->parallel_forest,
+                                     *tree,
+                                     tree_index,
+                                     cell,
+                                     p4est_coarse_cell,
+                                     new_needs);
+        }
+    }
+
+
+
+    template <int dim,int spacedim>
+    void
+    ExtrapolateImplementation<dim,spacedim>::
+    traverse_tree_recursively (const typename dealii::internal::p4est::types<dim>::forest    &forest,
+                               const typename dealii::internal::p4est::types<dim>::tree      &tree,
+                               const typename dealii::internal::p4est::types<dim>::locidx    &tree_index,
+                               const typename DoFHandler<dim,spacedim>::cell_iterator        &dealii_cell,
+                               const typename dealii::internal::p4est::types<dim>::quadrant  &p4est_cell,
+                               std::vector<CellData>                                         &new_needs)
+    {
+      // check if this cell exists in the local p4est
+      int idx = sc_array_bsearch(const_cast<sc_array_t *>(&tree.quadrants),
+                                 &p4est_cell,
+                                 dealii::internal::p4est::functions<dim>::quadrant_compare);
+
+      // if neither this cell nor one of it's children belongs to us, don't do anything
+      if (idx == -1 && (dealii::internal::p4est::functions<dim>::
+                        quadrant_overlaps_tree (const_cast<typename dealii::internal::p4est::types<dim>::tree *>(&tree),
+                                                &p4est_cell)
+                        == false))
+        return;
+
+      bool p4est_has_children = (idx == -1);
+
+      // this cell is part of a patch
+      // this process has to interpolate on
+      // if there is at least one locally
+      // owned child
+      bool locally_owned_children = false;
+      if (p4est_has_children)
+        {
+          Assert (dealii_cell->has_children (), ExcInternalError ());
+
+          for (unsigned int child_n=0; child_n<dealii_cell->n_children(); ++child_n)
+            if (dealii_cell->child(child_n)->active())
+              if (dealii_cell->child(child_n)->is_locally_owned())
+                {
+                  locally_owned_children=true;
+                  break;
+                }
+        }
+
+      // traverse the patch recursively
+      // to find cells needed from
+      // other processes
+      if (locally_owned_children)
+        {
+          traverse_patch_recursively (forest,
+                                      tree,
+                                      tree_index,
+                                      dealii_cell,
+                                      p4est_cell,
+                                      new_needs);
+        }
+
+      // traverse recursively over this tree
+      if (p4est_has_children)
+        {
+          typename dealii::internal::p4est::types<dim>::quadrant
+          p4est_child[GeometryInfo<dim>::max_children_per_cell];
+
+          dealii::internal::p4est::init_quadrant_children<dim> (p4est_cell,
+                                                                p4est_child);
+
+          for (unsigned int c=0; c<GeometryInfo<dim>::max_children_per_cell; ++c)
+            {
+              traverse_tree_recursively (forest,
+                                         tree,
+                                         tree_index,
+                                         dealii_cell->child(c),
+                                         p4est_child[c],
+                                         new_needs);
+            }
+        }
+    }
+
+
+
+    template <int dim,int spacedim>
+    void
+    ExtrapolateImplementation<dim,spacedim>::
+    traverse_patch_recursively (const typename dealii::internal::p4est::types<dim>::forest    &forest,
+                                const typename dealii::internal::p4est::types<dim>::tree      &tree,
+                                const typename dealii::internal::p4est::types<dim>::locidx    &tree_index,
+                                const typename DoFHandler<dim,spacedim>::cell_iterator        &dealii_cell,
+                                const typename dealii::internal::p4est::types<dim>::quadrant  &p4est_cell,
+                                std::vector<CellData>                                         &new_needs)
+    {
+      if (dealii_cell->has_children ())
+        {
+          Assert (dealii_cell->has_children (), ExcInternalError ());
+
+          typename dealii::internal::p4est::types<dim>::quadrant
+          p4est_child[GeometryInfo<dim>::max_children_per_cell];
+
+          dealii::internal::p4est::init_quadrant_children<dim> (p4est_cell,
+                                                                p4est_child);
+
+          for (unsigned int c=0; c<GeometryInfo<dim>::max_children_per_cell; ++c)
+            {
+              // check if this child
+              // is locally available
+              // in the p4est
+              if (dealii::internal::p4est::functions<dim>::
+                  quadrant_overlaps_tree (const_cast<typename dealii::internal::p4est::types<dim>::tree *>(&tree),
+                                          &p4est_child[c])
+                  == false)
+                {
+                  // this is a cell for which this process
+                  // needs data from another process
+                  // so add the cell to the list
+                  add_new_need (forest,
+                                tree_index,
+                                dealii_cell->child(c),
+                                p4est_child[c],
+                                new_needs);
+                }
+              else
+                {
+                  // at least some part of
+                  // the tree rooted in this
+                  // child is locally available
+                  traverse_patch_recursively (forest,
+                                              tree,
+                                              tree_index,
+                                              dealii_cell->child(c),
+                                              p4est_child[c],
+                                              new_needs);
+                }
+            }
+        }
+    }
+
+
+
+
+    template <int dim,int spacedim>
+    template <class InVector>
+    void
+    ExtrapolateImplementation<dim,spacedim>::
+    compute_cells (const DoFHandler<dim,spacedim> &dof2,
+                   const InVector                 &u,
+                   std::vector<CellData>          &cells_to_compute,
+                   std::vector<CellData>          &computed_cells,
+                   std::vector<CellData>          &new_needs)
+    {
+      const parallel::distributed::Triangulation< dim, spacedim > *tr
+        = (dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>
+           (&dof2.get_tria()));
+
+      Assert (tr != 0, ExcInternalError());
+
+      // collect in a set all trees this
+      // process has to compute cells on
+      std::set<unsigned int>  trees;
+      for (typename std::vector<CellData>::const_iterator it=cells_to_compute.begin(); it!=cells_to_compute.end(); ++it)
+        trees.insert (it->tree_index);
+
+      typename DoFHandler<dim,spacedim>::cell_iterator
+      cell=dof2.begin(0),
+      endc=dof2.end(0);
+      for (; cell!=endc; ++cell)
+        {
+          // check if this is a tree this process has to
+          // work on and that this tree is in the p4est
+          const unsigned int tree_index = tr->coarse_cell_to_p4est_tree_permutation[cell->index()];
+
+          if ((trees.find (tree_index) == trees.end ())
+              ||
+              (dealii::internal::p4est::tree_exists_locally<dim> (tr->parallel_forest,
+                                                                  tree_index)
+               == false))
+            continue;
+
+          typename dealii::internal::p4est::types<dim>::quadrant p4est_coarse_cell;
+          typename dealii::internal::p4est::types<dim>::tree *tree = tr->init_tree(cell->index());
+
+          dealii::internal::p4est::init_coarse_quadrant<dim>(p4est_coarse_cell);
+
+          compute_cells_in_tree_recursively (*tr->parallel_forest,
+                                             *tree,
+                                             tree_index,
+                                             cell,
+                                             p4est_coarse_cell,
+                                             u,
+                                             cells_to_compute,
+                                             computed_cells,
+                                             new_needs);
+        }
+    }
+
+
+
+    template <int dim,int spacedim>
+    template <class InVector>
+    void
+    ExtrapolateImplementation<dim,spacedim>::
+    compute_cells_in_tree_recursively (const typename dealii::internal::p4est::types<dim>::forest    &forest,
+                                       const typename dealii::internal::p4est::types<dim>::tree      &tree,
+                                       const typename dealii::internal::p4est::types<dim>::locidx    &tree_index,
+                                       const typename DoFHandler<dim,spacedim>::cell_iterator        &dealii_cell,
+                                       const typename dealii::internal::p4est::types<dim>::quadrant  &p4est_cell,
+                                       const InVector                                                &u,
+                                       std::vector<CellData>                                         &cells_to_compute,
+                                       std::vector<CellData>                                         &computed_cells,
+                                       std::vector<CellData>                                         &new_needs)
+    {
+      if (cells_to_compute.size () == 0)
+        return;
+
+      // check if this cell exists in the local p4est
+      int idx = sc_array_bsearch(const_cast<sc_array_t *>(&tree.quadrants),
+                                 &p4est_cell,
+                                 dealii::internal::p4est::functions<dim>::quadrant_compare);
+
+      // if neither this cell nor one one of it's children belongs to us, don't do anything
+      if (idx == -1 && (dealii::internal::p4est::functions<dim>::
+                        quadrant_overlaps_tree (const_cast<typename dealii::internal::p4est::types<dim>::tree *>(&tree),
+                                                &p4est_cell)
+                        == false))
+        return;
+
+      bool p4est_has_children = (idx == -1);
+
+      // check if this quadrant is in the list
+      CellData  cell_data;
+      cell_data.quadrant = p4est_cell;
+      int pos = cell_data_search (cell_data, cells_to_compute);
+      if (pos != -1)
+        {
+          std::vector<CellData>  tmp;
+          // compute dof values
+          get_interpolated_dof_values (forest,
+                                       tree,
+                                       tree_index,
+                                       dealii_cell,
+                                       p4est_cell,
+                                       u,
+                                       cells_to_compute[pos].dof_values,
+                                       tmp);
+          // there is no new cell_data
+          // this process needs for computing this cell,
+          // store cell_data in the list of
+          // computed cells and erase this cell
+          // from the list of cells to compute
+          if (tmp.size () == 0)
+            {
+              cell_data_insert (cells_to_compute[pos], computed_cells);
+              cells_to_compute.erase (cells_to_compute.begin () + pos);
+            }
+          else
+            {
+              for (unsigned int i=0; i<tmp.size (); ++i)
+                cell_data_insert (tmp[i], new_needs);
+            }
+        }
+
+      // search recursively over this tree
+      if (p4est_has_children)
+        {
+          typename dealii::internal::p4est::types<dim>::quadrant
+          p4est_child[GeometryInfo<dim>::max_children_per_cell];
+
+          dealii::internal::p4est::init_quadrant_children<dim> (p4est_cell,
+                                                                p4est_child);
+
+          for (unsigned int c=0; c<GeometryInfo<dim>::max_children_per_cell; ++c)
+            {
+              compute_cells_in_tree_recursively (forest,
+                                                 tree,
+                                                 tree_index,
+                                                 dealii_cell->child(c),
+                                                 p4est_child[c],
+                                                 u,
+                                                 cells_to_compute,
+                                                 computed_cells,
+                                                 new_needs);
+            }
+        }
+    }
+
+
+
+    template <int dim,int spacedim>
+    void
+    ExtrapolateImplementation<dim,spacedim>::
+    send_cells (const std::vector<CellData>  &cells_to_send,
+                std::vector<CellData>        &received_cells) const
+    {
+      std::vector<std::vector<char> > sendbuffers (cells_to_send.size());
+      std::vector<std::vector<char> >::iterator buffer = sendbuffers.begin();
+      std::vector<MPI_Request> requests (cells_to_send.size());
+      std::vector<unsigned int> destinations;
+
+      // send data
+      unsigned int idx=0;
+      for (typename std::vector<CellData>::const_iterator it=cells_to_send.begin();
+           it!=cells_to_send.end();
+           ++it, ++idx)
+        {
+          destinations.push_back (it->receiver);
+
+          it->pack_data (*buffer);
+          MPI_Isend (&(*buffer)[0], buffer->size(),
+                     MPI_BYTE,
+                     it->receiver,
+                     123,
+                     communicator,
+                     &requests[idx]);
+        }
+
+      Assert(destinations.size()==cells_to_send.size(), ExcInternalError());
+
+      std::vector<unsigned int> senders
+        = Utilities::MPI::compute_point_to_point_communication_pattern(communicator, destinations);
+
+      // receive data
+      std::vector<char> receive;
+      CellData  cell_data;
+      for (unsigned int index=0; index<senders.size (); ++index)
+        {
+          MPI_Status status;
+          int len;
+          MPI_Probe(MPI_ANY_SOURCE, 123, communicator, &status);
+          MPI_Get_count(&status, MPI_BYTE, &len);
+          receive.resize (len);
+
+          char *buf = &receive[0];
+          MPI_Recv (buf, len, MPI_BYTE, status.MPI_SOURCE, status.MPI_TAG, communicator, &status);
+
+          cell_data.unpack_data (receive);
+
+          // this process has to send this
+          // cell back to the sender
+          // the receiver is the old sender
+          cell_data.receiver = status.MPI_SOURCE;
+
+          received_cells.push_back (cell_data);
+        }
+
+      if (requests.size () > 0)
+        MPI_Waitall(requests.size(), &requests[0], MPI_STATUSES_IGNORE);
+
+      // finally sort the list of cells
+      std::sort (received_cells.begin (), received_cells.end ());
+    }
+
+
+
+    template <int dim,int spacedim>
+    void
+    ExtrapolateImplementation<dim,spacedim>::
+    add_new_need (const typename dealii::internal::p4est::types<dim>::forest    &forest,
+                  const typename dealii::internal::p4est::types<dim>::locidx    &tree_index,
+                  const typename DoFHandler<dim,spacedim>::cell_iterator        &dealii_cell,
+                  const typename dealii::internal::p4est::types<dim>::quadrant  &p4est_cell,
+                  std::vector<CellData>                                         &new_needs) const
+    {
+      const FiniteElement<dim,spacedim> &fe            = dealii_cell->get_dof_handler().get_fe();
+      const unsigned int                 dofs_per_cell = fe.dofs_per_cell;
+
+      Assert (&fe != 0, ExcNotInitialized());
+
+      CellData  cell_data (dofs_per_cell);
+      cell_data.quadrant = p4est_cell;
+      cell_data.tree_index = tree_index;
+      cell_data.receiver = dealii::internal::p4est::functions<dim>::
+                           comm_find_owner (const_cast<typename dealii::internal::p4est::types<dim>::forest *> (&forest),
+                                            tree_index,
+                                            &p4est_cell,
+                                            dealii_cell->level_subdomain_id ());
+
+      cell_data_insert (cell_data, new_needs);
+    }
+
+
+
+    template <int dim,int spacedim>
+    int
+    ExtrapolateImplementation<dim,spacedim>::
+    cell_data_search (const CellData               &cell_data,
+                      const std::vector<CellData>  &cells_list)
+    {
+      typename std::vector<CellData>::const_iterator  bound
+        = std::lower_bound (cells_list.begin(), cells_list.end(), cell_data);
+
+      if ((bound != cells_list.end ())
+          &&
+          !(cell_data < *bound))
+        return (int)(bound - cells_list.begin ());
+
+      return (-1);
+    }
+
+
+
+    template <int dim,int spacedim>
+    void
+    ExtrapolateImplementation<dim,spacedim>::
+    cell_data_insert (const CellData         &cell_data,
+                      std::vector<CellData>  &cells_list)
+    {
+      typename std::vector<CellData>::iterator  bound
+        = std::lower_bound (cells_list.begin(), cells_list.end(), cell_data);
+
+      if ((bound == cells_list.end ())
+          ||
+          (cell_data < *bound))
+        cells_list.insert (bound, 1, cell_data);
+    }
+
+
+
+    template <int dim,int spacedim>
+    template <class InVector>
+    void
+    ExtrapolateImplementation<dim,spacedim>::
+    compute_all_non_local_data (const DoFHandler<dim,spacedim> &dof2,
+                                const InVector                 &u)
+    {
+      std::vector<CellData>  cells_we_need,
+          cells_to_compute,
+          received_cells,
+          received_needs,
+          new_needs,
+          computed_cells,
+          cells_to_send;
+
+      // Compute all the cells needed
+      // from other processes.
+      compute_needs (dof2, cells_we_need);
+
+      // Send the cells needed to there
+      // owners and receive a list of cells other
+      // processes need from us.
+      send_cells (cells_we_need, received_needs);
+
+      // The list of received needs can contain
+      // some cells more than once because different
+      // processes may need data from the same cell.
+      // To compute data only once, generate a vector
+      // with unique entries and distribute the computed
+      // data afterwards back to a vector with correct
+      // receivers.
+      // Computing cell_data can cause a need for
+      // data from some new cells.
+      // If a cell is computed, send it back to
+      // their senders, maybe receive new needs and
+      // compute again, do not wait that all cells
+      // are computed or all needs are collected.
+      // Otherwise we can run into a deadlock if
+      // a cell needed from another process,
+      // itself needs some data from us.
+      unsigned int ready = 0;
+      do
+        {
+          for (unsigned int i=0; i<received_needs.size(); ++i)
+            cell_data_insert (received_needs[i], cells_to_compute);
+
+          compute_cells (dof2, u, cells_to_compute, computed_cells, new_needs);
+
+          // if there are no cells to compute and no new needs, stop
+          ready = Utilities::MPI::sum (new_needs.size () + cells_to_compute.size (), communicator);
+
+          for (typename std::vector<CellData>::const_iterator comp=computed_cells.begin ();
+               comp != computed_cells.end ();
+               ++comp)
+            {
+              // store computed cells
+              cell_data_insert (*comp, available_cells);
+
+              // and generate a vector
+              // of computed cells with
+              // correct receivers
+              // then delete this received
+              // need from the list
+              typename std::vector<CellData>::iterator recv=received_needs.begin();
+              while (recv != received_needs.end())
+                {
+                  if (dealii::internal::p4est::quadrant_is_equal<dim>(recv->quadrant, comp->quadrant))
+                    {
+                      recv->dof_values = comp->dof_values;
+                      cells_to_send.push_back (*recv);
+                      received_needs.erase (recv);
+                      recv = received_needs.begin();
+                    }
+                  else
+                    ++recv;
+                }
+            }
+
+          send_cells (cells_to_send, received_cells);
+
+          // strore received cell_data
+          for (typename std::vector<CellData>::const_iterator recv=received_cells.begin ();
+               recv != received_cells.end ();
+               ++recv)
+            {
+              cell_data_insert (*recv, available_cells);
+            }
+
+          // finally send and receive new
+          // needs and start a new round
+          send_cells (new_needs, received_needs);
+        }
+      while (ready != 0);
+    }
+
+
+
+    template <int dim,int spacedim>
+    template <class InVector,class OutVector>
+    void
+    ExtrapolateImplementation<dim,spacedim>::
+    extrapolate_parallel (const InVector &u2_relevant,
+                          const DoFHandler<dim,spacedim> &dof2,
+                          OutVector &u2)
+    {
+      const parallel::distributed::Triangulation< dim, spacedim > *tr
+        = (dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>
+           (&dof2.get_tria()));
+
+      Assert (tr != 0,
+              ExcMessage ("Extrapolate in parallel only works for parallel distributed triangulations!"));
+
+      communicator = tr->get_communicator ();
+
+      compute_all_non_local_data (dof2, u2_relevant);
+
+      // exclude dofs on more refined ghosted cells
+      const FiniteElement<dim,spacedim> &fe  = dof2.get_fe();
+      const unsigned int dofs_per_face = fe.dofs_per_face;
+      if (dofs_per_face > 0)
+        {
+          const unsigned int dofs_per_cell = fe.dofs_per_cell;
+          std::vector<types::global_dof_index> indices (dofs_per_cell);
+          typename DoFHandler<dim,spacedim>::active_cell_iterator
+          cell=dof2.begin_active(),
+          endc=dof2.end();
+          for (; cell!=endc; ++cell)
+            if (cell->is_ghost())
+              {
+                cell->get_dof_indices(indices);
+                for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+                  if (!cell->at_boundary(face))
+                    {
+                      const typename DoFHandler<dim,spacedim>::cell_iterator neighbor = cell->neighbor(face);
+                      if (neighbor->level() != cell->level())
+                        for (unsigned int i=0; i<dofs_per_face; ++i)
+                          {
+                            const types::global_dof_index index = indices[fe.face_to_cell_index(i,face)];;
+                            const bool index_stored
+                              = (dofs_on_refined_neighbors.find(index)!=dofs_on_refined_neighbors.end());
+                            const unsigned int level = index_stored ?
+                                                       std::max(cell->level(), dofs_on_refined_neighbors[index]) :
+                                                       cell->level();
+                            dofs_on_refined_neighbors[index] = level;
+                          }
+                    }
+              }
+        }
+
+      // after all dof values on
+      // not owned patch cells
+      // are computed, start
+      // the interpolation
+      u2 = 0;
+
+      std::queue<WorkPackage> queue;
+      {
+        typename DoFHandler<dim,spacedim>::cell_iterator
+        cell=dof2.begin(0),
+        endc=dof2.end(0);
+        for (; cell!=endc; ++cell)
+          {
+            if (dealii::internal::p4est::tree_exists_locally<dim> (tr->parallel_forest,
+                                                                   tr->coarse_cell_to_p4est_tree_permutation[cell->index()])
+                == false)
+              continue;
+
+            typename dealii::internal::p4est::types<dim>::quadrant p4est_coarse_cell;
+            const unsigned int tree_index = tr->coarse_cell_to_p4est_tree_permutation[cell->index()];
+            typename dealii::internal::p4est::types<dim>::tree *tree = tr->init_tree(cell->index());
+
+            dealii::internal::p4est::init_coarse_quadrant<dim>(p4est_coarse_cell);
+
+            const WorkPackage data (*tr->parallel_forest, *tree,tree_index,cell,p4est_coarse_cell);
+
+            queue.push(data);
+          }
+      }
+
+      while (!queue.empty())
+        {
+          const WorkPackage &data = queue.front();
+
+          const typename dealii::internal::p4est::types<dim>::forest    &forest = data.forest;
+          const typename dealii::internal::p4est::types<dim>::tree      &tree = data.tree;
+          const typename dealii::internal::p4est::types<dim>::locidx    &tree_index= data.tree_index;
+          const typename DoFHandler<dim,spacedim>::cell_iterator        &dealii_cell = data.dealii_cell;
+          const typename dealii::internal::p4est::types<dim>::quadrant  &p4est_cell = data.p4est_cell;
+
+          interpolate_recursively (forest, tree, tree_index, dealii_cell, p4est_cell, u2_relevant, u2);
+
+          // traverse recursively over this tree
+          if (dealii_cell->has_children())
+            {
+              Assert(dealii_cell->has_children(), ExcInternalError());
+              typename dealii::internal::p4est::types<dim>::quadrant
+              p4est_child[GeometryInfo<dim>::max_children_per_cell];
+
+              dealii::internal::p4est::init_quadrant_children<dim> (p4est_cell,
+                                                                    p4est_child);
+
+              for (unsigned int c=0; c<GeometryInfo<dim>::max_children_per_cell; ++c)
+                queue.push(WorkPackage(forest, tree, tree_index, dealii_cell->child(c), p4est_child[c]));
+            }
+          queue.pop();
+        }
+
+
+      u2.compress(VectorOperation::insert);
+    }
+#endif //DEAL_II_WITH_P4EST
+
+    namespace
+    {
+      template <int dim, class InVector, class OutVector, int spacedim>
+      void extrapolate_serial(const InVector &u3,
+                              const DoFHandler<dim,spacedim> &dof2,
+                              OutVector &u2)
+      {
+        const unsigned int dofs_per_cell  = dof2.get_fe().dofs_per_cell;
+        Vector<typename OutVector::value_type> dof_values(dofs_per_cell);
+
+        // then traverse grid bottom up
+        for (unsigned int level=0; level<dof2.get_triangulation().n_levels()-1; ++level)
+          {
+            typename DoFHandler<dim,spacedim>::cell_iterator cell=dof2.begin(level),
+                                                             endc=dof2.end(level);
+
+            for (; cell!=endc; ++cell)
+              if (!cell->active())
+                {
+                  // check whether this
+                  // cell has active
+                  // children
+                  bool active_children=false;
+                  for (unsigned int child_n=0; child_n<cell->n_children(); ++child_n)
+                    if (cell->child(child_n)->active())
+                      {
+                        active_children=true;
+                        break;
+                      }
+
+                  // if there are active
+                  // children, this process
+                  // has to work on this
+                  // cell. get the data
+                  // from the one vector
+                  // and set it on the
+                  // other
+                  if (active_children)
+                    {
+                      cell->get_interpolated_dof_values(u3, dof_values);
+                      cell->set_dof_values_by_interpolation(dof_values, u2);
+                    }
+                }
+          }
+      }
+
+
+
+#ifdef DEAL_II_WITH_P4EST
+      template <int dim, class InVector, class OutVector, int spacedim>
+      void extrapolate_parallel(const InVector &u2_relevant,
+                                const DoFHandler<dim,spacedim> &dof2,
+                                OutVector &u2)
+      {
+        if (dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>(&dof2.get_tria()) != 0)
+          {
+            ExtrapolateImplementation<dim,spacedim>  implementation;
+            implementation.extrapolate_parallel (u2_relevant, dof2, u2);
+          }
+        else
+          {
+            extrapolate_serial (u2_relevant, dof2, u2);
+          }
+      }
+#endif //DEAL_II_WITH_P4EST
+    }
+  }
+
+  template <int dim, class InVector, class OutVector, int spacedim>
+  void extrapolate(const DoFHandler<dim,spacedim> &dof1,
+                   const InVector &u1,
+                   const DoFHandler<dim,spacedim> &dof2,
+                   OutVector &u2)
+  {
+    ConstraintMatrix dummy;
+    dummy.close();
+    extrapolate(dof1, u1, dof2, dummy, u2);
+  }
+
+
+
+  template <int dim, class InVector, class OutVector, int spacedim>
+  void extrapolate(const DoFHandler<dim,spacedim> &dof1,
+                   const InVector &u1,
+                   const DoFHandler<dim,spacedim> &dof2,
+                   const ConstraintMatrix &constraints,
+                   OutVector &u2)
+  {
+    Assert(dof1.get_fe().n_components() == dof2.get_fe().n_components(),
+           ExcDimensionMismatch(dof1.get_fe().n_components(), dof2.get_fe().n_components()));
+    Assert(&dof1.get_triangulation()==&dof2.get_triangulation(), ExcTriangulationMismatch());
+    Assert(u1.size()==dof1.n_dofs(), ExcDimensionMismatch(u1.size(), dof1.n_dofs()));
+    Assert(u2.size()==dof2.n_dofs(), ExcDimensionMismatch(u2.size(), dof2.n_dofs()));
+
+    // make sure that each cell on the
+    // coarsest level is at least once
+    // refined, otherwise, these cells
+    // can't be treated and would
+    // generate a bogus result
+    {
+      typename DoFHandler<dim,spacedim>::cell_iterator cell = dof2.begin(0),
+                                                       endc = dof2.end(0);
+      for (; cell!=endc; ++cell)
+        Assert (cell->has_children(), ExcGridNotRefinedAtLeastOnce());
+    }
+
+    OutVector u3;
+    u3.reinit(u2);
+    interpolate(dof1, u1, dof2, constraints, u3);
+
+    internal::extrapolate_serial (u3, dof2, u2);
+
+    constraints.distribute(u2);
+  }
+
+
+#ifdef DEAL_II_WITH_P4EST
+  template <int dim, class VectorType, int spacedim>
+  void extrapolate_parallel(const DoFHandler<dim,spacedim> &dof1,
+                            const VectorType &u1,
+                            const DoFHandler<dim,spacedim> &dof2,
+                            const ConstraintMatrix &constraints2,
+                            VectorType &u2)
+  {
+    IndexSet  dof2_locally_owned_dofs = dof2.locally_owned_dofs();
+    IndexSet  dof2_locally_relevant_dofs;
+    DoFTools::extract_locally_relevant_dofs (dof2,
+                                             dof2_locally_relevant_dofs);
+
+    VectorType u3 (dof2_locally_owned_dofs, u1.get_mpi_communicator ());
+    interpolate (dof1, u1, dof2, constraints2, u3);
+    VectorType u3_relevant (dof2_locally_owned_dofs,
+                            dof2_locally_relevant_dofs,
+                            u1.get_mpi_communicator ());
+    u3_relevant = u3;
+
+    internal::extrapolate_parallel (u3_relevant, dof2, u2);
+
+    constraints2.distribute(u2);
+  }
+
+
+
+  template <int dim, int spacedim, typename Number>
+  void extrapolate_parallel(const DoFHandler<dim,spacedim> &dof1,
+                            const LinearAlgebra::distributed::Vector<Number> &u1,
+                            const DoFHandler<dim,spacedim> &dof2,
+                            const ConstraintMatrix &constraints2,
+                            LinearAlgebra::distributed::Vector<Number> &u2)
+  {
+    IndexSet  dof2_locally_owned_dofs = dof2.locally_owned_dofs();
+    IndexSet  dof2_locally_relevant_dofs;
+    DoFTools::extract_locally_relevant_dofs (dof2,
+                                             dof2_locally_relevant_dofs);
+
+    LinearAlgebra::distributed::Vector<Number> u3(dof2_locally_owned_dofs,
+                                                  dof2_locally_relevant_dofs,
+                                                  u1.get_mpi_communicator ());
+    interpolate (dof1, u1, dof2, constraints2, u3);
+    u3.update_ghost_values();
+
+    internal::extrapolate_parallel (u3, dof2, u2);
+
+    constraints2.distribute(u2);
+  }
+#endif //DEAL_II_WITH_P4EST
+
+} // end of namespace FETools
+
+
+/*-------------- Explicit Instantiations -------------------------------*/
+#include "fe_tools_extrapolate.inst"
+
+
+/*----------------------------   fe_tools.cc     ---------------------------*/
+
+DEAL_II_NAMESPACE_CLOSE
diff --git a/source/fe/fe_tools_extrapolate.inst.in b/source/fe/fe_tools_extrapolate.inst.in
new file mode 100755 (executable)
index 0000000..dbcb937
--- /dev/null
@@ -0,0 +1,90 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 1998 - 2014 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS; VEC : SERIAL_VECTORS)
+{
+    namespace FETools
+    \{
+#if deal_II_dimension == deal_II_space_dimension
+    template
+    void extrapolate<deal_II_dimension>
+    (const DoFHandler<deal_II_dimension> &, const VEC &,
+     const DoFHandler<deal_II_dimension> &, VEC &);
+
+    template
+    void extrapolate<deal_II_dimension>
+    (const DoFHandler<deal_II_dimension> &, const VEC &,
+     const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
+     VEC &);
+#endif
+    \}
+}
+
+
+
+for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS)
+{
+#ifdef DEAL_II_WITH_P4EST
+    namespace FETools
+    \{
+#if deal_II_dimension == deal_II_space_dimension
+
+#ifdef DEAL_II_WITH_PETSC
+    template
+    void extrapolate_parallel<deal_II_dimension, PETScWrappers::MPI::Vector, deal_II_space_dimension>
+    (const DoFHandler<deal_II_dimension,deal_II_space_dimension> &,
+     const PETScWrappers::MPI::Vector &,
+     const DoFHandler<deal_II_dimension,deal_II_space_dimension> &,
+     const ConstraintMatrix &,
+     PETScWrappers::MPI::Vector &);
+#endif
+
+#ifdef DEAL_II_WITH_TRILINOS
+    template
+    void extrapolate_parallel <deal_II_dimension, TrilinosWrappers::MPI::Vector, deal_II_space_dimension>
+    (const DoFHandler<deal_II_dimension,deal_II_space_dimension> &,
+     const TrilinosWrappers::MPI::Vector &,
+     const DoFHandler<deal_II_dimension,deal_II_space_dimension> &,
+     const ConstraintMatrix &,
+     TrilinosWrappers::MPI::Vector &);
+#endif
+
+#endif
+    \}
+#endif //DEAL_II_WITH_P4EST
+}
+
+
+
+for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS; scalar : REAL_SCALARS)
+{
+#ifdef DEAL_II_WITH_P4EST
+    namespace FETools
+    \{
+#if deal_II_dimension == deal_II_space_dimension
+    template
+    void extrapolate_parallel <deal_II_dimension, deal_II_space_dimension, scalar>
+    (const DoFHandler<deal_II_dimension,deal_II_space_dimension> &,
+     const LinearAlgebra::distributed::Vector<scalar> &,
+     const DoFHandler<deal_II_dimension,deal_II_space_dimension> &,
+     const ConstraintMatrix &,
+     LinearAlgebra::distributed::Vector<scalar> &);
+#endif
+    \}
+#endif //DEAL_II_WITH_P4EST
+}
+
old mode 100644 (file)
new mode 100755 (executable)
index 0b1d5cf..b82ece6
@@ -1,6 +1,6 @@
 // ---------------------------------------------------------------------
 //
-// Copyright (C) 2000 - 2015 by the deal.II authors
+// Copyright (C) 2000 - 2016 by the deal.II authors
 //
 // This file is part of the deal.II library.
 //
@@ -44,9 +44,7 @@
 
 #include <deal.II/base/index_set.h>
 
-#include <deal.II/distributed/p4est_wrappers.h>
 #include <iostream>
-#include <queue>
 
 
 DEAL_II_NAMESPACE_OPEN
@@ -718,1527 +716,6 @@ namespace FETools
         ++cell2;
       }
   }
-
-
-  namespace internal
-  {
-#ifdef DEAL_II_WITH_P4EST
-    // Implementation of the @p extrapolate function
-    // on parallel distributed grids.
-    template <int dim,int spacedim>
-    class ExtrapolateImplementation
-    {
-    public:
-
-      ExtrapolateImplementation ();
-
-      ~ExtrapolateImplementation ();
-
-      template <class InVector,class OutVector>
-      void extrapolate_parallel (const InVector &u2_relevant,
-                                 const DoFHandler<dim,spacedim> &dof2,
-                                 OutVector &u2);
-
-    private:
-
-      // A structure holding all data to
-      // set dofs recursively on cells of arbitrary level
-      struct WorkPackage
-      {
-        const typename dealii::internal::p4est::types<dim>::forest    forest;
-        const typename dealii::internal::p4est::types<dim>::tree      tree;
-        const typename dealii::internal::p4est::types<dim>::locidx    tree_index;
-        const typename DoFHandler<dim,spacedim>::cell_iterator        dealii_cell;
-        const typename dealii::internal::p4est::types<dim>::quadrant  p4est_cell;
-
-        WorkPackage(const typename dealii::internal::p4est::types<dim>::forest    &forest_,
-                    const typename dealii::internal::p4est::types<dim>::tree      &tree_,
-                    const typename dealii::internal::p4est::types<dim>::locidx    &tree_index_,
-                    const typename DoFHandler<dim,spacedim>::cell_iterator        &dealii_cell_,
-                    const typename dealii::internal::p4est::types<dim>::quadrant  &p4est_cell_)
-          :
-          forest(forest_),
-          tree(tree_),
-          tree_index(tree_index_),
-          dealii_cell(dealii_cell_),
-          p4est_cell(p4est_cell_)
-        {}
-      };
-
-
-      // A structure holding all data
-      // of cells needed from other processes
-      // for the extrapolate algorithm.
-      struct CellData
-      {
-        CellData ();
-
-        CellData (const unsigned int dofs_per_cell);
-
-        Vector<double>  dof_values;
-
-        unsigned int tree_index;
-
-        typename dealii::internal::p4est::types<dim>::quadrant  quadrant;
-
-        int receiver;
-
-        bool operator < (const CellData &rhs) const
-        {
-          if (dealii::internal::p4est::functions<dim>::quadrant_compare (&quadrant, &rhs.quadrant) < 0)
-            return true;
-
-          return false;
-        }
-
-        unsigned int bytes_for_buffer () const
-        {
-          return (sizeof(unsigned int) +                                              // dofs_per_cell
-                  dof_values.size() * sizeof(double) +                                // dof_values
-                  sizeof(unsigned int) +                                              // tree_index
-                  sizeof(typename dealii::internal::p4est::types<dim>::quadrant));    // quadrant
-        }
-
-        void pack_data (std::vector<char> &buffer) const
-        {
-          buffer.resize(bytes_for_buffer());
-
-          char *ptr = &buffer[0];
-
-          unsigned int n_dofs = dof_values.size ();
-          std::memcpy(ptr, &n_dofs, sizeof(unsigned int));
-          ptr += sizeof(unsigned int);
-
-          std::memcpy(ptr,dof_values.begin(),n_dofs*sizeof(double));
-          ptr += n_dofs*sizeof(double);
-
-          std::memcpy(ptr,&tree_index,sizeof(unsigned int));
-          ptr += sizeof(unsigned int);
-
-          std::memcpy(ptr,&quadrant,sizeof(typename dealii::internal::p4est::types<dim>::quadrant));
-          ptr += sizeof(typename dealii::internal::p4est::types<dim>::quadrant);
-
-          Assert (ptr == &buffer[0]+buffer.size(),
-                  ExcInternalError());
-        }
-
-        void unpack_data (const std::vector<char> &buffer)
-        {
-          const char *ptr = &buffer[0];
-          unsigned int n_dofs;
-          memcpy(&n_dofs, ptr, sizeof(unsigned int));
-          ptr += sizeof(unsigned int);
-
-          dof_values.reinit(n_dofs);
-          std::memcpy(dof_values.begin(),ptr,n_dofs * sizeof(double));
-          ptr += n_dofs * sizeof(double);
-
-          std::memcpy(&tree_index,ptr,sizeof(unsigned int));
-          ptr += sizeof(unsigned int);
-
-          std::memcpy(&quadrant,ptr,sizeof(typename dealii::internal::p4est::types<dim>::quadrant));
-          ptr += sizeof(typename dealii::internal::p4est::types<dim>::quadrant);
-
-          Assert (ptr == &buffer[0]+buffer.size(),
-                  ExcInternalError());
-        }
-      };
-
-      // Problem: The function extrapolates a polynomial
-      // function from a finer mesh of size $h$ to a polynmial
-      // function of higher degree but on a coarser mesh of
-      // size $2h$. Therefor the mesh has to consist of patches
-      // of four (2d) or eight (3d) cells and the function requires
-      // that the mesh is refined globally at least once.
-      // The algorithm starts on the coarsest level of the grid,
-      // loops over all cells and if a cell has at least one active child,
-      // dof values are set via
-      //   cell->get_interpolated_dof_values(u_input, dof_values)
-      //   cell->set_dof_values_by_interpolation(dof_values, u_output)
-      // both *_interpolation_* functions traverse recursively
-      // over all children down to the active cells
-      // and get/set dof values by interpolation.
-      //
-      // On distributed parallel grids the problem arises, that
-      // if a cell has at least one active child, there is no guarantee
-      // that all children of this cell belong to this process.
-      // There might be children which are owned by other processes
-      // and the algorithm needs to find and has to get the
-      // dof values from these processes and so on.
-      //
-      // Algorithm:
-      // 1) Loop over all coarse cells
-      //      From each coarse cell traverse down the tree
-      //      of refined cells and search for active children
-      //      If there is an active child, check, if all
-      //      other children down to the finest level are part
-      //      of this process, if not, add the cell to the list
-      //      of needs.
-      // 2) Send the list of needs to other processes
-      //      Each process has a list of needs and after
-      //      the loop over all coarse cells (all trees)
-      //      is finished, this list can be send to other processes.
-      // 3) Compute the needs required by other processes
-      //      While computing the values required from other
-      //      processes there can arise new needs and they are
-      //      stored in a list again.
-      // 4) Send the computed values and the list of new needs around
-      //
-      // This procedure has to be repeated until no process needs any
-      // new need and all needs are computed, but there are at most the
-      // "number of grid levels" rounds of sending/receiving cell data.
-      //
-      // Then each process has all data needed for extrapolation.
-
-      // driver function sending/receiving all values from
-      // different processes
-      template <class InVector>
-      void compute_all_non_local_data (const DoFHandler<dim,spacedim> &dof2,
-                                       const InVector                 &u);
-
-      // traverse recursively over
-      // the cells of this tree and
-      // interpolate over patches which
-      // are part of this process
-      template <class InVector,class OutVector>
-      void interpolate_recursively (const typename dealii::internal::p4est::types<dim>::forest    &forest,
-                                    const typename dealii::internal::p4est::types<dim>::tree      &tree,
-                                    const typename dealii::internal::p4est::types<dim>::locidx    &tree_index,
-                                    const typename DoFHandler<dim,spacedim>::cell_iterator        &dealii_cell,
-                                    const typename dealii::internal::p4est::types<dim>::quadrant  &p4est_cell,
-                                    const InVector                                                &u1,
-                                    OutVector                                                     &u2);
-
-      // get dof values for this
-      // cell by interpolation
-      // if a child is reached which
-      // is not part of this process
-      // a new need is created
-      template <class InVector,typename number>
-      void get_interpolated_dof_values (const typename dealii::internal::p4est::types<dim>::forest    &forest,
-                                        const typename dealii::internal::p4est::types<dim>::tree      &tree,
-                                        const typename dealii::internal::p4est::types<dim>::locidx    &tree_index,
-                                        const typename DoFHandler<dim,spacedim>::cell_iterator        &dealii_cell,
-                                        const typename dealii::internal::p4est::types<dim>::quadrant  &p4est_cell,
-                                        const InVector                                                &u,
-                                        Vector<number>                                                &interpolated_values,
-                                        std::vector<CellData>                                         &new_needs);
-
-      // set dof values for this
-      // cell by interpolation
-      template <class OutVector,typename number>
-      void set_dof_values_by_interpolation (const typename DoFHandler<dim,spacedim>::cell_iterator        &dealii_cell,
-                                            const typename dealii::internal::p4est::types<dim>::quadrant  &p4est_cell,
-                                            const Vector<number>                                          &interpolated_values,
-                                            OutVector                                                     &u);
-
-      // compute all cell_data
-      // needed from other processes
-      // to interpolate the part of
-      // this process
-      void compute_needs (const DoFHandler<dim,spacedim> &dof2,
-                          std::vector<CellData>          &new_needs);
-
-      // traverse over the tree
-      // and look for patches this
-      // process has to work on
-      void traverse_tree_recursively (const typename dealii::internal::p4est::types<dim>::forest    &forest,
-                                      const typename dealii::internal::p4est::types<dim>::tree      &tree,
-                                      const typename dealii::internal::p4est::types<dim>::locidx    &tree_index,
-                                      const typename DoFHandler<dim,spacedim>::cell_iterator        &dealii_cell,
-                                      const typename dealii::internal::p4est::types<dim>::quadrant  &p4est_cell,
-                                      std::vector<CellData>                                         &new_needs);
-
-      // traverse recursively
-      // over a patch and look
-      // for cells needed from
-      // other processes for interpolation
-      void traverse_patch_recursively (const typename dealii::internal::p4est::types<dim>::forest    &forest,
-                                       const typename dealii::internal::p4est::types<dim>::tree      &tree,
-                                       const typename dealii::internal::p4est::types<dim>::locidx    &tree_index,
-                                       const typename DoFHandler<dim,spacedim>::cell_iterator        &dealii_cell,
-                                       const typename dealii::internal::p4est::types<dim>::quadrant  &p4est_cell,
-                                       std::vector<CellData>                                         &new_needs);
-
-      // compute dof values of all
-      // cells collected in cells_to_compute
-      // computed cells are deleted
-      // from cells_to_compute and
-      // stored in computed_cells
-      // during computation there can
-      // be new cells needed from
-      // other processes, these cells
-      // are stored in new_needs
-      template <class InVector>
-      void compute_cells (const DoFHandler<dim,spacedim> &dof2,
-                          const InVector                 &u,
-                          std::vector<CellData>          &cells_to_compute,
-                          std::vector<CellData>          &computed_cells,
-                          std::vector<CellData>          &new_needs);
-
-      // traverse over the tree
-      // and compute cells
-      template <class InVector>
-      void compute_cells_in_tree_recursively (const typename dealii::internal::p4est::types<dim>::forest    &forest,
-                                              const typename dealii::internal::p4est::types<dim>::tree      &tree,
-                                              const typename dealii::internal::p4est::types<dim>::locidx    &tree_index,
-                                              const typename DoFHandler<dim,spacedim>::cell_iterator        &dealii_cell,
-                                              const typename dealii::internal::p4est::types<dim>::quadrant  &p4est_cell,
-                                              const InVector                                                &u,
-                                              std::vector<CellData>                                         &cells_to_compute,
-                                              std::vector<CellData>                                         &computed_cells,
-                                              std::vector<CellData>                                         &new_needs);
-
-      // send all cell_data in the vector
-      // cells_to_send to their receivers
-      // and receives a vector of cell_data
-      void send_cells (const std::vector<CellData>  &cells_to_send,
-                       std::vector<CellData>  &received_cells) const;
-
-      // add new cell_data to
-      // the ordered list new_needs
-      // uses cell_data_insert
-      void add_new_need (const typename dealii::internal::p4est::types<dim>::forest    &forest,
-                         const typename dealii::internal::p4est::types<dim>::locidx    &tree_index,
-                         const typename DoFHandler<dim,spacedim>::cell_iterator        &dealii_cell,
-                         const typename dealii::internal::p4est::types<dim>::quadrant  &p4est_cell,
-                         std::vector<CellData>                                         &new_needs) const;
-
-      // binary search in cells_list
-      // assume that cells_list
-      // is ordered
-      static int cell_data_search (const CellData               &cell_data,
-                                   const std::vector<CellData>  &cells_list);
-
-      // insert cell_data into a sorted
-      // vector cells_list at the correct
-      // position if cell_data
-      // not exists already in cells_list
-      static void cell_data_insert (const CellData         &cell_data,
-                                    std::vector<CellData>  &cells_list);
-
-      MPI_Comm  communicator;
-
-      // a vector of all cells this process has
-      // computed or received data
-      std::vector<CellData>  available_cells;
-
-      // stores the indices of dofs on more refined ghosted cells along
-      // with the maximum level
-      std::map<types::global_dof_index, int> dofs_on_refined_neighbors;
-    };
-
-    template <>
-    class ExtrapolateImplementation<1,1>
-    {
-    public:
-
-      ExtrapolateImplementation ()
-      {
-        AssertThrow(false, ExcNotImplemented())
-      }
-
-      ~ExtrapolateImplementation ()
-      {}
-
-      template <class InVector,class OutVector>
-      void extrapolate_parallel (const InVector &/*u2_relevant*/,
-                                 const DoFHandler<1,1> &/*dof2*/,
-                                 OutVector &/*u2*/)
-      {}
-    };
-
-    template <>
-    class ExtrapolateImplementation<1,2>
-    {
-    public:
-
-      ExtrapolateImplementation ()
-      {
-        AssertThrow(false, ExcNotImplemented())
-      }
-
-      ~ExtrapolateImplementation ()
-      {}
-
-      template <class InVector,class OutVector>
-      void extrapolate_parallel (const InVector &/*u2_relevant*/,
-                                 const DoFHandler<1,2> &/*dof2*/,
-                                 OutVector &/*u2*/)
-      {}
-    };
-
-    template <>
-    class ExtrapolateImplementation<1,3>
-    {
-    public:
-
-      ExtrapolateImplementation ()
-      {
-        AssertThrow(false, ExcNotImplemented())
-      }
-
-      ~ExtrapolateImplementation ()
-      {}
-
-      template <class InVector,class OutVector>
-      void extrapolate_parallel (const InVector &/*u2_relevant*/,
-                                 const DoFHandler<1,3> &/*dof2*/,
-                                 OutVector &/*u2*/)
-      {}
-    };
-
-
-
-    template <int dim,int spacedim>
-    ExtrapolateImplementation<dim,spacedim>::
-    ExtrapolateImplementation ()
-    {}
-
-
-
-    template <int dim,int spacedim>
-    ExtrapolateImplementation<dim,spacedim>::
-    ~ExtrapolateImplementation ()
-    {}
-
-
-
-    template <int dim,int spacedim>
-    ExtrapolateImplementation<dim,spacedim>::
-    CellData::CellData ()
-      : tree_index (0),
-        receiver (0)
-    {}
-
-
-
-    template <int dim,int spacedim>
-    ExtrapolateImplementation<dim,spacedim>::
-    CellData::CellData (const unsigned int dofs_per_cell)
-      : tree_index (0),
-        receiver (0)
-    {
-      dof_values.reinit (dofs_per_cell);
-    }
-
-
-
-    template <int dim,int spacedim>
-    template <class InVector,class OutVector>
-    void
-    ExtrapolateImplementation<dim,spacedim>::
-    interpolate_recursively (const typename dealii::internal::p4est::types<dim>::forest    &forest,
-                             const typename dealii::internal::p4est::types<dim>::tree      &tree,
-                             const typename dealii::internal::p4est::types<dim>::locidx    &tree_index,
-                             const typename DoFHandler<dim,spacedim>::cell_iterator        &dealii_cell,
-                             const typename dealii::internal::p4est::types<dim>::quadrant  &p4est_cell,
-                             const InVector                                                &u1,
-                             OutVector                                                     &u2)
-    {
-      // check if this cell exists in the local p4est
-      int idx = sc_array_bsearch(const_cast<sc_array_t *>(&tree.quadrants),
-                                 &p4est_cell,
-                                 dealii::internal::p4est::functions<dim>::quadrant_compare);
-
-      // if neither this cell nor one of it's children belongs to us, don't do anything
-      if (idx == -1 && (dealii::internal::p4est::functions<dim>::
-                        quadrant_overlaps_tree (const_cast<typename dealii::internal::p4est::types<dim>::tree *>(&tree),
-                                                &p4est_cell)
-                        == false))
-        return;
-
-      bool p4est_has_children = (idx == -1);
-
-      bool locally_owned_children = false;
-      if (p4est_has_children)
-        {
-          Assert (dealii_cell->has_children (), ExcInternalError ());
-
-          // check if at least one child is locally owned on our process
-          for (unsigned int child_n=0; child_n<dealii_cell->n_children(); ++child_n)
-            if (dealii_cell->child(child_n)->active())
-              if (dealii_cell->child(child_n)->is_locally_owned())
-                {
-                  locally_owned_children=true;
-                  break;
-                }
-        }
-
-      if (locally_owned_children)
-        {
-          const FiniteElement<dim,spacedim> &fe            = dealii_cell->get_dof_handler().get_fe();
-          const unsigned int                 dofs_per_cell = fe.dofs_per_cell;
-
-          Vector<typename OutVector::value_type> interpolated_values(dofs_per_cell);
-
-          std::vector<CellData>  new_needs;
-          get_interpolated_dof_values (forest,
-                                       tree,
-                                       tree_index,
-                                       dealii_cell,
-                                       p4est_cell,
-                                       u1,
-                                       interpolated_values,
-                                       new_needs);
-
-          // at this point of
-          // the procedure no new
-          // needs should come up
-          Assert (new_needs.size () == 0,
-                  ExcInternalError ());
-
-          set_dof_values_by_interpolation (dealii_cell,
-                                           p4est_cell,
-                                           interpolated_values,
-                                           u2);
-        }
-
-
-
-    }
-
-
-
-    template <int dim,int spacedim>
-    template <class InVector,typename number>
-    void
-    ExtrapolateImplementation<dim,spacedim>::
-    get_interpolated_dof_values (const typename dealii::internal::p4est::types<dim>::forest    &forest,
-                                 const typename dealii::internal::p4est::types<dim>::tree      &tree,
-                                 const typename dealii::internal::p4est::types<dim>::locidx    &tree_index,
-                                 const typename DoFHandler<dim,spacedim>::cell_iterator        &dealii_cell,
-                                 const typename dealii::internal::p4est::types<dim>::quadrant  &p4est_cell,
-                                 const InVector                                                &u,
-                                 Vector<number>                                                &interpolated_values,
-                                 std::vector<CellData>                                         &new_needs)
-    {
-      if (!dealii_cell->has_children ())
-        {
-          if (dealii_cell->is_locally_owned ())
-            {
-              // if this is one of our cells,
-              // get dof values from input vector
-              dealii_cell->get_dof_values (u, interpolated_values);
-            }
-          else
-            {
-              add_new_need (forest,
-                            tree_index,
-                            dealii_cell,
-                            p4est_cell,
-                            new_needs);
-            }
-        }
-      else
-        {
-          const FiniteElement<dim,spacedim> &fe            = dealii_cell->get_dof_handler().get_fe();
-          const unsigned int                 dofs_per_cell = fe.dofs_per_cell;
-
-          Assert (&fe != 0,
-                  ExcNotInitialized());
-          Assert (interpolated_values.size() == dofs_per_cell,
-                  ExcDimensionMismatch(interpolated_values.size(), dofs_per_cell));
-          Assert (u.size() == dealii_cell->get_dof_handler().n_dofs(),
-                  ExcDimensionMismatch(u.size(), dealii_cell->get_dof_handler().n_dofs()));
-
-          Vector<number> tmp1(dofs_per_cell);
-          Vector<number> tmp2(dofs_per_cell);
-
-          interpolated_values = 0;
-          std::vector<bool> restriction_is_additive (dofs_per_cell);
-          for (unsigned int i=0; i<dofs_per_cell; ++i)
-            restriction_is_additive[i] = fe.restriction_is_additive(i);
-
-          typename dealii::internal::p4est::types<dim>::quadrant
-          p4est_child[GeometryInfo<dim>::max_children_per_cell];
-
-          dealii::internal::p4est::init_quadrant_children<dim> (p4est_cell,
-                                                                p4est_child);
-
-          bool found_child = true;
-          for (unsigned int c=0; c<GeometryInfo<dim>::max_children_per_cell; ++c)
-            {
-              if (dealii::internal::p4est::functions<dim>::
-                  quadrant_overlaps_tree (const_cast<typename dealii::internal::p4est::types<dim>::tree *>(&tree),
-                                          &p4est_child[c])
-                  == false)
-                {
-                  // this is a cell this process needs
-                  // data from another process
-
-                  // check if this cell is
-                  // available from other
-                  // computed patches
-                  CellData  cell_data;
-                  cell_data.quadrant = p4est_child[c];
-                  int pos = cell_data_search (cell_data, available_cells);
-
-                  if (pos == -1)
-                    {
-                      // data is not available
-                      // create a new need
-                      found_child = false;
-
-                      add_new_need (forest,
-                                    tree_index,
-                                    dealii_cell->child(c),
-                                    p4est_child[c],
-                                    new_needs);
-                    }
-                  else
-                    {
-                      Assert (available_cells[pos].dof_values.size() == dofs_per_cell,
-                              ExcDimensionMismatch(available_cells[pos].dof_values.size(), dofs_per_cell));
-
-                      tmp1 = available_cells[pos].dof_values;
-                    }
-                }
-              else
-                {
-                  // get the values from the present child, if necessary by
-                  // interpolation either from its own children or
-                  // by interpolating from the finite element on an active
-                  // child to the finite element space requested here
-                  get_interpolated_dof_values (forest,
-                                               tree,
-                                               tree_index,
-                                               dealii_cell->child(c),
-                                               p4est_child[c],
-                                               u,
-                                               tmp1,
-                                               new_needs);
-                }
-
-              if (found_child)
-                {
-                  // interpolate these to the mother cell
-                  fe.get_restriction_matrix(c, dealii_cell->refinement_case()).vmult (tmp2, tmp1);
-
-                  // and add up or set them in the output vector
-                  for (unsigned int i=0; i<dofs_per_cell; ++i)
-                    if (tmp2(i) != number())
-                      {
-                        if (restriction_is_additive[i])
-                          interpolated_values(i) += tmp2(i);
-                        else
-                          interpolated_values(i) = tmp2(i);
-                      }
-                }
-            }
-
-          if (found_child == false)
-            interpolated_values = 0;
-        }
-    }
-
-
-
-    template <int dim,int spacedim>
-    template <class OutVector,typename number>
-    void
-    ExtrapolateImplementation<dim,spacedim>::
-    set_dof_values_by_interpolation (const typename DoFHandler<dim,spacedim>::cell_iterator        &dealii_cell,
-                                     const typename dealii::internal::p4est::types<dim>::quadrant  &p4est_cell,
-                                     const Vector<number>                                          &local_values,
-                                     OutVector                                                     &u)
-    {
-      const FiniteElement<dim,spacedim> &fe            = dealii_cell->get_dof_handler().get_fe();
-      const unsigned int                 dofs_per_cell = fe.dofs_per_cell;
-
-      if (!dealii_cell->has_children ())
-        {
-          if (dealii_cell->is_locally_owned ())
-            {
-              std::vector<types::global_dof_index> indices(dofs_per_cell);
-              dealii_cell->get_dof_indices(indices);
-              for (unsigned int j=0; j<dofs_per_cell; ++j)
-                {
-                  // don't set this dof if there is a more refined ghosted neighbor setting this dof.
-                  const bool on_refined_neighbor
-                    = (dofs_on_refined_neighbors.find(indices[j])!=dofs_on_refined_neighbors.end());
-                  if (!(on_refined_neighbor && dofs_on_refined_neighbors[indices[j]]>dealii_cell->level()))
-                    u(indices[j]) = local_values(j);
-                }
-            }
-        }
-      else
-        {
-          Assert (&fe != 0,
-                  ExcNotInitialized());
-          Assert (local_values.size() == dofs_per_cell,
-                  ExcDimensionMismatch(local_values.size(), dofs_per_cell));
-          Assert (u.size() == dealii_cell->get_dof_handler().n_dofs(),
-                  ExcDimensionMismatch(u.size(), dealii_cell->get_dof_handler().n_dofs()));
-
-          Vector<number> tmp(dofs_per_cell);
-
-          typename dealii::internal::p4est::types<dim>::quadrant
-          p4est_child[GeometryInfo<dim>::max_children_per_cell];
-
-          dealii::internal::p4est::init_quadrant_children<dim> (p4est_cell,
-                                                                p4est_child);
-
-          for (unsigned int c=0; c<GeometryInfo<dim>::max_children_per_cell; ++c)
-            {
-              if (tmp.size() > 0)
-                fe.get_prolongation_matrix(c, dealii_cell->refinement_case())
-                .vmult (tmp, local_values);
-
-              set_dof_values_by_interpolation (dealii_cell->child(c),
-                                               p4est_child[c],
-                                               tmp,
-                                               u);
-            }
-        }
-    }
-
-
-
-    template <int dim,int spacedim>
-    void
-    ExtrapolateImplementation<dim,spacedim>::
-    compute_needs (const DoFHandler<dim,spacedim> &dof2,
-                   std::vector<CellData>          &new_needs)
-    {
-      const parallel::distributed::Triangulation< dim, spacedim > *tr
-        = (dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>
-           (&dof2.get_tria()));
-
-      Assert (tr != 0, ExcInternalError());
-
-      typename DoFHandler<dim,spacedim>::cell_iterator
-      cell=dof2.begin(0),
-      endc=dof2.end(0);
-      for (; cell!=endc; ++cell)
-        {
-          if (dealii::internal::p4est::tree_exists_locally<dim> (tr->parallel_forest,
-                                                                 tr->coarse_cell_to_p4est_tree_permutation[cell->index()])
-              == false)
-            continue;
-
-          typename dealii::internal::p4est::types<dim>::quadrant p4est_coarse_cell;
-          const unsigned int tree_index = tr->coarse_cell_to_p4est_tree_permutation[cell->index()];
-          typename dealii::internal::p4est::types<dim>::tree *tree = tr->init_tree(cell->index());
-
-          dealii::internal::p4est::init_coarse_quadrant<dim>(p4est_coarse_cell);
-
-          // make sure that each cell on the
-          // coarsest level is at least once
-          // refined, otherwise, these cells
-          // can't be treated and would
-          // generate a bogus result
-          {
-            int idx = sc_array_bsearch(const_cast<sc_array_t *>(&tree->quadrants),
-                                       &p4est_coarse_cell,
-                                       dealii::internal::p4est::functions<dim>::quadrant_compare);
-
-            AssertThrow (idx == -1, ExcGridNotRefinedAtLeastOnce ());
-          }
-
-          traverse_tree_recursively (*tr->parallel_forest,
-                                     *tree,
-                                     tree_index,
-                                     cell,
-                                     p4est_coarse_cell,
-                                     new_needs);
-        }
-    }
-
-
-
-    template <int dim,int spacedim>
-    void
-    ExtrapolateImplementation<dim,spacedim>::
-    traverse_tree_recursively (const typename dealii::internal::p4est::types<dim>::forest    &forest,
-                               const typename dealii::internal::p4est::types<dim>::tree      &tree,
-                               const typename dealii::internal::p4est::types<dim>::locidx    &tree_index,
-                               const typename DoFHandler<dim,spacedim>::cell_iterator        &dealii_cell,
-                               const typename dealii::internal::p4est::types<dim>::quadrant  &p4est_cell,
-                               std::vector<CellData>                                         &new_needs)
-    {
-      // check if this cell exists in the local p4est
-      int idx = sc_array_bsearch(const_cast<sc_array_t *>(&tree.quadrants),
-                                 &p4est_cell,
-                                 dealii::internal::p4est::functions<dim>::quadrant_compare);
-
-      // if neither this cell nor one of it's children belongs to us, don't do anything
-      if (idx == -1 && (dealii::internal::p4est::functions<dim>::
-                        quadrant_overlaps_tree (const_cast<typename dealii::internal::p4est::types<dim>::tree *>(&tree),
-                                                &p4est_cell)
-                        == false))
-        return;
-
-      bool p4est_has_children = (idx == -1);
-
-      // this cell is part of a patch
-      // this process has to interpolate on
-      // if there is at least one locally
-      // owned child
-      bool locally_owned_children = false;
-      if (p4est_has_children)
-        {
-          Assert (dealii_cell->has_children (), ExcInternalError ());
-
-          for (unsigned int child_n=0; child_n<dealii_cell->n_children(); ++child_n)
-            if (dealii_cell->child(child_n)->active())
-              if (dealii_cell->child(child_n)->is_locally_owned())
-                {
-                  locally_owned_children=true;
-                  break;
-                }
-        }
-
-      // traverse the patch recursively
-      // to find cells needed from
-      // other processes
-      if (locally_owned_children)
-        {
-          traverse_patch_recursively (forest,
-                                      tree,
-                                      tree_index,
-                                      dealii_cell,
-                                      p4est_cell,
-                                      new_needs);
-        }
-
-      // traverse recursively over this tree
-      if (p4est_has_children)
-        {
-          typename dealii::internal::p4est::types<dim>::quadrant
-          p4est_child[GeometryInfo<dim>::max_children_per_cell];
-
-          dealii::internal::p4est::init_quadrant_children<dim> (p4est_cell,
-                                                                p4est_child);
-
-          for (unsigned int c=0; c<GeometryInfo<dim>::max_children_per_cell; ++c)
-            {
-              traverse_tree_recursively (forest,
-                                         tree,
-                                         tree_index,
-                                         dealii_cell->child(c),
-                                         p4est_child[c],
-                                         new_needs);
-            }
-        }
-    }
-
-
-
-    template <int dim,int spacedim>
-    void
-    ExtrapolateImplementation<dim,spacedim>::
-    traverse_patch_recursively (const typename dealii::internal::p4est::types<dim>::forest    &forest,
-                                const typename dealii::internal::p4est::types<dim>::tree      &tree,
-                                const typename dealii::internal::p4est::types<dim>::locidx    &tree_index,
-                                const typename DoFHandler<dim,spacedim>::cell_iterator        &dealii_cell,
-                                const typename dealii::internal::p4est::types<dim>::quadrant  &p4est_cell,
-                                std::vector<CellData>                                         &new_needs)
-    {
-      if (dealii_cell->has_children ())
-        {
-          Assert (dealii_cell->has_children (), ExcInternalError ());
-
-          typename dealii::internal::p4est::types<dim>::quadrant
-          p4est_child[GeometryInfo<dim>::max_children_per_cell];
-
-          dealii::internal::p4est::init_quadrant_children<dim> (p4est_cell,
-                                                                p4est_child);
-
-          for (unsigned int c=0; c<GeometryInfo<dim>::max_children_per_cell; ++c)
-            {
-              // check if this child
-              // is locally available
-              // in the p4est
-              if (dealii::internal::p4est::functions<dim>::
-                  quadrant_overlaps_tree (const_cast<typename dealii::internal::p4est::types<dim>::tree *>(&tree),
-                                          &p4est_child[c])
-                  == false)
-                {
-                  // this is a cell for which this process
-                  // needs data from another process
-                  // so add the cell to the list
-                  add_new_need (forest,
-                                tree_index,
-                                dealii_cell->child(c),
-                                p4est_child[c],
-                                new_needs);
-                }
-              else
-                {
-                  // at least some part of
-                  // the tree rooted in this
-                  // child is locally available
-                  traverse_patch_recursively (forest,
-                                              tree,
-                                              tree_index,
-                                              dealii_cell->child(c),
-                                              p4est_child[c],
-                                              new_needs);
-                }
-            }
-        }
-    }
-
-
-
-
-    template <int dim,int spacedim>
-    template <class InVector>
-    void
-    ExtrapolateImplementation<dim,spacedim>::
-    compute_cells (const DoFHandler<dim,spacedim> &dof2,
-                   const InVector                 &u,
-                   std::vector<CellData>          &cells_to_compute,
-                   std::vector<CellData>          &computed_cells,
-                   std::vector<CellData>          &new_needs)
-    {
-      const parallel::distributed::Triangulation< dim, spacedim > *tr
-        = (dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>
-           (&dof2.get_tria()));
-
-      Assert (tr != 0, ExcInternalError());
-
-      // collect in a set all trees this
-      // process has to compute cells on
-      std::set<unsigned int>  trees;
-      for (typename std::vector<CellData>::const_iterator it=cells_to_compute.begin(); it!=cells_to_compute.end(); ++it)
-        trees.insert (it->tree_index);
-
-      typename DoFHandler<dim,spacedim>::cell_iterator
-      cell=dof2.begin(0),
-      endc=dof2.end(0);
-      for (; cell!=endc; ++cell)
-        {
-          // check if this is a tree this process has to
-          // work on and that this tree is in the p4est
-          const unsigned int tree_index = tr->coarse_cell_to_p4est_tree_permutation[cell->index()];
-
-          if ((trees.find (tree_index) == trees.end ())
-              ||
-              (dealii::internal::p4est::tree_exists_locally<dim> (tr->parallel_forest,
-                                                                  tree_index)
-               == false))
-            continue;
-
-          typename dealii::internal::p4est::types<dim>::quadrant p4est_coarse_cell;
-          typename dealii::internal::p4est::types<dim>::tree *tree = tr->init_tree(cell->index());
-
-          dealii::internal::p4est::init_coarse_quadrant<dim>(p4est_coarse_cell);
-
-          compute_cells_in_tree_recursively (*tr->parallel_forest,
-                                             *tree,
-                                             tree_index,
-                                             cell,
-                                             p4est_coarse_cell,
-                                             u,
-                                             cells_to_compute,
-                                             computed_cells,
-                                             new_needs);
-        }
-    }
-
-
-
-    template <int dim,int spacedim>
-    template <class InVector>
-    void
-    ExtrapolateImplementation<dim,spacedim>::
-    compute_cells_in_tree_recursively (const typename dealii::internal::p4est::types<dim>::forest    &forest,
-                                       const typename dealii::internal::p4est::types<dim>::tree      &tree,
-                                       const typename dealii::internal::p4est::types<dim>::locidx    &tree_index,
-                                       const typename DoFHandler<dim,spacedim>::cell_iterator        &dealii_cell,
-                                       const typename dealii::internal::p4est::types<dim>::quadrant  &p4est_cell,
-                                       const InVector                                                &u,
-                                       std::vector<CellData>                                         &cells_to_compute,
-                                       std::vector<CellData>                                         &computed_cells,
-                                       std::vector<CellData>                                         &new_needs)
-    {
-      if (cells_to_compute.size () == 0)
-        return;
-
-      // check if this cell exists in the local p4est
-      int idx = sc_array_bsearch(const_cast<sc_array_t *>(&tree.quadrants),
-                                 &p4est_cell,
-                                 dealii::internal::p4est::functions<dim>::quadrant_compare);
-
-      // if neither this cell nor one one of it's children belongs to us, don't do anything
-      if (idx == -1 && (dealii::internal::p4est::functions<dim>::
-                        quadrant_overlaps_tree (const_cast<typename dealii::internal::p4est::types<dim>::tree *>(&tree),
-                                                &p4est_cell)
-                        == false))
-        return;
-
-      bool p4est_has_children = (idx == -1);
-
-      // check if this quadrant is in the list
-      CellData  cell_data;
-      cell_data.quadrant = p4est_cell;
-      int pos = cell_data_search (cell_data, cells_to_compute);
-      if (pos != -1)
-        {
-          std::vector<CellData>  tmp;
-          // compute dof values
-          get_interpolated_dof_values (forest,
-                                       tree,
-                                       tree_index,
-                                       dealii_cell,
-                                       p4est_cell,
-                                       u,
-                                       cells_to_compute[pos].dof_values,
-                                       tmp);
-          // there is no new cell_data
-          // this process needs for computing this cell,
-          // store cell_data in the list of
-          // computed cells and erase this cell
-          // from the list of cells to compute
-          if (tmp.size () == 0)
-            {
-              cell_data_insert (cells_to_compute[pos], computed_cells);
-              cells_to_compute.erase (cells_to_compute.begin () + pos);
-            }
-          else
-            {
-              for (unsigned int i=0; i<tmp.size (); ++i)
-                cell_data_insert (tmp[i], new_needs);
-            }
-        }
-
-      // search recursively over this tree
-      if (p4est_has_children)
-        {
-          typename dealii::internal::p4est::types<dim>::quadrant
-          p4est_child[GeometryInfo<dim>::max_children_per_cell];
-
-          dealii::internal::p4est::init_quadrant_children<dim> (p4est_cell,
-                                                                p4est_child);
-
-          for (unsigned int c=0; c<GeometryInfo<dim>::max_children_per_cell; ++c)
-            {
-              compute_cells_in_tree_recursively (forest,
-                                                 tree,
-                                                 tree_index,
-                                                 dealii_cell->child(c),
-                                                 p4est_child[c],
-                                                 u,
-                                                 cells_to_compute,
-                                                 computed_cells,
-                                                 new_needs);
-            }
-        }
-    }
-
-
-
-    template <int dim,int spacedim>
-    void
-    ExtrapolateImplementation<dim,spacedim>::
-    send_cells (const std::vector<CellData>  &cells_to_send,
-                std::vector<CellData>        &received_cells) const
-    {
-      std::vector<std::vector<char> > sendbuffers (cells_to_send.size());
-      std::vector<std::vector<char> >::iterator buffer = sendbuffers.begin();
-      std::vector<MPI_Request> requests (cells_to_send.size());
-      std::vector<unsigned int> destinations;
-
-      // send data
-      unsigned int idx=0;
-      for (typename std::vector<CellData>::const_iterator it=cells_to_send.begin();
-           it!=cells_to_send.end();
-           ++it, ++idx)
-        {
-          destinations.push_back (it->receiver);
-
-          it->pack_data (*buffer);
-          MPI_Isend (&(*buffer)[0], buffer->size(),
-                     MPI_BYTE,
-                     it->receiver,
-                     123,
-                     communicator,
-                     &requests[idx]);
-        }
-
-      Assert(destinations.size()==cells_to_send.size(), ExcInternalError());
-
-      std::vector<unsigned int> senders
-        = Utilities::MPI::compute_point_to_point_communication_pattern(communicator, destinations);
-
-      // receive data
-      std::vector<char> receive;
-      CellData  cell_data;
-      for (unsigned int index=0; index<senders.size (); ++index)
-        {
-          MPI_Status status;
-          int len;
-          MPI_Probe(MPI_ANY_SOURCE, 123, communicator, &status);
-          MPI_Get_count(&status, MPI_BYTE, &len);
-          receive.resize (len);
-
-          char *buf = &receive[0];
-          MPI_Recv (buf, len, MPI_BYTE, status.MPI_SOURCE, status.MPI_TAG, communicator, &status);
-
-          cell_data.unpack_data (receive);
-
-          // this process has to send this
-          // cell back to the sender
-          // the receiver is the old sender
-          cell_data.receiver = status.MPI_SOURCE;
-
-          received_cells.push_back (cell_data);
-        }
-
-      if (requests.size () > 0)
-        MPI_Waitall(requests.size(), &requests[0], MPI_STATUSES_IGNORE);
-
-      // finally sort the list of cells
-      std::sort (received_cells.begin (), received_cells.end ());
-    }
-
-
-
-    template <int dim,int spacedim>
-    void
-    ExtrapolateImplementation<dim,spacedim>::
-    add_new_need (const typename dealii::internal::p4est::types<dim>::forest    &forest,
-                  const typename dealii::internal::p4est::types<dim>::locidx    &tree_index,
-                  const typename DoFHandler<dim,spacedim>::cell_iterator        &dealii_cell,
-                  const typename dealii::internal::p4est::types<dim>::quadrant  &p4est_cell,
-                  std::vector<CellData>                                         &new_needs) const
-    {
-      const FiniteElement<dim,spacedim> &fe            = dealii_cell->get_dof_handler().get_fe();
-      const unsigned int                 dofs_per_cell = fe.dofs_per_cell;
-
-      Assert (&fe != 0, ExcNotInitialized());
-
-      CellData  cell_data (dofs_per_cell);
-      cell_data.quadrant = p4est_cell;
-      cell_data.tree_index = tree_index;
-      cell_data.receiver = dealii::internal::p4est::functions<dim>::
-                           comm_find_owner (const_cast<typename dealii::internal::p4est::types<dim>::forest *> (&forest),
-                                            tree_index,
-                                            &p4est_cell,
-                                            dealii_cell->level_subdomain_id ());
-
-      cell_data_insert (cell_data, new_needs);
-    }
-
-
-
-    template <int dim,int spacedim>
-    int
-    ExtrapolateImplementation<dim,spacedim>::
-    cell_data_search (const CellData               &cell_data,
-                      const std::vector<CellData>  &cells_list)
-    {
-      typename std::vector<CellData>::const_iterator  bound
-        = std::lower_bound (cells_list.begin(), cells_list.end(), cell_data);
-
-      if ((bound != cells_list.end ())
-          &&
-          !(cell_data < *bound))
-        return (int)(bound - cells_list.begin ());
-
-      return (-1);
-    }
-
-
-
-    template <int dim,int spacedim>
-    void
-    ExtrapolateImplementation<dim,spacedim>::
-    cell_data_insert (const CellData         &cell_data,
-                      std::vector<CellData>  &cells_list)
-    {
-      typename std::vector<CellData>::iterator  bound
-        = std::lower_bound (cells_list.begin(), cells_list.end(), cell_data);
-
-      if ((bound == cells_list.end ())
-          ||
-          (cell_data < *bound))
-        cells_list.insert (bound, 1, cell_data);
-    }
-
-
-
-    template <int dim,int spacedim>
-    template <class InVector>
-    void
-    ExtrapolateImplementation<dim,spacedim>::
-    compute_all_non_local_data (const DoFHandler<dim,spacedim> &dof2,
-                                const InVector                 &u)
-    {
-      std::vector<CellData>  cells_we_need,
-          cells_to_compute,
-          received_cells,
-          received_needs,
-          new_needs,
-          computed_cells,
-          cells_to_send;
-
-      // Compute all the cells needed
-      // from other processes.
-      compute_needs (dof2, cells_we_need);
-
-      // Send the cells needed to there
-      // owners and receive a list of cells other
-      // processes need from us.
-      send_cells (cells_we_need, received_needs);
-
-      // The list of received needs can contain
-      // some cells more than once because different
-      // processes may need data from the same cell.
-      // To compute data only once, generate a vector
-      // with unique entries and distribute the computed
-      // data afterwards back to a vector with correct
-      // receivers.
-      // Computing cell_data can cause a need for
-      // data from some new cells.
-      // If a cell is computed, send it back to
-      // their senders, maybe receive new needs and
-      // compute again, do not wait that all cells
-      // are computed or all needs are collected.
-      // Otherwise we can run into a deadlock if
-      // a cell needed from another process,
-      // itself needs some data from us.
-      unsigned int ready = 0;
-      do
-        {
-          for (unsigned int i=0; i<received_needs.size(); ++i)
-            cell_data_insert (received_needs[i], cells_to_compute);
-
-          compute_cells (dof2, u, cells_to_compute, computed_cells, new_needs);
-
-          // if there are no cells to compute and no new needs, stop
-          ready = Utilities::MPI::sum (new_needs.size () + cells_to_compute.size (), communicator);
-
-          for (typename std::vector<CellData>::const_iterator comp=computed_cells.begin ();
-               comp != computed_cells.end ();
-               ++comp)
-            {
-              // store computed cells
-              cell_data_insert (*comp, available_cells);
-
-              // and generate a vector
-              // of computed cells with
-              // correct receivers
-              // then delete this received
-              // need from the list
-              typename std::vector<CellData>::iterator recv=received_needs.begin();
-              while (recv != received_needs.end())
-                {
-                  if (dealii::internal::p4est::quadrant_is_equal<dim>(recv->quadrant, comp->quadrant))
-                    {
-                      recv->dof_values = comp->dof_values;
-                      cells_to_send.push_back (*recv);
-                      received_needs.erase (recv);
-                      recv = received_needs.begin();
-                    }
-                  else
-                    ++recv;
-                }
-            }
-
-          send_cells (cells_to_send, received_cells);
-
-          // strore received cell_data
-          for (typename std::vector<CellData>::const_iterator recv=received_cells.begin ();
-               recv != received_cells.end ();
-               ++recv)
-            {
-              cell_data_insert (*recv, available_cells);
-            }
-
-          // finally send and receive new
-          // needs and start a new round
-          send_cells (new_needs, received_needs);
-        }
-      while (ready != 0);
-    }
-
-
-
-    template <int dim,int spacedim>
-    template <class InVector,class OutVector>
-    void
-    ExtrapolateImplementation<dim,spacedim>::
-    extrapolate_parallel (const InVector &u2_relevant,
-                          const DoFHandler<dim,spacedim> &dof2,
-                          OutVector &u2)
-    {
-      const parallel::distributed::Triangulation< dim, spacedim > *tr
-        = (dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>
-           (&dof2.get_tria()));
-
-      Assert (tr != 0,
-              ExcMessage ("Extrapolate in parallel only works for parallel distributed triangulations!"));
-
-      communicator = tr->get_communicator ();
-
-      compute_all_non_local_data (dof2, u2_relevant);
-
-      // exclude dofs on more refined ghosted cells
-      const FiniteElement<dim,spacedim> &fe  = dof2.get_fe();
-      const unsigned int dofs_per_face = fe.dofs_per_face;
-      if (dofs_per_face > 0)
-        {
-          const unsigned int dofs_per_cell = fe.dofs_per_cell;
-          std::vector<types::global_dof_index> indices (dofs_per_cell);
-          typename DoFHandler<dim,spacedim>::active_cell_iterator
-          cell=dof2.begin_active(),
-          endc=dof2.end();
-          for (; cell!=endc; ++cell)
-            if (cell->is_ghost())
-              {
-                cell->get_dof_indices(indices);
-                for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
-                  if (!cell->at_boundary(face))
-                    {
-                      const typename DoFHandler<dim,spacedim>::cell_iterator neighbor = cell->neighbor(face);
-                      if (neighbor->level() != cell->level())
-                        for (unsigned int i=0; i<dofs_per_face; ++i)
-                          {
-                            const types::global_dof_index index = indices[fe.face_to_cell_index(i,face)];;
-                            const bool index_stored
-                              = (dofs_on_refined_neighbors.find(index)!=dofs_on_refined_neighbors.end());
-                            const unsigned int level = index_stored ?
-                                                       std::max(cell->level(), dofs_on_refined_neighbors[index]) :
-                                                       cell->level();
-                            dofs_on_refined_neighbors[index] = level;
-                          }
-                    }
-              }
-        }
-
-      // after all dof values on
-      // not owned patch cells
-      // are computed, start
-      // the interpolation
-      u2 = 0;
-
-      std::queue<WorkPackage> queue;
-      {
-        typename DoFHandler<dim,spacedim>::cell_iterator
-        cell=dof2.begin(0),
-        endc=dof2.end(0);
-        for (; cell!=endc; ++cell)
-          {
-            if (dealii::internal::p4est::tree_exists_locally<dim> (tr->parallel_forest,
-                                                                   tr->coarse_cell_to_p4est_tree_permutation[cell->index()])
-                == false)
-              continue;
-
-            typename dealii::internal::p4est::types<dim>::quadrant p4est_coarse_cell;
-            const unsigned int tree_index = tr->coarse_cell_to_p4est_tree_permutation[cell->index()];
-            typename dealii::internal::p4est::types<dim>::tree *tree = tr->init_tree(cell->index());
-
-            dealii::internal::p4est::init_coarse_quadrant<dim>(p4est_coarse_cell);
-
-            const WorkPackage data (*tr->parallel_forest, *tree,tree_index,cell,p4est_coarse_cell);
-
-            queue.push(data);
-          }
-      }
-
-      while (!queue.empty())
-        {
-          const WorkPackage &data = queue.front();
-
-          const typename dealii::internal::p4est::types<dim>::forest    &forest = data.forest;
-          const typename dealii::internal::p4est::types<dim>::tree      &tree = data.tree;
-          const typename dealii::internal::p4est::types<dim>::locidx    &tree_index= data.tree_index;
-          const typename DoFHandler<dim,spacedim>::cell_iterator        &dealii_cell = data.dealii_cell;
-          const typename dealii::internal::p4est::types<dim>::quadrant  &p4est_cell = data.p4est_cell;
-
-          interpolate_recursively (forest, tree, tree_index, dealii_cell, p4est_cell, u2_relevant, u2);
-
-          // traverse recursively over this tree
-          if (dealii_cell->has_children())
-            {
-              Assert(dealii_cell->has_children(), ExcInternalError());
-              typename dealii::internal::p4est::types<dim>::quadrant
-              p4est_child[GeometryInfo<dim>::max_children_per_cell];
-
-              dealii::internal::p4est::init_quadrant_children<dim> (p4est_cell,
-                                                                    p4est_child);
-
-              for (unsigned int c=0; c<GeometryInfo<dim>::max_children_per_cell; ++c)
-                queue.push(WorkPackage(forest, tree, tree_index, dealii_cell->child(c), p4est_child[c]));
-            }
-          queue.pop();
-        }
-
-
-      u2.compress(VectorOperation::insert);
-    }
-#endif //DEAL_II_WITH_P4EST
-
-    namespace
-    {
-      template <int dim, class InVector, class OutVector, int spacedim>
-      void extrapolate_serial(const InVector &u3,
-                              const DoFHandler<dim,spacedim> &dof2,
-                              OutVector &u2)
-      {
-        const unsigned int dofs_per_cell  = dof2.get_fe().dofs_per_cell;
-        Vector<typename OutVector::value_type> dof_values(dofs_per_cell);
-
-        // then traverse grid bottom up
-        for (unsigned int level=0; level<dof2.get_triangulation().n_levels()-1; ++level)
-          {
-            typename DoFHandler<dim,spacedim>::cell_iterator cell=dof2.begin(level),
-                                                             endc=dof2.end(level);
-
-            for (; cell!=endc; ++cell)
-              if (!cell->active())
-                {
-                  // check whether this
-                  // cell has active
-                  // children
-                  bool active_children=false;
-                  for (unsigned int child_n=0; child_n<cell->n_children(); ++child_n)
-                    if (cell->child(child_n)->active())
-                      {
-                        active_children=true;
-                        break;
-                      }
-
-                  // if there are active
-                  // children, this process
-                  // has to work on this
-                  // cell. get the data
-                  // from the one vector
-                  // and set it on the
-                  // other
-                  if (active_children)
-                    {
-                      cell->get_interpolated_dof_values(u3, dof_values);
-                      cell->set_dof_values_by_interpolation(dof_values, u2);
-                    }
-                }
-          }
-      }
-
-
-
-#ifdef DEAL_II_WITH_P4EST
-      template <int dim, class InVector, class OutVector, int spacedim>
-      void extrapolate_parallel(const InVector &u2_relevant,
-                                const DoFHandler<dim,spacedim> &dof2,
-                                OutVector &u2)
-      {
-        if (dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>(&dof2.get_tria()) != 0)
-          {
-            ExtrapolateImplementation<dim,spacedim>  implementation;
-            implementation.extrapolate_parallel (u2_relevant, dof2, u2);
-          }
-        else
-          {
-            extrapolate_serial (u2_relevant, dof2, u2);
-          }
-      }
-#endif //DEAL_II_WITH_P4EST
-    }
-  }
-
-  template <int dim, class InVector, class OutVector, int spacedim>
-  void extrapolate(const DoFHandler<dim,spacedim> &dof1,
-                   const InVector &u1,
-                   const DoFHandler<dim,spacedim> &dof2,
-                   OutVector &u2)
-  {
-    ConstraintMatrix dummy;
-    dummy.close();
-    extrapolate(dof1, u1, dof2, dummy, u2);
-  }
-
-
-
-  template <int dim, class InVector, class OutVector, int spacedim>
-  void extrapolate(const DoFHandler<dim,spacedim> &dof1,
-                   const InVector &u1,
-                   const DoFHandler<dim,spacedim> &dof2,
-                   const ConstraintMatrix &constraints,
-                   OutVector &u2)
-  {
-    Assert(dof1.get_fe().n_components() == dof2.get_fe().n_components(),
-           ExcDimensionMismatch(dof1.get_fe().n_components(), dof2.get_fe().n_components()));
-    Assert(&dof1.get_triangulation()==&dof2.get_triangulation(), ExcTriangulationMismatch());
-    Assert(u1.size()==dof1.n_dofs(), ExcDimensionMismatch(u1.size(), dof1.n_dofs()));
-    Assert(u2.size()==dof2.n_dofs(), ExcDimensionMismatch(u2.size(), dof2.n_dofs()));
-
-    // make sure that each cell on the
-    // coarsest level is at least once
-    // refined, otherwise, these cells
-    // can't be treated and would
-    // generate a bogus result
-    {
-      typename DoFHandler<dim,spacedim>::cell_iterator cell = dof2.begin(0),
-                                                       endc = dof2.end(0);
-      for (; cell!=endc; ++cell)
-        Assert (cell->has_children(), ExcGridNotRefinedAtLeastOnce());
-    }
-
-    OutVector u3;
-    u3.reinit(u2);
-    interpolate(dof1, u1, dof2, constraints, u3);
-
-    internal::extrapolate_serial (u3, dof2, u2);
-
-    constraints.distribute(u2);
-  }
-
-
-#ifdef DEAL_II_WITH_P4EST
-  template <int dim, class VectorType, int spacedim>
-  void extrapolate_parallel(const DoFHandler<dim,spacedim> &dof1,
-                            const VectorType &u1,
-                            const DoFHandler<dim,spacedim> &dof2,
-                            const ConstraintMatrix &constraints2,
-                            VectorType &u2)
-  {
-    IndexSet  dof2_locally_owned_dofs = dof2.locally_owned_dofs();
-    IndexSet  dof2_locally_relevant_dofs;
-    DoFTools::extract_locally_relevant_dofs (dof2,
-                                             dof2_locally_relevant_dofs);
-
-    VectorType u3 (dof2_locally_owned_dofs, u1.get_mpi_communicator ());
-    interpolate (dof1, u1, dof2, constraints2, u3);
-    VectorType u3_relevant (dof2_locally_owned_dofs,
-                            dof2_locally_relevant_dofs,
-                            u1.get_mpi_communicator ());
-    u3_relevant = u3;
-
-    internal::extrapolate_parallel (u3_relevant, dof2, u2);
-
-    constraints2.distribute(u2);
-  }
-
-
-
-  template <int dim, int spacedim, typename Number>
-  void extrapolate_parallel(const DoFHandler<dim,spacedim> &dof1,
-                            const LinearAlgebra::distributed::Vector<Number> &u1,
-                            const DoFHandler<dim,spacedim> &dof2,
-                            const ConstraintMatrix &constraints2,
-                            LinearAlgebra::distributed::Vector<Number> &u2)
-  {
-    IndexSet  dof2_locally_owned_dofs = dof2.locally_owned_dofs();
-    IndexSet  dof2_locally_relevant_dofs;
-    DoFTools::extract_locally_relevant_dofs (dof2,
-                                             dof2_locally_relevant_dofs);
-
-    LinearAlgebra::distributed::Vector<Number> u3(dof2_locally_owned_dofs,
-                                                  dof2_locally_relevant_dofs,
-                                                  u1.get_mpi_communicator ());
-    interpolate (dof1, u1, dof2, constraints2, u3);
-    u3.update_ghost_values();
-
-    internal::extrapolate_parallel (u3, dof2, u2);
-
-    constraints2.distribute(u2);
-  }
-#endif //DEAL_II_WITH_P4EST
-
 } // end of namespace FETools
 
 
old mode 100644 (file)
new mode 100755 (executable)
index 24fcd7d..d051fdf
@@ -1,6 +1,6 @@
 // ---------------------------------------------------------------------
 //
-// Copyright (C) 1998 - 2014 by the deal.II authors
+// Copyright (C) 1998 - 2016 by the deal.II authors
 //
 // This file is part of the deal.II library.
 //
@@ -97,67 +97,6 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS
     void project_dg<deal_II_dimension>
     (const DoFHandler<deal_II_dimension> &, const VEC &,
      const DoFHandler<deal_II_dimension> &, VEC &);
-
-    template
-    void extrapolate<deal_II_dimension>
-    (const DoFHandler<deal_II_dimension> &, const VEC &,
-     const DoFHandler<deal_II_dimension> &, VEC &);
-
-    template
-    void extrapolate<deal_II_dimension>
-    (const DoFHandler<deal_II_dimension> &, const VEC &,
-     const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
-     VEC &);
-#endif
-    \}
-}
-
-for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS)
-{
-#ifdef DEAL_II_WITH_P4EST
-    namespace FETools
-    \{
-#if deal_II_dimension == deal_II_space_dimension
-
-#ifdef DEAL_II_WITH_PETSC
-    template
-    void extrapolate_parallel<deal_II_dimension, PETScWrappers::MPI::Vector, deal_II_space_dimension>
-    (const DoFHandler<deal_II_dimension,deal_II_space_dimension> &,
-     const PETScWrappers::MPI::Vector &,
-     const DoFHandler<deal_II_dimension,deal_II_space_dimension> &,
-     const ConstraintMatrix &,
-     PETScWrappers::MPI::Vector &);
-#endif
-
-#ifdef DEAL_II_WITH_TRILINOS
-    template
-    void extrapolate_parallel <deal_II_dimension, TrilinosWrappers::MPI::Vector, deal_II_space_dimension>
-    (const DoFHandler<deal_II_dimension,deal_II_space_dimension> &,
-     const TrilinosWrappers::MPI::Vector &,
-     const DoFHandler<deal_II_dimension,deal_II_space_dimension> &,
-     const ConstraintMatrix &,
-     TrilinosWrappers::MPI::Vector &);
-#endif
-
-#endif
-    \}
-#endif //DEAL_II_WITH_P4EST
-}
-
-for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS; scalar : REAL_SCALARS)
-{
-#ifdef DEAL_II_WITH_P4EST
-    namespace FETools
-    \{
-#if deal_II_dimension == deal_II_space_dimension
-    template
-    void extrapolate_parallel <deal_II_dimension, deal_II_space_dimension, scalar>
-    (const DoFHandler<deal_II_dimension,deal_II_space_dimension> &,
-     const LinearAlgebra::distributed::Vector<scalar> &,
-     const DoFHandler<deal_II_dimension,deal_II_space_dimension> &,
-     const ConstraintMatrix &,
-     LinearAlgebra::distributed::Vector<scalar> &);
 #endif
     \}
-#endif //DEAL_II_WITH_P4EST
 }
diff --git a/tests/mpi/fe_tools_extrapolate_03.mpirun=10.output b/tests/mpi/fe_tools_extrapolate_03.mpirun=10.output
new file mode 100644 (file)
index 0000000..f406bba
--- /dev/null
@@ -0,0 +1,115 @@
+
+DEAL:0::Checking Q1 against Q1 in 2d:
+DEAL:0::23.00000000 3.944933460 1.000000000
+DEAL:0::23.00000000 3.944933460 1.000000000
+DEAL:0::global_error_before: 0
+DEAL:0::global_error_after: 0
+DEAL:0::OK
+DEAL:0::Checking Q1 against Q1 in 3d:
+DEAL:0::162.3750000 12.71625436 1.500000000
+DEAL:0::162.3750000 12.71625436 1.500000000
+DEAL:0::global_error_before: 0
+DEAL:0::global_error_after: 0
+DEAL:0::OK
+DEAL:0::Checking Q1 against Q2 in 2d:
+DEAL:0::9.125000000 1.682283604 0.5000000000
+DEAL:0::26.88281250 2.700093063 0.5000000000
+DEAL:0::global_error_before: 0.01934481330
+DEAL:0::global_error_after: 0.001208379635
+DEAL:0::OK
+DEAL:0::Checking Q1 against Q2 in 3d:
+DEAL:0::66.28125000 5.415995551 0.7500000000
+DEAL:0::342.1987305 11.49092805 0.7500000000
+DEAL:0::global_error_before: 0.03152756518
+DEAL:0::global_error_after: 0.001607487444
+DEAL:0::OK
+DEAL:0::Checking Q2 against Q1 in 2d:
+DEAL:0::26.84375000 2.699035020 0.5000000000
+DEAL:0::10.25000000 1.831793861 0.5000000000
+DEAL:0::global_error_before: 0
+DEAL:0::global_error_after: 0.07284222390
+DEAL:0::OK
+DEAL:0::Checking Q2 against Q1 in 3d:
+DEAL:0::341.5859375 11.47797042 0.7500000000
+DEAL:0::75.09375000 5.992915740 0.7500000000
+DEAL:0::global_error_before: 0
+DEAL:0::global_error_after: 0.1190491594
+DEAL:0::OK
+DEAL:0::Checking Q2 against Q2 in 2d:
+DEAL:0::26.84375000 2.699035020 0.5000000000
+DEAL:0::26.84375000 2.699035020 0.5000000000
+DEAL:0::global_error_before: 0
+DEAL:0::global_error_after: 0
+DEAL:0::OK
+DEAL:0::Checking Q2 against Q2 in 3d:
+DEAL:0::341.5859375 11.47797042 0.7500000000
+DEAL:0::341.5859375 11.47797042 0.7500000000
+DEAL:0::global_error_before: 0
+DEAL:0::global_error_after: 0
+DEAL:0::OK
+DEAL:0::Checking DGQ0 against DGQ0 in 2d:
+DEAL:0::56.00000000 10.58300524 2.000000000
+DEAL:0::56.00000000 10.58300524 2.000000000
+DEAL:0::global_error_before: 0
+DEAL:0::global_error_after: 0
+DEAL:0::OK
+DEAL:0::Checking DGQ0 against DGQ0 in 3d:
+DEAL:0::276.0000000 28.77498914 3.000000000
+DEAL:0::276.0000000 28.77498914 3.000000000
+DEAL:0::global_error_before: 0
+DEAL:0::global_error_after: 0
+DEAL:0::OK
+DEAL:0::Checking DGQ0 against DGQ1 in 2d:
+DEAL:0::14.00000000 2.828427125 0.8750000000
+DEAL:0::56.00000000 5.678908346 0.9375000000
+DEAL:0::global_error_before: 0.09199750902
+DEAL:0::global_error_after: 0.04599875451
+DEAL:0::OK
+DEAL:0::Checking DGQ0 against DGQ1 in 3d:
+DEAL:0::72.50000000 7.854139036 1.312500000
+DEAL:0::580.6250000 22.31268054 1.406250000
+DEAL:0::global_error_before: 0.1220351512
+DEAL:0::global_error_after: 0.06040494959
+DEAL:0::OK
+DEAL:0::Checking DGQ1 against DGQ0 in 2d:
+DEAL:0::56.00000000 5.830951895 1.000000000
+DEAL:0::14.00000000 2.738612788 0.7500000000
+DEAL:0::global_error_before: 0
+DEAL:0::global_error_after: 0.1839950180
+DEAL:0::OK
+DEAL:0::Checking DGQ1 against DGQ0 in 3d:
+DEAL:0::580.0000000 22.78157150 1.500000000
+DEAL:0::73.00000000 7.697402159 1.125000000
+DEAL:0::global_error_before: 0
+DEAL:0::global_error_after: 0.2440703024
+DEAL:0::OK
+DEAL:0::Checking DGQ1 against DGQ1 in 2d:
+DEAL:0::56.00000000 5.830951895 1.000000000
+DEAL:0::56.00000000 5.830951895 1.000000000
+DEAL:0::global_error_before: 0
+DEAL:0::global_error_after: 0
+DEAL:0::OK
+DEAL:0::Checking DGQ1 against DGQ1 in 3d:
+DEAL:0::580.0000000 22.78157150 1.500000000
+DEAL:0::580.0000000 22.78157150 1.500000000
+DEAL:0::global_error_before: 0
+DEAL:0::global_error_after: 0
+DEAL:0::OK
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
old mode 100644 (file)
new mode 100755 (executable)
index 353fea4..5cd255f
@@ -73,7 +73,8 @@ parallel::distributed::Triangulation<dim> *make_tria ()
       cell = tria->begin_active();
       cell->set_refine_flag();
       ++cell;
-      cell->set_refine_flag();
+      if (cell != tria->end())
+        cell->set_refine_flag();
 
       tria->execute_coarsening_and_refinement ();
     }
@@ -98,8 +99,6 @@ template <unsigned int dim, typename VectorType>
 void
 output_vector (const VectorType &v, const std::string &output_name, const DoFHandler<dim> &dof_handler)
 {
-  //v.print(deallog.get_file_stream());
-
   DataOut<dim> data_out;
   data_out.attach_dof_handler (dof_handler);
   data_out.add_data_vector (v, output_name, DataOut<dim>::type_dof_data);

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.