Added a testcase for the hanging_node_constraints in the hp-case and

modified test step-3c to work properly in the hp-setting.

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

diff --git a/tests/hp/Makefile b/tests/hp/Makefile
index f96dfe52f41ebc3913c72aa28400bafd66f34591..908f7e67097d3890e45851e94a50a9ee714abdeb 100644 (file)
--- a/tests/hp/Makefile
+++ b/tests/hp/Makefile
@@ -32,7 +32,8 @@ tests_x = hp_dof_handler \
          step-* \
          n_active_fe_indices \
          integrate_difference \
-         hp_*_dof_identities_*
+         hp_*_dof_identities_* \
+         hp_hanging_nodes_*
 
 # from above list of regular expressions, generate the real set of
 # tests

diff --git a/tests/hp/hp_hanging_nodes_01.cc b/tests/hp/hp_hanging_nodes_01.cc
new file mode 100644 (file)
index 0000000..3698604
--- /dev/null
+++ b/tests/hp/hp_hanging_nodes_01.cc
@@ -0,0 +1,178 @@
+/*    $Id$       */
+/*    Version: $Name$                                          */
+/*                                                                */
+/*    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2006 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.                        */
+
+
+#include <grid/tria.h>
+#include <dofs/dof_handler.h>
+#include <dofs/hp_dof_handler.h>
+#include <grid/grid_generator.h>
+#include <grid/grid_out.h>
+
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <dofs/dof_accessor.h>
+
+#include <fe/fe_q.h>
+
+#include <dofs/dof_tools.h>
+#include <dofs/dof_constraints.h>
+
+#include <fe/fe_values.h>
+#include <fe/fe_collection.h>
+#include <base/quadrature_lib.h>
+
+#include <base/function.h>
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+
+#include <lac/vector.h>
+#include <lac/full_matrix.h>
+#include <lac/sparse_matrix.h>
+#include <lac/solver_cg.h>
+#include <lac/precondition.h>
+
+#include <numerics/data_out.h>
+#include <fstream>
+#include <iostream>
+
+
+// This a test for the hp capable version of the make_hanging_node_constraints
+// method. It uses a triangulation with one refined element beside an
+// unrefined element to create the constraints for this configuration.
+
+template <int dim>
+int generate_grid (Triangulation<dim> &tria)
+{
+  Point<dim> p1,
+    p2;
+  std::vector<unsigned int> sub_div;
+
+                                  // Define a rectangular shape
+  for (unsigned int d=0; d < dim; ++d)
+    {
+      p1(d) = 0;
+      p2(d) = (d == 0) ? 2.0 : 1.0;
+      sub_div.push_back ( (d == 0) ? 2 : 1);
+    }  
+  GridGenerator::subdivided_hyper_rectangle (tria, sub_div, p1, p2, true);
+
+                                  // Refine the first cell.
+  tria.begin_active ()->set_refine_flag ();
+  tria.execute_coarsening_and_refinement ();
+  
+  return (0);
+}
+
+
+
+template <int dim>
+void test_constraints (hp::FECollection<dim> &fe_coll)
+{
+  Triangulation<dim> tria;
+
+                                  // Setup a rectangular domain
+                                  // where one cell is h-refined,
+                                  // while the other cell is
+                                  // unrefined. Furthermore every cell
+                                  // gets a different active_fe_index.
+                                  // This should serve as a testcase
+                                  // for the hanging node constraints.
+  generate_grid (tria);
+
+                                  // Now assign increasing
+                                  // active_fe_indices to
+                                  // the different cells.
+  hp::DoFHandler<dim> dof_handler (tria);
+  typename hp::DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active (),
+                                           endc = dof_handler.end ();
+  unsigned int fe_indx = 0;
+  for (; cell != endc; ++cell)
+    {
+      cell->set_active_fe_index (fe_indx);
+      ++fe_indx;
+    }
+  
+                                  // Distribute DoFs;
+  dof_handler.distribute_dofs (fe_coll);
+  deallog << "DoFs: " << dof_handler.n_dofs () << std::endl;
+  
+                                  // Create the constraints.
+  ConstraintMatrix constraint_matrix;
+  
+  DoFTools::make_hanging_node_constraints (dof_handler,
+                                          constraint_matrix);
+
+                                  // Output the constraints
+  constraint_matrix.print (deallog.get_file_stream ());  
+}
+
+
+template <int dim>
+void test_constraints_old (FiniteElement<dim> &fe)
+{
+  Triangulation<dim> tria;
+
+                                  // Setup a rectangular domain
+                                  // where one cell is h-refined,
+                                  // while the other cell is
+                                  // unrefined. Furthermore every cell
+                                  // gets a different active_fe_index.
+                                  // This should serve as a testcase
+                                  // for the hanging node constraints.
+  generate_grid (tria);
+
+                                  // Now assign increasing
+                                  // active_fe_indices to
+                                  // the different cells.
+  DoFHandler<dim> dof_handler (tria);
+  typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active (),
+                                                endc = dof_handler.end ();
+  
+                                  // Distribute DoFs;
+  dof_handler.distribute_dofs (fe);
+  deallog << "DoFs: " << dof_handler.n_dofs () << std::endl;
+  
+                                  // Create the constraints.
+  ConstraintMatrix constraint_matrix;
+  
+  DoFTools::make_hanging_node_constraints (dof_handler,
+                                          constraint_matrix);
+
+                                  // Output the constraints
+  constraint_matrix.print (deallog.get_file_stream ());
+}
+
+int main () 
+{
+  std::ofstream logfile("hp_hanging_nodes_01/output");
+  logfile.precision(2);
+  
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+  
+  FE_Q<2> fe_1 (1);
+  FE_Q<2> fe_2 (2);
+  FE_Q<2> fe_3 (3);
+
+  hp::FECollection<2> fe_coll2;
+  fe_coll2.push_back (fe_3);
+  fe_coll2.push_back (fe_2);
+  fe_coll2.push_back (fe_2);
+
+  fe_coll2.push_back (fe_2);
+  fe_coll2.push_back (fe_3); 
+  
+  test_constraints<2> (fe_coll2);
+
+  test_constraints_old<2> (fe_1);
+
+  return 0;
+}

