From fcb42ced3f4219e123b752b67ebb9db681f9e17b Mon Sep 17 00:00:00 2001
From: bangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Thu, 9 Jan 2014 00:18:50 +0000
Subject: [PATCH] Fix the same bug in one more place.

git-svn-id: https://svn.dealii.org/trunk@32182 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/source/fe/fe_values.cc            |  22 ++--
 tests/hp/do_function_laplacians_01.cc     | 120 ++++++++++++++++++++++
 tests/hp/do_function_laplacians_01.output |   2 +
 3 files changed, 135 insertions(+), 9 deletions(-)
 create mode 100644 tests/hp/do_function_laplacians_01.cc
 create mode 100644 tests/hp/do_function_laplacians_01.output

diff --git a/deal.II/source/fe/fe_values.cc b/deal.II/source/fe/fe_values.cc
index c7a95ede57..aaf9fe0a0a 100644
--- a/deal.II/source/fe/fe_values.cc
+++ b/deal.II/source/fe/fe_values.cc
@@ -2385,8 +2385,7 @@ namespace internal
       return;
 
 
-    const unsigned int n_quadrature_points = dofs_per_cell > 0 ?
-                                             shape_derivatives[0].size() : 0;
+    const unsigned int n_quadrature_points = shape_derivatives[0].size();
     const unsigned int n_components = fe.n_components();
 
     // Assert that we can write all components into the result vectors
@@ -2495,9 +2494,19 @@ namespace internal
                           const bool quadrature_points_fastest  = false,
                           const unsigned int component_multiple = 1)
   {
+    // initialize with zero
+    for (unsigned int i=0; i<laplacians.size(); ++i)
+      std::fill_n (laplacians[i].begin(), laplacians[i].size(),
+                   typename VectorType::value_type());
+
+    // see if there the current cell has DoFs at all, and if not
+    // then there is nothing else to do.
     const unsigned int dofs_per_cell = fe.dofs_per_cell;
-    const unsigned int n_quadrature_points = dofs_per_cell > 0 ?
-                                             shape_hessians[0].size() : 0;
+    if (dofs_per_cell == 0)
+      return;
+
+
+    const unsigned int n_quadrature_points = shape_hessians[0].size();
     const unsigned int n_components = fe.n_components();
 
     // Assert that we can write all components into the result vectors
@@ -2515,11 +2524,6 @@ namespace internal
           AssertDimension (laplacians[i].size(), result_components);
       }
 
-    // initialize with zero
-    for (unsigned int i=0; i<laplacians.size(); ++i)
-      std::fill_n (laplacians[i].begin(), laplacians[i].size(),
-                   typename VectorType::value_type());
-
     // add up contributions of trial functions.  now check whether the shape
     // function is primitive or not. if it is, then set its only non-zero
     // component, otherwise loop over components
diff --git a/tests/hp/do_function_laplacians_01.cc b/tests/hp/do_function_laplacians_01.cc
new file mode 100644
index 0000000000..3cd68e0644
--- /dev/null
+++ b/tests/hp/do_function_laplacians_01.cc
@@ -0,0 +1,120 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2011 - 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// testcase by Minh Do-Quang: a case where SolutionTransfer got into trouble
+// in a couple of places when using FE_Nothing and FESystem.
+
+#include "../tests.h"
+#include <fstream>
+#include <sstream>
+#include <iostream>
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/hp/fe_collection.h>
+
+#include <deal.II/numerics/solution_transfer.h>
+
+#include <deal.II/numerics/data_out.h>
+
+using namespace dealii;
+
+int main()
+{
+  std::ofstream logfile("output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  Triangulation<2> triangulation;
+  GridGenerator::hyper_cube (triangulation);
+  triangulation.refine_global (1);
+
+
+  hp::FECollection<2>     fe_collection;
+
+  fe_collection.push_back(FESystem<2>(FE_Q<2>(1), 1, FE_Q<2>(1), 1));
+  fe_collection.push_back(FESystem<2>(FE_Nothing<2>(), 1, FE_Nothing<2>(), 1));
+
+  hp::DoFHandler<2> dof_handler(triangulation);
+
+  // Assign FE to cells
+  hp::DoFHandler<2>::active_cell_iterator cell;
+  hp::DoFHandler<2>::active_cell_iterator endc = dof_handler.end();      
+
+
+  cell = dof_handler.begin_active();
+  cell->set_active_fe_index(1);
+  cell++;
+  cell->set_active_fe_index(1);
+  cell++;
+  cell->set_active_fe_index(0);
+  cell++;
+  cell->set_active_fe_index(0);
+
+
+  dof_handler.distribute_dofs (fe_collection);
+
+  // Init solution
+  Vector<double> solution(dof_handler.n_dofs());
+  solution = 1.0;
+
+
+  SolutionTransfer<2, Vector<double>, hp::DoFHandler<2> > solultion_trans(dof_handler);
+  solultion_trans.prepare_for_coarsening_and_refinement(solution);
+
+  triangulation.execute_coarsening_and_refinement ();
+  dof_handler.distribute_dofs (fe_collection);
+  
+  Vector<double> new_solution(dof_handler.n_dofs());
+  solultion_trans.interpolate(solution, new_solution);
+
+  hp::QCollection<2> q;
+  q.push_back(QMidpoint<2>());
+  q.push_back(QMidpoint<2>());
+  hp::FEValues<2> x_fe_values (fe_collection, q, update_hessians);
+  for (typename hp::DoFHandler<2>::active_cell_iterator  cell = dof_handler.begin_active();
+       cell != dof_handler.end(); ++cell)
+    {
+      x_fe_values.reinit (cell);
+      std::vector<Vector<double> >
+	derivatives (q[0].size(), Vector<double>(cell->get_fe().n_components()));
+      
+      x_fe_values.get_present_fe_values().get_function_laplacians (new_solution,
+								   derivatives);
+    }
+  
+  // we are good if we made it to here
+  deallog << "OK" << std::endl;
+}
+
+
+
diff --git a/tests/hp/do_function_laplacians_01.output b/tests/hp/do_function_laplacians_01.output
new file mode 100644
index 0000000000..0fd8fc12f0
--- /dev/null
+++ b/tests/hp/do_function_laplacians_01.output
@@ -0,0 +1,2 @@
+
+DEAL::OK
-- 
2.39.5