From: Denis Davydov <davydden@gmail.com>
Date: Mon, 13 Apr 2015 09:19:11 +0000 (+0200)
Subject: added a SolutionTransfer test for FE_Q_Hierarchical.
X-Git-Tag: v8.3.0-rc1~284^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=8f4f6198a0c421ec6559d8dbf255b6d6fb722e81;p=dealii.git

added a SolutionTransfer test for FE_Q_Hierarchical.

Also cleaned up unnecessary Assert in the class.
---

diff --git a/source/fe/fe_q_hierarchical.cc b/source/fe/fe_q_hierarchical.cc
index d4a4017ee2..90e6d1ed71 100644
--- a/source/fe/fe_q_hierarchical.cc
+++ b/source/fe/fe_q_hierarchical.cc
@@ -164,10 +164,6 @@ FE_Q_Hierarchical<dim>::get_interpolation_matrix(const FiniteElement< dim> &sour
               ExcDimensionMismatch (matrix.m(),
                                     source_fe->dofs_per_cell));
 
-      // make sure that we interpolate to a richer element.
-      Assert (this->dofs_per_cell >= source_fe->dofs_per_cell,
-              ExcMessage("Interpolation to higher polynomial degree is advised."));
-
       // Recall that DoFs are renumbered in the following order:
       // vertices, lines, quads, hexes.
       // As we deal with hierarchical FE, interpolation matrix is rather easy:
diff --git a/tests/hp/solution_transfer_12.cc b/tests/hp/solution_transfer_12.cc
new file mode 100644
index 0000000000..c2e29e3fda
--- /dev/null
+++ b/tests/hp/solution_transfer_12.cc
@@ -0,0 +1,238 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+#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/vector_tools.h>
+#include <deal.II/numerics/solution_transfer.h>
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/fe/fe_q_hierarchical.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
+// This is a copy-paste of solution_transfer.cc test,
+// with the two exceptions:
+// - FE_Q_Hierarchical
+// - manual interpolation as there is no general interpolation implemented
+//   for this FE currently, thus VectorTools::interpolate can't be used.
+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> func;
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria);
+  tria.refine_global(5-dim);
+  const unsigned int max_degree = 6-dim;
+  hp::FECollection<dim> fe_q;
+  for (unsigned int deg=1; deg<=max_degree; ++deg)
+    {
+      fe_q.push_back (FE_Q_Hierarchical<dim>(deg));
+    }
+  hp::DoFHandler<dim> q_dof_handler(tria);
+  Vector<double> q_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(),
+    endc=q_dof_handler.end();
+    for (; cell!=endc; ++cell, ++counter)
+      {
+        if (counter < 15)
+          cell->set_active_fe_index(1);
+        else
+          cell->set_active_fe_index(Testing::rand()%max_degree);
+      }
+  }
+
+  q_dof_handler.distribute_dofs (fe_q);
+  q_solution.reinit(q_dof_handler.n_dofs());
+
+  ConstraintMatrix cm;
+  cm.close();
+
+  // interpolation is not straight forward for FE_Q_Hierarchical.
+  // Since we know that we do linear function, just get those DoFs
+  // which correspond to vertices and assign a value of shape function there
+  // while keeping other modes zero.
+  {
+    q_solution = 0.;
+
+    hp::QCollection<dim> quad;
+    quad.push_back (QTrapez<dim> ());
+    hp::FEValues<dim> hp_fe_val (fe_q, quad, update_values |
+                                 update_quadrature_points);
+    std::vector<unsigned int> local_dof_indices;
+    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();
+        local_dof_indices.resize(fe_val.dofs_per_cell);
+
+        cell->get_dof_indices(local_dof_indices);
+
+        Assert (local_dof_indices.size() >= fe_val.n_quadrature_points,
+                ExcInternalError());
+
+        for (unsigned int q=0; q<fe_val.n_quadrature_points; ++q)
+          {
+            q_solution[local_dof_indices[q]] = func.value(fe_val.quadrature_point(q),0);
+          }
+      }
+
+  }
+
+  SolutionTransfer<dim,Vector<double>,hp::DoFHandler<dim> > q_soltrans(q_dof_handler);
+
+
+  // test b): do some coarsening and
+  // refinement
+  q_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 (Testing::rand() % 3 == 0)
+        cell->set_refine_flag();
+      else if (Testing::rand() % 3 == 3)
+        cell->set_coarsen_flag();
+    }
+
+  Vector<double> q_old_solution=q_solution;
+  tria.prepare_coarsening_and_refinement();
+  q_soltrans.prepare_for_coarsening_and_refinement(q_old_solution);
+  tria.execute_coarsening_and_refinement();
+
+  counter = 0;
+  {
+    typename hp::DoFHandler<dim>::active_cell_iterator
+    cell = q_dof_handler.begin_active(),
+    endc = q_dof_handler.end();
+    for (; cell!=endc; ++cell, ++counter)
+      {
+        if (counter > 20 && counter < 90)
+          cell->set_active_fe_index(0);
+        else
+          cell->set_active_fe_index(Testing::rand()%max_degree);
+      }
+  }
+
+  q_dof_handler.distribute_dofs (fe_q);
+  q_solution.reinit(q_dof_handler.n_dofs());
+  q_soltrans.interpolate(q_old_solution, q_solution);
+
+  // check correctness by comparing the values
+  // on points of QGauss of order 2.
+  {
+    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_Hierarchical: " << error << std::endl;
+  }
+}
+
+
+int main()
+{
+  std::ofstream logfile("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_12.output b/tests/hp/solution_transfer_12.output
new file mode 100644
index 0000000000..adc30e4418
--- /dev/null
+++ b/tests/hp/solution_transfer_12.output
@@ -0,0 +1,7 @@
+
+DEAL::   1D solution transfer
+DEAL::Error in interpolating hp FE_Q_Hierarchical: 0
+DEAL::   2D solution transfer
+DEAL::Error in interpolating hp FE_Q_Hierarchical: 0
+DEAL::   3D solution transfer
+DEAL::Error in interpolating hp FE_Q_Hierarchical: 0