From 698081da3b321d4197dce32ed716847132330430 Mon Sep 17 00:00:00 2001
From: Wolfgang Bangerth <bangerth@math.tamu.edu>
Date: Wed, 3 May 2006 02:26:26 +0000
Subject: [PATCH] Add another test

git-svn-id: https://svn.dealii.org/trunk@13009 0785d39b-7218-0410-832d-ea1e28bc413d
---
 tests/hp/step-5.cc              | 339 ++++++++++++++++++++++++++++++++
 tests/hp/step-5/circle-grid.inp |  46 +++++
 tests/hp/step-5/cmp/generic     |  43 ++++
 3 files changed, 428 insertions(+)
 create mode 100644 tests/hp/step-5.cc
 create mode 100644 tests/hp/step-5/circle-grid.inp
 create mode 100644 tests/hp/step-5/cmp/generic

diff --git a/tests/hp/step-5.cc b/tests/hp/step-5.cc
new file mode 100644
index 0000000000..de9ac2b34f
--- /dev/null
+++ b/tests/hp/step-5.cc
@@ -0,0 +1,339 @@
+//----------------------------  step-5.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2005, 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.
+//
+//----------------------------  step-5.cc  ---------------------------
+
+
+// a hp-ified version of step-5
+
+
+#include <base/logstream.h>
+#include <fstream>
+std::ofstream logfile("step-5/output");
+
+
+#include <base/quadrature_lib.h>
+#include <base/function.h>
+#include <base/logstream.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 <grid/tria.h>
+#include <dofs/hp_dof_handler.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+#include <fe/fe_q.h>
+#include <fe/hp_fe_values.h>
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <numerics/data_out.h>
+
+#include <grid/grid_in.h>
+
+#include <grid/tria_boundary_lib.h>
+
+#include <fstream>
+#include <sstream>
+
+
+
+template <int dim>
+class LaplaceProblem 
+{
+  public:
+    LaplaceProblem ();
+    void run ();
+    
+  private:
+    void setup_system ();
+    void assemble_system ();
+    void solve ();
+    void output_results (const unsigned int cycle) const;
+
+    Triangulation<dim>   triangulation;
+    hp::FECollection<dim>            fe;
+    hp::DoFHandler<dim>      dof_handler;
+
+    SparsityPattern      sparsity_pattern;
+    SparseMatrix<double> system_matrix;
+
+    Vector<double>       solution;
+    Vector<double>       system_rhs;
+};
+
+
+
+template <int dim>
+class Coefficient : public Function<dim> 
+{
+  public:
+    Coefficient ()  : Function<dim>() {};
+    
+    virtual double value (const Point<dim>   &p,
+			  const unsigned int  component = 0) const;
+    
+    virtual void value_list (const std::vector<Point<dim> > &points,
+			     std::vector<double>            &values,
+			     const unsigned int              component = 0) const;
+};
+
+
+
+template <int dim>
+double Coefficient<dim>::value (const Point<dim> &p,
+				const unsigned int /*component*/) const 
+{
+  if (p.square() < 0.5*0.5)
+    return 20;
+  else
+    return 1;
+}
+
+
+
+
+template <int dim>
+void Coefficient<dim>::value_list (const std::vector<Point<dim> > &points,
+				   std::vector<double>            &values,
+				   const unsigned int              component) const 
+{
+  Assert (values.size() == points.size(), 
+	  ExcDimensionMismatch (values.size(), points.size()));
+  
+  Assert (component == 0, 
+	  ExcIndexRange (component, 0, 1));
+
+  const unsigned int n_points = points.size();
+
+  for (unsigned int i=0; i<n_points; ++i)
+    {
+      if (points[i].square() < 0.5*0.5)
+	values[i] = 20;
+      else
+	values[i] = 1;
+    }
+}
+
+
+
+
+template <int dim>
+LaplaceProblem<dim>::LaplaceProblem () :
+                fe (FE_Q<dim>(1)),
+		dof_handler (triangulation)
+{}
+
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::setup_system ()
+{
+  dof_handler.distribute_dofs (fe);
+
+  deallog << "   Number of degrees of freedom: "
+	    << dof_handler.n_dofs()
+	    << std::endl;
+
+  sparsity_pattern.reinit (dof_handler.n_dofs(),
+			   dof_handler.n_dofs(),
+			   dof_handler.max_couplings_between_dofs());
+  DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
+  sparsity_pattern.compress();
+
+  system_matrix.reinit (sparsity_pattern);
+
+  solution.reinit (dof_handler.n_dofs());
+  system_rhs.reinit (dof_handler.n_dofs());
+}
+
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::assemble_system () 
+{  
+  hp::QCollection<dim>  quadrature_formula(QGauss<dim>(2));
+
+  hp::FEValues<dim> x_fe_values (fe, quadrature_formula, 
+			   update_values    |  update_gradients |
+			   update_q_points  |  update_JxW_values);
+
+  const unsigned int   dofs_per_cell = fe[0].dofs_per_cell;
+  const unsigned int   n_q_points    = quadrature_formula[0].n_quadrature_points;
+
+  FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
+  Vector<double>       cell_rhs (dofs_per_cell);
+
+  std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+
+  const Coefficient<dim> coefficient;
+  std::vector<double>    coefficient_values (n_q_points);
+
+  typename hp::DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+  for (; cell!=endc; ++cell)
+    {
+      cell_matrix = 0;
+      cell_rhs = 0;
+
+      x_fe_values.reinit (cell);
+      const FEValues<2> &fe_values = x_fe_values.get_present_fe_values();
+
+      coefficient.value_list (fe_values.get_quadrature_points(),
+			      coefficient_values);
+      
+      for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+	for (unsigned int i=0; i<dofs_per_cell; ++i)
+	  {
+	    for (unsigned int j=0; j<dofs_per_cell; ++j)
+	      cell_matrix(i,j) += (coefficient_values[q_point] *
+				   fe_values.shape_grad(i,q_point) *
+				   fe_values.shape_grad(j,q_point) *
+				   fe_values.JxW(q_point));
+
+	    cell_rhs(i) += (fe_values.shape_value(i,q_point) *
+			    1.0 *
+			    fe_values.JxW(q_point));
+	  }
+
+
+      cell->get_dof_indices (local_dof_indices);
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+	{
+	  for (unsigned int j=0; j<dofs_per_cell; ++j)
+	    system_matrix.add (local_dof_indices[i],
+			       local_dof_indices[j],
+			       cell_matrix(i,j));
+	  
+	  system_rhs(local_dof_indices[i]) += cell_rhs(i);
+	}
+    }
+
+  std::map<unsigned int,double> boundary_values;
+  VectorTools::interpolate_boundary_values (dof_handler,
+					    0,
+					    ZeroFunction<dim>(),
+					    boundary_values);
+  MatrixTools::apply_boundary_values (boundary_values,
+				      system_matrix,
+				      solution,
+				      system_rhs);
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::solve () 
+{
+  SolverControl           solver_control (1000, 1e-12);
+  SolverCG<>              cg (solver_control);
+
+  PreconditionSSOR<> preconditioner;
+  preconditioner.initialize(system_matrix, 1.2);
+
+  cg.solve (system_matrix, solution, system_rhs,
+	    preconditioner);
+
+  deallog << "   " << solver_control.last_step()
+	    << " CG iterations needed to obtain convergence."
+	    << std::endl;
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::output_results (const unsigned int cycle) const
+{
+  DataOut<dim,hp::DoFHandler<dim> > data_out;
+
+  data_out.attach_dof_handler (dof_handler);
+  data_out.add_data_vector (solution, "solution");
+
+  data_out.build_patches ();
+
+  DataOutBase::EpsFlags eps_flags;
+  eps_flags.z_scaling = 4;
+  eps_flags.azimut_angle = 40;
+  eps_flags.turn_angle   = 10;
+  data_out.set_flags (eps_flags);
+
+  data_out.write_eps (deallog.get_file_stream());
+}
+
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::run () 
+{
+  for (unsigned int cycle=0; cycle<6; ++cycle)
+    {
+      deallog << "Cycle " << cycle << ':' << std::endl;
+
+      if (cycle != 0)
+	triangulation.refine_global (1);
+      else
+	{
+	  GridIn<dim> grid_in;
+	  grid_in.attach_triangulation (triangulation);
+	  std::ifstream input_file("step-5/circle-grid.inp");
+	  Assert (dim==2, ExcInternalError());
+	  
+	  grid_in.read_ucd (input_file);
+
+	  static const HyperBallBoundary<dim> boundary;
+	  triangulation.set_boundary (0, boundary);
+	}
+
+      deallog << "   Number of active cells: "
+		<< triangulation.n_active_cells()
+		<< std::endl
+		<< "   Total number of cells: "
+		<< triangulation.n_cells()
+		<< std::endl;
+
+      setup_system ();
+      assemble_system ();
+      solve ();
+      output_results (cycle);
+    }
+}
+
+
+
+int main () 
+{
+  logfile.precision(2);
+  
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);  
+
+  deallog.depth_console (0);
+
+  LaplaceProblem<2> laplace_problem_2d;
+  laplace_problem_2d.run ();
+
+/*  
+  Coefficient<2>    coefficient;
+  std::vector<Point<2> > points (2);
+  std::vector<double>    coefficient_values (1);
+  coefficient.value_list (points, coefficient_values);
+*/
+  
+  return 0;
+}
diff --git a/tests/hp/step-5/circle-grid.inp b/tests/hp/step-5/circle-grid.inp
new file mode 100644
index 0000000000..f28a7a238a
--- /dev/null
+++ b/tests/hp/step-5/circle-grid.inp
@@ -0,0 +1,46 @@
+25 20 0 0 0
+1  -0.7071 -0.7071 0
+2  0.7071 -0.7071 0
+3  -0.2668 -0.2668 0
+4  0.2668 -0.2668 0
+5  -0.2668 0.2668 0
+6  0.2668 0.2668 0
+7  -0.7071 0.7071 0
+8  0.7071 0.7071 0
+9  0 -1 0
+10  0.5 -0.5 0
+11  0 -0.3139 0
+12  -0.5 -0.5 0
+13  0 -0.6621 0
+14  -0.3139 0 0
+15  -0.5 0.5 0
+16  -1 0 0
+17  -0.6621 0 0
+18  0.3139 0 0
+19  0 0.3139 0
+20  0 0 0
+21  1 0 0
+22  0.5 0.5 0
+23  0.6621 0 0
+24  0 1 0
+25  0 0.6621 0
+1 0 quad    1 9 13 12 
+2 0 quad    9 2 10 13 
+3 0 quad    13 10 4 11 
+4 0 quad    12 13 11 3 
+5 0 quad    1 12 17 16 
+6 0 quad    12 3 14 17 
+7 0 quad    17 14 5 15 
+8 0 quad    16 17 15 7 
+9 0 quad    3 11 20 14 
+10 0 quad    11 4 18 20 
+11 0 quad    20 18 6 19 
+12 0 quad    14 20 19 5 
+13 0 quad    2 21 23 10 
+14 0 quad    21 8 22 23 
+15 0 quad    23 22 6 18 
+16 0 quad    10 23 18 4 
+17 0 quad    7 15 25 24 
+18 0 quad    15 5 19 25 
+19 0 quad    25 19 6 22 
+20 0 quad    24 25 22 8 
diff --git a/tests/hp/step-5/cmp/generic b/tests/hp/step-5/cmp/generic
new file mode 100644
index 0000000000..309c3cb828
--- /dev/null
+++ b/tests/hp/step-5/cmp/generic
@@ -0,0 +1,43 @@
+
+DEAL::Cycle 0:
+DEAL::   Number of active cells: 20
+DEAL::   Total number of cells: 20
+DEAL::   Number of degrees of freedom: 25
+DEAL:cg::Starting value 0.51
+DEAL:cg::Convergence step 13 value 0
+DEAL::   13 CG iterations needed to obtain convergence.
+DEAL::Cycle 1:
+DEAL::   Number of active cells: 80
+DEAL::   Total number of cells: 100
+DEAL::   Number of degrees of freedom: 89
+DEAL:cg::Starting value 0.31
+DEAL:cg::Convergence step 18 value 0
+DEAL::   18 CG iterations needed to obtain convergence.
+DEAL::Cycle 2:
+DEAL::   Number of active cells: 320
+DEAL::   Total number of cells: 420
+DEAL::   Number of degrees of freedom: 337
+DEAL:cg::Starting value 0.17
+DEAL:cg::Convergence step 29 value 0
+DEAL::   29 CG iterations needed to obtain convergence.
+DEAL::Cycle 3:
+DEAL::   Number of active cells: 1280
+DEAL::   Total number of cells: 1700
+DEAL::   Number of degrees of freedom: 1313
+DEAL:cg::Starting value 0.092
+DEAL:cg::Convergence step 52 value 0
+DEAL::   52 CG iterations needed to obtain convergence.
+DEAL::Cycle 4:
+DEAL::   Number of active cells: 5120
+DEAL::   Total number of cells: 6820
+DEAL::   Number of degrees of freedom: 5185
+DEAL:cg::Starting value 0.047
+DEAL:cg::Convergence step 95 value 0
+DEAL::   95 CG iterations needed to obtain convergence.
+DEAL::Cycle 5:
+DEAL::   Number of active cells: 20480
+DEAL::   Total number of cells: 27300
+DEAL::   Number of degrees of freedom: 20609
+DEAL:cg::Starting value 0.024
+DEAL:cg::Convergence step 182 value 0
+DEAL::   182 CG iterations needed to obtain convergence.
-- 
2.39.5