]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Solution transfer for hp elements.
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 10 Mar 2012 11:22:23 +0000 (11:22 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 10 Mar 2012 11:22:23 +0000 (11:22 +0000)
git-svn-id: https://svn.dealii.org/trunk@25237 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/include/deal.II/numerics/solution_transfer.h
deal.II/source/numerics/solution_transfer.cc
tests/hp/solution_transfer.cc [new file with mode: 0644]
tests/hp/solution_transfer/cmp/generic [new file with mode: 0644]

index 03e8d52645ef03f7b4e06b7e3355589d100481f4..c16b1b6e5f9745ab1dada629ba84c5737aac72e4 100644 (file)
@@ -152,6 +152,12 @@ enabled due to a missing include file in file
 <h3>Specific improvements</h3>
 
 <ol>
+<li> Extended:
+SolutionTransfer is now also able to transfer solutions between hp::DoFHandler where
+the finite element index changes during refinement.
+<br>
+(Katharina Kormann, Martin Kronbichler, 2012/03/10)
+
 <li> Changed:
 A new method to determine an initial guess for the Newton method was coded
 in MappingQ::transform_real_to_unit_cell.
index a214daaeb603abd51f2272a3109f18d59539c65d..1e9a1c31f6a135a4cd64b676cccd90a7b1fefe4d 100644 (file)
@@ -273,7 +273,7 @@ class SolutionTransfer
                                     /**
                                      * This function
                                      * interpolates the discrete functions
-                                     * that are stored in @p all_out onto
+                                     * that are stored in @p all_in onto
                                      * the refined and/or coarsenend grid.
                                      * It assumes the vectors in @p all_in
                                      * denote the same vectors
@@ -423,17 +423,29 @@ class SolutionTransfer
                                      * <tt>vector<double> dof_values</tt> (if the
                                      * children of this cell will be deleted)
                                      * is needed, hence one @p user_pointer should
-                                     * be sufficient, but to allow some errorchecks
+                                     * be sufficient, but to allow some error checks
                                      * and to preserve the user from making
                                      * user errors the @p user_pointer will be
                                      * 'multiplied' by this structure.
                                      */
     struct Pointerstruct {
-       Pointerstruct();
+      Pointerstruct() : indices_ptr(0), dof_values_ptr(0), active_fe_index(0) {};
+      Pointerstruct(std::vector<unsigned int> *indices_ptr_in,
+                   const unsigned int active_fe_index_in = 0)
+       : 
+       indices_ptr(indices_ptr_in),
+       dof_values_ptr (0),
+       active_fe_index(active_fe_index_in) {};
+      Pointerstruct(std::vector<Vector<typename VECTOR::value_type> > *dof_values_ptr_in,
+                   const unsigned int active_fe_index_in = 0) :
+       indices_ptr (0),
+       dof_values_ptr(dof_values_ptr_in),
+       active_fe_index(active_fe_index_in) {};
        std::size_t memory_consumption () const;
 
        std::vector<unsigned int>    *indices_ptr;
        std::vector<Vector<typename VECTOR::value_type> > *dof_values_ptr;
+        unsigned int active_fe_index;
     };
 
                                     /**
index eeca9f945eefe939be2a0dfacb06108476b9171b..a88ae40eda6ce0666563afb7437cdd4098495b7a 100644 (file)
@@ -40,6 +40,7 @@ SolutionTransfer<dim, VECTOR, DH>::SolutionTransfer(const DH &dof)
 {}
 
 
+
 template<int dim, typename VECTOR, class DH>
 SolutionTransfer<dim, VECTOR, DH>::~SolutionTransfer()
 {
@@ -47,13 +48,6 @@ SolutionTransfer<dim, VECTOR, DH>::~SolutionTransfer()
 }
 
 
-template<int dim, typename VECTOR, class DH>
-SolutionTransfer<dim, VECTOR, DH>::Pointerstruct::Pointerstruct():
-               indices_ptr(0),
-               dof_values_ptr(0)
-{}
-
-
 
 template<int dim, typename VECTOR, class DH>
 void SolutionTransfer<dim, VECTOR, DH>::clear ()
@@ -66,6 +60,7 @@ void SolutionTransfer<dim, VECTOR, DH>::clear ()
 }
 
 
+
 template<int dim, typename VECTOR, class DH>
 void SolutionTransfer<dim, VECTOR, DH>::prepare_for_pure_refinement()
 {
@@ -101,6 +96,7 @@ void SolutionTransfer<dim, VECTOR, DH>::prepare_for_pure_refinement()
 }
 
 
+
 template<int dim, typename VECTOR, class DH>
 void
 SolutionTransfer<dim, VECTOR, DH>::refine_interpolate(const VECTOR &in,
@@ -155,6 +151,51 @@ SolutionTransfer<dim, VECTOR, DH>::refine_interpolate(const VECTOR &in,
 }
 
 
+
+namespace internal
+{
+  template <typename DH>
+  void extract_interpolation_matrices (const DH&,
+                                      Table<2,FullMatrix<double> > &)
+  {}
+
+  template <int dim, int spacedim>
+  void extract_interpolation_matrices (const dealii::hp::DoFHandler<dim,spacedim> &dof,
+                                      Table<2,FullMatrix<double> > &matrices)
+  {
+    const dealii::hp::FECollection<dim,spacedim> &fe = dof.get_fe();
+    matrices.reinit (fe.size(), fe.size());
+    for (unsigned int i=0; i<fe.size(); ++i)
+      for (unsigned int j=0; j<fe.size(); ++j)
+       if (i != j)
+         {
+           matrices(i,j).reinit (fe[i].dofs_per_cell, fe[j].dofs_per_cell);
+           fe[i].get_interpolation_matrix (fe[j], matrices(i,j));
+         }
+  }
+
+
+  template <int dim, int spacedim>
+  void restriction_additive (const FiniteElement<dim,spacedim> &,
+                            std::vector<std::vector<bool> > &)
+  {}
+
+  template <int dim, int spacedim>
+  void restriction_additive (const dealii::hp::FECollection<dim,spacedim> &fe,
+                            std::vector<std::vector<bool> > &restriction_is_additive)
+  {
+    restriction_is_additive.resize (fe.size());
+    for (unsigned int f=0; f<fe.size(); ++f)
+      {
+       restriction_is_additive[f].resize (fe[f].dofs_per_cell);
+       for (unsigned int i=0; i<fe[f].dofs_per_cell; ++i)
+         restriction_is_additive[f][i] = fe[f].restriction_is_additive(i);
+      }
+  }
+}
+
+
+
 template<int dim, typename VECTOR, class DH>
 void
 SolutionTransfer<dim, VECTOR, DH>::
@@ -218,6 +259,13 @@ prepare_for_coarsening_and_refinement(const std::vector<VECTOR> &all_in)
      std::vector<Vector<typename VECTOR::value_type> > (in_size))
     .swap(dof_values_on_cell);
 
+  typename VECTOR::value_type zero_val = typename VECTOR::value_type();
+  Table<2,FullMatrix<double> > interpolation_hp;
+  std::vector<std::vector<bool> > restriction_is_additive;
+  internal::extract_interpolation_matrices (*dof_handler, interpolation_hp);
+  internal::restriction_additive (dof_handler->get_fe(), restriction_is_additive);
+  Vector<typename VECTOR::value_type> tmp, tmp2;
+
                                   // we need counters for
                                   // the 'to_stay_or_refine' cells 'n_sr' and
                                   // the 'coarsen_fathers' cells 'n_cf',
@@ -234,7 +282,8 @@ prepare_for_coarsening_and_refinement(const std::vector<VECTOR> &all_in)
                                           // dof indices and later
                                           // interpolating to the children
          cell->get_dof_indices(indices_on_cell[n_sr]);
-         cell_map[std::make_pair(cell->level(), cell->index())].indices_ptr=&indices_on_cell[n_sr];
+         cell_map[std::make_pair(cell->level(), cell->index())]
+           = Pointerstruct(&indices_on_cell[n_sr],cell->active_fe_index());
          ++n_sr;
        }
       else if (cell->has_children() && cell->child(0)->coarsen_flag_set())
@@ -249,17 +298,76 @@ prepare_for_coarsening_and_refinement(const std::vector<VECTOR> &all_in)
          const unsigned int dofs_per_cell=cell->get_fe().dofs_per_cell;
 
          std::vector<Vector<typename VECTOR::value_type> >(in_size,
-                                                               Vector<typename VECTOR::value_type>(dofs_per_cell))
+                                                           Vector<typename VECTOR::value_type>(dofs_per_cell))
            .swap(dof_values_on_cell[n_cf]);
 
+         unsigned int fe_index = cell->active_fe_index();
+         unsigned int most_general_child = 0;
+         bool different_elements = false;
+         for (unsigned int child=0; child<cell->n_children(); ++child)
+           {
+             if (cell->child(child)->active_fe_index() != fe_index)
+               different_elements = true;
+                               // take FE index from the child with most
+                               // degrees of freedom locally
+             if (cell->child(child)->get_fe().dofs_per_cell >
+                 cell->child(most_general_child)->get_fe().dofs_per_cell)
+               most_general_child = child;
+           }
+         if (different_elements == true)
+           fe_index = cell->child(most_general_child)->active_fe_index();
+
          for (unsigned int j=0; j<in_size; ++j)
            {
                                               // store the data of each of
                                               // the input vectors
-             cell->get_interpolated_dof_values(all_in[j],
-                                               dof_values_on_cell[n_cf][j]);
+             if (different_elements == false)
+               cell->get_interpolated_dof_values(all_in[j],
+                                                 dof_values_on_cell[n_cf][j]);
+             else
+               {
+                               // if we have different elements, first
+                               // interpolate the children's contribution to
+                               // the most general FE on the child
+                               // level. Then we manually write the
+                               // interpolation operation to the coarser
+                               // level
+                 dof_values_on_cell[n_cf][j].reinit (cell->child(most_general_child)->get_fe().dofs_per_cell);
+                 const unsigned int fe_ind_general = 
+                   cell->child(most_general_child)->active_fe_index();
+                 for (unsigned int child=0; child<cell->n_children(); ++child)
+                   {
+                     tmp.reinit (cell->child(child)->get_fe().dofs_per_cell,
+                                 true);
+                     cell->child(child)->get_dof_values (all_in[j],
+                                                         tmp);
+                     const unsigned int child_ind =
+                       cell->child(child)->active_fe_index();
+                     if (child_ind != fe_ind_general)
+                       {
+                         tmp2.reinit (cell->child(most_general_child)->get_fe().dofs_per_cell,
+                                      true);
+                         interpolation_hp (fe_ind_general, child_ind).vmult (tmp2, tmp);
+                       }
+                     else
+                       tmp2.swap (tmp);
+
+                               // now do the interpolation operation
+                     const unsigned int dofs_per_cell = tmp2.size();
+                     tmp.reinit (dofs_per_cell, true);
+                     cell->child(most_general_child)->get_fe().
+                       get_restriction_matrix(child, cell->refinement_case()).vmult (tmp, tmp2);
+                     for (unsigned int i=0; i<dofs_per_cell; ++i)
+                       if (restriction_is_additive[fe_ind_general][i])
+                         dof_values_on_cell[n_cf][j](i) += tmp(i);
+                       else
+                         if (tmp(i) != zero_val)
+                           dof_values_on_cell[n_cf][j](i) = tmp(i);
+                   }                 
+               }
            }
-         cell_map[std::make_pair(cell->level(), cell->index())].dof_values_ptr=&dof_values_on_cell[n_cf];
+         cell_map[std::make_pair(cell->level(), cell->index())]
+           = Pointerstruct(&dof_values_on_cell[n_cf], fe_index);
          ++n_cf;
        }
     }
@@ -270,6 +378,7 @@ prepare_for_coarsening_and_refinement(const std::vector<VECTOR> &all_in)
 }
 
 
+
 template<int dim, typename VECTOR, class DH>
 void
 SolutionTransfer<dim, VECTOR, DH>::prepare_for_coarsening_and_refinement(const VECTOR &in)
@@ -279,6 +388,7 @@ SolutionTransfer<dim, VECTOR, DH>::prepare_for_coarsening_and_refinement(const V
 }
 
 
+
 template<int dim, typename VECTOR, class DH>
 void SolutionTransfer<dim, VECTOR, DH>::
 interpolate (const std::vector<VECTOR> &all_in,
@@ -306,6 +416,10 @@ interpolate (const std::vector<VECTOR> &all_in,
     pointerstruct,
     cell_map_end=cell_map.end();
 
+  Table<2,FullMatrix<double> > interpolation_hp;
+  internal::extract_interpolation_matrices (*dof_handler, interpolation_hp);
+  Vector<typename VECTOR::value_type> tmp, tmp2;
+
   typename DH::cell_iterator cell = dof_handler->begin(),
                             endc = dof_handler->end();
   for (; cell!=endc; ++cell)
@@ -328,28 +442,91 @@ interpolate (const std::vector<VECTOR> &all_in,
            {
              Assert (valuesptr == 0,
                      ExcInternalError());
-                                              // make sure that the size of the
-                                              // stored indices is the same as
-                                              // dofs_per_cell. this is kind of
-                                              // a test if we use the same fe
-                                              // in the hp case. to really do
-                                              // that test we would have to
-                                              // store the fe_index of all
-                                              // cells
-             Assert(dofs_per_cell==(*indexptr).size(),
-                    ExcNumberOfDoFsPerCellHasChanged());
+
                                               // get the values of
                                               // each of the input
                                               // data vectors on this
                                               // cell and prolong it
                                               // to its children
+             unsigned int in_size = indexptr->size();
              local_values.reinit(dofs_per_cell, true);
-             for (unsigned int j=0; j<size; ++j)
+             for (unsigned int j=0; j<size; ++j)
                {
-                 for (unsigned int i=0; i<dofs_per_cell; ++i)
-                   local_values(i)=all_in[j]((*indexptr)[i]);
-                 cell->set_dof_values_by_interpolation(local_values,
-                                                       all_out[j]);
+                               // check the FE index of the new element. if
+                               // we have children and not all of the
+                               // children have the same FE index, need to
+                               // manually implement
+                               // set_dof_values_by_interpolation
+                 unsigned int new_fe_index = cell->active_fe_index();
+                 bool different_elements = false;
+                 if (cell->has_children())
+                   {
+                     new_fe_index = cell->child(0)->active_fe_index();
+                     for (unsigned int child=1; child<cell->n_children(); ++child)
+                       if (cell->child(child)->active_fe_index() != new_fe_index)
+                         {
+                           different_elements = true;
+                           break;
+                         }
+                   }
+                 if (different_elements == true ||
+                     new_fe_index != pointerstruct->second.active_fe_index)
+                   {
+                     const unsigned int old_index =
+                       pointerstruct->second.active_fe_index;
+                     tmp.reinit (in_size, true);
+                     for (unsigned int i=0; i<in_size; ++i)
+                       tmp(i)=all_in[j]((*indexptr)[i]);
+                     if (different_elements == false)
+                       {
+                         AssertDimension (tmp.size(),
+                                          interpolation_hp(new_fe_index,old_index).n());
+                         local_values.reinit (cell->has_children() ?
+                                              cell->child(0)->get_fe().dofs_per_cell
+                                              : cell->get_fe().dofs_per_cell, true);
+                         AssertDimension (local_values.size(),
+                                          interpolation_hp(new_fe_index,old_index).m());
+                               // simple case where all children have the
+                               // same FE index: just interpolate to their FE
+                               // first and then use the standard routines
+                         interpolation_hp(new_fe_index,old_index).vmult (local_values, tmp);
+                       }
+
+                     if (cell->has_children() == false)
+                       cell->set_dof_values (local_values, all_out[j]);
+                     else
+                     for (unsigned int child=0; child<cell->n_children(); ++child)
+                       {
+                         if (different_elements == true)
+                           {
+                             const unsigned int c_index =
+                               cell->child(child)->active_fe_index();
+                             if (c_index != old_index)
+                               {
+                                 AssertDimension (tmp.size(),
+                                                  interpolation_hp(c_index,old_index).n());
+                                 local_values.reinit(cell->child(child)->get_fe().dofs_per_cell, true);
+                                 AssertDimension (local_values.size(),
+                                                  interpolation_hp(c_index,old_index).m());
+                                 interpolation_hp(c_index,old_index).vmult (local_values, tmp);
+                               }
+                             else
+                               local_values = tmp;
+                           }
+                         tmp2.reinit (cell->child(child)->get_fe().dofs_per_cell, true);
+                         cell->child(child)->get_fe().get_prolongation_matrix(child, cell->refinement_case())
+                           .vmult (tmp2, local_values);
+                         cell->child(child)->set_dof_values(tmp2,all_out[j]);
+                       }
+                   }
+                 else
+                   {
+                     AssertDimension (dofs_per_cell, indexptr->size());
+                     for (unsigned int i=0; i<dofs_per_cell; ++i)
+                       local_values(i)=all_in[j]((*indexptr)[i]);
+                     cell->set_dof_values_by_interpolation(local_values,
+                                                           all_out[j]);
+                   }
                }
            }
                                           // children of cell were
@@ -379,11 +556,25 @@ interpolate (const std::vector<VECTOR> &all_in,
                                                   // test we would have to
                                                   // store the fe_index of all
                                                   // cells
-                 Assert(dofs_per_cell==(*valuesptr)[j].size(),
-                        ExcNumberOfDoFsPerCellHasChanged());
+                 const Vector<typename VECTOR::value_type>* data = 0; 
+                 const unsigned int active_fe_index = cell->active_fe_index();
+                 if (active_fe_index != pointerstruct->second.active_fe_index)
+                   {
+                     const unsigned int old_index = pointerstruct->second.active_fe_index;
+                     tmp.reinit (dofs_per_cell, true);
+                     AssertDimension ((*valuesptr)[j].size(),
+                                      interpolation_hp(active_fe_index,old_index).n());
+                     AssertDimension (tmp.size(),
+                                      interpolation_hp(active_fe_index,old_index).m());
+                     interpolation_hp(active_fe_index,old_index).vmult (tmp, (*valuesptr)[j]);
+                     data = &tmp;
+                   }
+                 else
+                   data = &(*valuesptr)[j];
+                 
 
                  for (unsigned int i=0; i<dofs_per_cell; ++i)
-                   all_out[j](dofs[i])=(*valuesptr)[j](i);
+                   all_out[j](dofs[i])=(*data)(i);
                }
            }
                                           // undefined status
diff --git a/tests/hp/solution_transfer.cc b/tests/hp/solution_transfer.cc
new file mode 100644 (file)
index 0000000..01d9fb5
--- /dev/null
@@ -0,0 +1,246 @@
+//----------------------------  solution_transfer.cc  ---------------------------
+//    fe_data_test.cc,v 1.14 2003/11/28 11:52:35 guido Exp
+//    Version:
+//
+//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2005, 2006, 2007, 2008, 2012 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  solution_transfer.cc  ---------------------------
+
+
+#include "../tests.h"
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/numerics/vectors.h>
+#include <deal.II/numerics/solution_transfer.h>
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/mapping_q1.h>
+#include <deal.II/hp/fe_values.h>
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/hp/fe_collection.h>
+#include <fstream>
+#include <iostream>
+#include <vector>
+
+// a linear function that should be transferred exactly with Q1 and Q2
+// elements
+template<int dim>
+class MyFunction : public Function<dim>
+{
+  public:
+    MyFunction () : Function<dim>() {};
+
+    virtual double value (const Point<dim>   &p,
+                         const unsigned int) const
+      {
+       double f=0.25 + 2 * p[0];
+       if (dim>1)
+         f+=0.772 * p[1];
+       if (dim>2)
+         f-=3.112 * p[2];
+       return f;
+      };
+};
+
+
+template <int dim>
+void transfer(std::ostream &out)
+{
+  MyFunction<dim> function;
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria);
+  tria.refine_global(5-dim);
+  const unsigned int max_degree = 6-dim;
+  hp::FECollection<dim> fe_q;
+  hp::FECollection<dim> fe_dgq;
+  for (unsigned int deg=1; deg<=max_degree; ++deg)
+    {
+      fe_q.push_back (FE_Q<dim>(deg));
+      fe_dgq.push_back (FE_DGQ<dim>(deg));
+    }
+  hp::DoFHandler<dim> q_dof_handler(tria);
+  hp::DoFHandler<dim> dgq_dof_handler(tria);
+  Vector<double> q_solution;
+  Vector<double> dgq_solution;
+  MappingQ1<dim> mapping;
+
+                               // refine a few cells
+  typename Triangulation<dim>::active_cell_iterator cell=tria.begin_active(),
+                                                   endc=tria.end();
+  ++cell;
+  ++cell;
+  for (; cell!=endc; ++cell)
+    cell->set_refine_flag();
+  tria.prepare_coarsening_and_refinement();
+  tria.execute_coarsening_and_refinement();
+
+                               // randomly assign FE orders
+  unsigned int counter = 0;
+  {
+    typename hp::DoFHandler<dim>::active_cell_iterator
+      cell=q_dof_handler.begin_active(),
+      celldg=dgq_dof_handler.begin_active(),
+      endc=q_dof_handler.end();
+    for (; cell!=endc; ++cell, ++celldg, ++counter)
+      {
+       if (counter < 15)
+         cell->set_active_fe_index(1);
+       else
+         cell->set_active_fe_index(rand()%max_degree);
+       if (counter < 15)
+         celldg->set_active_fe_index(1);
+       else
+         celldg->set_active_fe_index(rand()%max_degree);
+      }
+  }
+
+  q_dof_handler.distribute_dofs (fe_q);
+  q_solution.reinit(q_dof_handler.n_dofs());
+
+  dgq_dof_handler.distribute_dofs (fe_dgq);
+  dgq_solution.reinit(dgq_dof_handler.n_dofs());
+
+  ConstraintMatrix cm;
+  cm.close();
+  VectorTools::interpolate (mapping, q_dof_handler, function, q_solution);
+  VectorTools::interpolate (mapping, dgq_dof_handler, function, dgq_solution);
+
+  SolutionTransfer<dim,Vector<double>,hp::DoFHandler<dim> > q_soltrans(q_dof_handler);
+  SolutionTransfer<dim,Vector<double>,hp::DoFHandler<dim> > dgq_soltrans(dgq_dof_handler);
+
+
+                               // test b): do some coarsening and
+                               // refinement
+  q_soltrans.clear();
+  dgq_soltrans.clear();
+
+  counter = 0;
+  cell=tria.begin_active();
+  endc=tria.end();
+  for (; cell!=endc; ++cell, ++counter)
+    {
+      if (counter > 120)
+       cell->set_coarsen_flag();
+      else if (rand() % 3 == 0)
+       cell->set_refine_flag();
+      else if (rand() % 3 == 3)
+       cell->set_coarsen_flag();
+    }
+
+  Vector<double> q_old_solution=q_solution,
+    dgq_old_solution=dgq_solution;
+  tria.prepare_coarsening_and_refinement();
+  q_soltrans.prepare_for_coarsening_and_refinement(q_old_solution);
+  dgq_soltrans.prepare_for_coarsening_and_refinement(dgq_old_solution);
+  tria.execute_coarsening_and_refinement();
+
+  counter = 0;
+  {
+    typename hp::DoFHandler<dim>::active_cell_iterator
+      cell = q_dof_handler.begin_active(),
+      celldg = dgq_dof_handler.begin_active(),
+      endc = q_dof_handler.end();
+    for (; cell!=endc; ++cell, ++celldg, ++counter)
+      {
+       if (counter > 20 && counter < 90)
+         cell->set_active_fe_index(0);
+       else
+         cell->set_active_fe_index(rand()%max_degree);
+       if (counter > 20 && counter < 90)
+         celldg->set_active_fe_index(0);
+       else
+         celldg->set_active_fe_index(rand()%max_degree);
+      }
+  }
+
+  q_dof_handler.distribute_dofs (fe_q);
+  dgq_dof_handler.distribute_dofs (fe_dgq);
+  q_solution.reinit(q_dof_handler.n_dofs());
+  dgq_solution.reinit(dgq_dof_handler.n_dofs());
+  q_soltrans.interpolate(q_old_solution, q_solution);
+  dgq_soltrans.interpolate(dgq_old_solution, dgq_solution);
+
+                               // check correctness by comparing the values
+                               // on points of QGauss of order 2.
+  MyFunction<dim> func;
+  {
+    double error = 0;
+    const hp::QCollection<dim> quad (QGauss<dim> (2));
+    hp::FEValues<dim> hp_fe_val (fe_q, quad, update_values |
+                                update_quadrature_points);
+    std::vector<double> vals (quad[0].size());
+    typename hp::DoFHandler<dim>::active_cell_iterator
+      cell = q_dof_handler.begin_active(),
+      endc = q_dof_handler.end();
+    for (; cell!=endc; ++cell)
+      {
+       hp_fe_val.reinit (cell, 0);
+       const FEValues<dim> &fe_val = hp_fe_val.get_present_fe_values();
+       fe_val.get_function_values (q_solution, vals);
+       for (unsigned int q=0; q<fe_val.n_quadrature_points; ++q)
+         {
+           error += std::fabs(func.value(fe_val.quadrature_point(q),0)-
+                              vals[q]);
+         }
+      }
+    deallog << "Error in interpolating hp FE_Q: " << error << std::endl;
+  }
+  {
+    double error = 0;
+    const hp::QCollection<dim> quad (QGauss<dim> (2));
+    hp::FEValues<dim> hp_fe_val (fe_dgq, quad, update_values |
+                                update_quadrature_points);
+    std::vector<double> vals (quad[0].size());
+    typename hp::DoFHandler<dim>::active_cell_iterator
+      celldg = dgq_dof_handler.begin_active(),
+      endc =  dgq_dof_handler.end();
+    for (; celldg!=endc; ++celldg)
+      {
+       hp_fe_val.reinit (celldg, 0);
+       const FEValues<dim> &fe_val = hp_fe_val.get_present_fe_values();
+       fe_val.get_function_values (dgq_solution, vals);
+       for (unsigned int q=0; q<fe_val.n_quadrature_points; ++q)
+         {
+           error += std::fabs(func.value(fe_val.quadrature_point(q),0)-
+                              vals[q]);
+         }
+      }
+    deallog << "Error in interpolating hp FE_DGQ: " << error << std::endl;
+  }
+}
+
+
+int main()
+{
+  std::ofstream logfile("solution_transfer/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  deallog << "   1D solution transfer" << std::endl;
+  transfer<1>(logfile);
+
+  deallog << "   2D solution transfer" << std::endl;
+  transfer<2>(logfile);
+
+  deallog << "   3D solution transfer" << std::endl;
+  transfer<3>(logfile);
+}
+
+
+
diff --git a/tests/hp/solution_transfer/cmp/generic b/tests/hp/solution_transfer/cmp/generic
new file mode 100644 (file)
index 0000000..7af228e
--- /dev/null
@@ -0,0 +1,10 @@
+
+DEAL::   1D solution transfer
+DEAL::Error in interpolating hp FE_Q: 0
+DEAL::Error in interpolating hp FE_DGQ: 0
+DEAL::   2D solution transfer
+DEAL::Error in interpolating hp FE_Q: 0
+DEAL::Error in interpolating hp FE_DGQ: 0
+DEAL::   3D solution transfer
+DEAL::Error in interpolating hp FE_Q: 0
+DEAL::Error in interpolating hp FE_DGQ: 0

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.