diff --git a/tests/hp/step-3c.cc b/tests/hp/step-3c.cc
index 96cc6bb5836b1149f8c06fb40a71520ff58abb10..b2737bffc2121fe764114749296b2993c05d72ef 100644 (file)
--- a/tests/hp/step-3c.cc
+++ b/tests/hp/step-3c.cc
@@ -144,12 +144,15 @@ void LaplaceProblem::make_grid_and_dofs ()
 
 void LaplaceProblem::assemble_system () 
 {
-  hp::QCollection<2>  quadrature_formula(QGauss<2>(6));
+  hp::QCollection<2>  quadrature_formula;
+
+  for (unsigned int p = 0; p < fe.size (); ++p)
+    quadrature_formula.push_back (QGauss<2> (p + 2));
+
   hp::FEValues<2> x_fe_values (fe, quadrature_formula, 
                               update_values | update_gradients | update_JxW_values);
   
   const unsigned int   max_dofs_per_cell = fe.max_dofs_per_cell ();
-  const unsigned int   n_q_points    = quadrature_formula[0].n_quadrature_points;
 
   FullMatrix<double>   cell_matrix (max_dofs_per_cell, max_dofs_per_cell);
   Vector<double>       cell_rhs (max_dofs_per_cell);
@@ -169,7 +172,7 @@ void LaplaceProblem::assemble_system ()
       cell_rhs = 0;
 
       const unsigned int dofs_per_cell = cell->get_fe ().dofs_per_cell;
-      
+      const unsigned int   n_q_points    = quadrature_formula[cell->active_fe_index ()].n_quadrature_points;
       for (unsigned int i=0; i<dofs_per_cell; ++i)
        for (unsigned int j=0; j<dofs_per_cell; ++j)
          for (unsigned int q_point=0; q_point<n_q_points; ++q_point)