Fix a testcase by Minh Do-Quang: DataOut also got with FE_Nothing elements.

git-svn-id: https://svn.dealii.org/trunk@32179 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h
index 444e2a9d142539372f0581d77bced8415327ad56..2b1f3d07b5e98f27f0d03f6cbb2c7340db2451ee 100644 (file)
--- a/deal.II/doc/news/changes.h
+++ b/deal.II/doc/news/changes.h
@@ -85,9 +85,16 @@ inconvenience this causes.
 
 <ol>
 
+  <li> Fixed: DataOut got confused in some situations where one uses FE_Nothing.
+  This is now fixed.
+  <br>
+  (Minh Do-Quang, Wolfgang Bangerth, 2014/01/08)
+  </li>
+
   <li> Fixed: FESystem::get_interpolation_matrix, a function that is among
   other places used by SolutionTransfer, had a bug that prevented it from
   running correctly in some situations where one uses FE_Nothing.
+  This is now fixed.
   <br>
   (Minh Do-Quang, Wolfgang Bangerth, 2014/01/08)
   </li>

diff --git a/deal.II/source/fe/fe_values.cc b/deal.II/source/fe/fe_values.cc
index 066c79c4c870fa4b8439a46357df6d2d1e2baf98..7e00704767472b3541ffc32ffe4f5550d131ca65 100644 (file)
--- a/deal.II/source/fe/fe_values.cc
+++ b/deal.II/source/fe/fe_values.cc
@@ -1,7 +1,7 @@
 // ---------------------------------------------------------------------
 // $Id$
 //
-// Copyright (C) 1998 - 2013 by the deal.II authors
+// Copyright (C) 1998 - 2014 by the deal.II authors
 //
 // This file is part of the deal.II library.
 //
@@ -2242,9 +2242,18 @@ namespace internal
                       const bool quadrature_points_fastest  = false,
                       const unsigned int component_multiple = 1)
   {
+    // initialize with zero
+    for (unsigned int i=0; i<values.size(); ++i)
+      std::fill_n (values[i].begin(), values[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_values.n_cols() : 0;
+    if (dofs_per_cell == 0)
+      return;
+    
+    const unsigned int n_quadrature_points = shape_values.n_cols();
     const unsigned int n_components = fe.n_components();
 
     // Assert that we can write all components into the result vectors
@@ -2262,11 +2271,6 @@ namespace internal
           AssertDimension (values[i].size(), result_components);
       }
 
-    // initialize with zero
-    for (unsigned int i=0; i<values.size(); ++i)
-      std::fill_n (values[i].begin(), values[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/solution_transfer_05.cc b/tests/hp/solution_transfer_05.cc
new file mode 100644 (file)
index 0000000..a7d47b2
--- /dev/null
+++ b/tests/hp/solution_transfer_05.cc
@@ -0,0 +1,112 @@
+// ---------------------------------------------------------------------
+// $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/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);
+
+  // a follow-up error to the one fixed with _04 was that DataOut also got
+  // itself confused
+  DataOut<2, hp::DoFHandler<2> > data_out2;
+  data_out2.attach_dof_handler (dof_handler);
+  data_out2.add_data_vector (new_solution, "Solution");
+  data_out2.build_patches();
+  data_out2.write_vtu(deallog.get_file_stream());
+  
+  // we are good if we made it to here
+  deallog << "OK" << std::endl;
+}
+
+
+

diff --git a/tests/hp/solution_transfer_05.output b/tests/hp/solution_transfer_05.output
new file mode 100644 (file)
index 0000000..437080e
--- /dev/null
+++ b/tests/hp/solution_transfer_05.output
@@ -0,0 +1,36 @@
+
+<?xml version="1.0" ?> 
+<!-- 
+# vtk DataFile Version 3.0
+#This file was generated 
+-->
+<VTKFile type="UnstructuredGrid" version="0.1" compressor="vtkZLibDataCompressor" byte_order="LittleEndian">
+<UnstructuredGrid>
+<Piece NumberOfPoints="16" NumberOfCells="4" >
+  <Points>
+    <DataArray type="Float64" NumberOfComponents="3" format="binary">
+AQAAAIABAACAAQAAMwAAAA==eNpjYMAHHtgzkCQP4+MSRwcf7PGbi67vgz1+84h1Dy77Yeo+2FPmHlzmfMAQBwBtNxtp
+    </DataArray>
+  </Points>
+
+  <Cells>
+    <DataArray type="Int32" Name="connectivity" format="binary">
+AQAAAEAAAABAAAAAKAAAAA==eNoNwwkSwBAQALB1VIvi/7+VzCQiIlnMVh9fm5/d6fB3edxeCvAAeQ==
+    </DataArray>
+    <DataArray type="Int32" Name="offsets" format="binary">
+AQAAABAAAAAQAAAAEwAAAA==eNpjYWBg4ABiHiAWAGIAAVAAKQ==
+    </DataArray>
+    <DataArray type="UInt8" Name="types" format="binary">
+AQAAAAQAAAAEAAAADAAAAA==eNrj5OTkBAAAXgAl
+    </DataArray>
+  </Cells>
+  <PointData Scalars="scalars">
+    <DataArray type="Float64" Name="Solution_0" format="binary">
+AQAAAIAAAACAAAAAEAAAAA==eNpjYKAG+GBPLg0AGqcJeQ==    </DataArray>
+    <DataArray type="Float64" Name="Solution_1" format="binary">
+AQAAAIAAAACAAAAAEAAAAA==eNpjYKAG+GBPLg0AGqcJeQ==    </DataArray>
+  </PointData>
+ </Piece>
+ </UnstructuredGrid>
+</VTKFile>
+DEAL::OK