From 34fb3288ec1c6a91e17bf3b87cde2e06d718878e Mon Sep 17 00:00:00 2001
From: wolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Sun, 2 Jan 2000 20:34:42 +0000
Subject: [PATCH] Add step-6.

git-svn-id: https://svn.dealii.org/trunk@2141 0785d39b-7218-0410-832d-ea1e28bc413d
---
 .../examples/step-by-step/step-6/Makefile     |  94 +++++
 .../examples/step-by-step/step-6/step-6.cc    | 349 ++++++++++++++++++
 deal.II/examples/step-6/Makefile              |  94 +++++
 deal.II/examples/step-6/step-6.cc             | 349 ++++++++++++++++++
 4 files changed, 886 insertions(+)
 create mode 100644 deal.II/deal.II/Attic/examples/step-by-step/step-6/Makefile
 create mode 100644 deal.II/deal.II/Attic/examples/step-by-step/step-6/step-6.cc
 create mode 100644 deal.II/examples/step-6/Makefile
 create mode 100644 deal.II/examples/step-6/step-6.cc

diff --git a/deal.II/deal.II/Attic/examples/step-by-step/step-6/Makefile b/deal.II/deal.II/Attic/examples/step-by-step/step-6/Makefile
new file mode 100644
index 0000000000..dc2798f839
--- /dev/null
+++ b/deal.II/deal.II/Attic/examples/step-by-step/step-6/Makefile
@@ -0,0 +1,94 @@
+# $Id$
+# Copyright W. Bangerth, University of Heidelberg, 1999
+
+# Template for makefiles for the examples subdirectory. In principle,
+# everything should be done automatically if you set the target file
+# here correctly. We get deduce it from the files in the present
+# directory:
+target   = $(basename $(shell echo step-*.cc))
+
+# All dependencies between files should be updated by the included
+# file Makefile.dep if necessary. Object files are compiled into
+# the archives ./Obj.a and ./Obj.g.a. By default, the debug version
+# is used to link. It you don't like that, change the following
+# variable to "off"
+debug-mode = on
+
+
+
+###############################################################################
+# Internals
+
+#deal include base path
+D = ../../../..
+
+include $D/common/Make.global_options
+
+# get lists of files we need
+
+
+# list of libraries needed to link with
+libs     = -ldeal_II_2d  -llac -lbase
+libs.g   = -ldeal_II_2d.g -llac.g -lbase.g
+
+
+# check whether we use debug mode or not
+ifeq ($(debug-mode),on)
+  libraries = $(target).go $(libs.g)
+  flags     = $(CXXFLAGS.g)
+else
+  libraries = $(target).go $(libs)
+  flags     = $(CXXFLAGS.o)
+endif
+
+
+
+# make rule for the target. $^ is the object file $(target).g?o
+$(target) : $(libraries)
+	@echo ============================ Linking $@
+	@$(CXX) $(flags) -o $@ $^
+
+# rule how to run the program
+run: $(target)
+	@echo ============================ Running $@
+	@./$(target)
+
+
+# rule to make object files
+%.go : %.cc
+	@echo ============================ Compiling with debugging information:   $<
+	@$(CXX) $(CXXFLAGS.g) -c $< -o $@
+%.o : %.cc
+	@echo ============================ Compiling with optimization:   $<
+	@$(CXX) $(CXXFLAGS) -c $< -o $@
+
+
+clean:
+	-rm -f *.o *.go *~ Makefile.dep $(target) *gmv *gnuplot *gpl *eps
+
+
+
+.PHONY: clean
+
+
+# Rule to generate the dependency file. This file is
+# automagically remade whenever needed, i.e. whenever
+# one of the cc-/h-files changed. Make detects whether
+# to remake this file upon inclusion at the bottom
+# of this file.
+#
+# Since the script prefixes the output names by lib/g?o, we have to
+# strip that again (the script was written for the main libraries and
+# large projects where object files are put into subdirs)
+Makefile.dep: $(target).cc Makefile \
+              $(shell echo $D/base/include/base/*.h \
+                           $D/lac/include/lac/*.h  \
+                           $D/deal.II/include/*/*.h)
+	@echo ============================ Remaking Makefile
+	@perl $D/common/scripts/make_dependencies.pl  $(INCLUDE) $(target).cc \
+		| perl -pi -e 's!lib/g?o/!!g;' \
+		> Makefile.dep
+
+
+include Makefile.dep
+
diff --git a/deal.II/deal.II/Attic/examples/step-by-step/step-6/step-6.cc b/deal.II/deal.II/Attic/examples/step-by-step/step-6/step-6.cc
new file mode 100644
index 0000000000..b87329d1dc
--- /dev/null
+++ b/deal.II/deal.II/Attic/examples/step-by-step/step-6/step-6.cc
@@ -0,0 +1,349 @@
+/* $Id$ */
+/* Author: Wolfgang Bangerth, University of Heidelberg, 2000 */
+
+
+				 // still unfinished
+#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/vector_memory.h>
+#include <lac/precondition.h>
+#include <grid/tria.h>
+#include <dofs/dof_handler.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+#include <fe/fe_lib.lagrange.h>
+#include <fe/fe_values.h>
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <numerics/data_out.h>
+
+				 // out statt in
+#include <grid/grid_out.h>
+
+#include <grid/tria_boundary_lib.h>
+
+				 //...
+#include <dofs/dof_constraints.h>
+#include <numerics/error_estimator.h>
+
+#include <fstream>
+
+
+
+template <int dim>
+class LaplaceProblem 
+{
+  public:
+    LaplaceProblem ();
+    void run ();
+    
+  private:
+    void setup_system ();
+    void assemble_system ();
+    void solve ();
+    void refine_grid ();
+    void output_results (const unsigned int cycle) const;
+
+    Triangulation<dim>   triangulation;
+    FEQ1<dim>            fe;
+    DoFHandler<dim>      dof_handler;
+
+    ConstraintMatrix     hanging_node_constraints;
+
+    SparseMatrixStruct   sparsity_pattern;
+    SparseMatrix<double> system_matrix;
+
+    Vector<double>       solution;
+    Vector<double>       system_rhs;
+};
+
+
+
+template <int dim>
+class Coefficient : public Function<dim> 
+{
+  public:
+    virtual double value (const Point<dim>   &p,
+			  const unsigned int  component = 0) const;
+    
+    virtual void value_list (const vector<Point<dim> > &points,
+			     vector<double>            &values,
+			     const unsigned int         component = 0) const;
+};
+
+
+
+template <int dim>
+double Coefficient<dim>::value (const Point<dim> &p,
+				const unsigned int) const 
+{
+  if (p.square() < 0.5*0.5)
+    return 20;
+  else
+    return 1;
+};
+
+
+
+template <int dim>
+void Coefficient<dim>::value_list (const vector<Point<dim> > &points,
+				   vector<double>            &values,
+				   const unsigned int component) const 
+{
+  const unsigned int n_points = points.size();
+
+  Assert (values.size() == n_points, 
+	  ExcVectorHasWrongSize (values.size(), n_points));
+  
+  Assert (component == 0, 
+	  ExcWrongComponent (component, 1));
+  
+  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 () :
+		dof_handler (triangulation)
+{};
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::setup_system ()
+{
+  dof_handler.distribute_dofs (fe);
+
+				   //...
+  hanging_node_constraints.clear ();
+  DoFTools::make_hanging_node_constraints (dof_handler,
+					   hanging_node_constraints);
+  hanging_node_constraints.close ();
+
+  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);
+
+				   //...
+  hanging_node_constraints.condense (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 () 
+{  
+  const Coefficient<dim> coefficient;
+
+  QGauss3<dim>  quadrature_formula;
+
+  FEValues<dim> fe_values (fe, quadrature_formula, 
+			   UpdateFlags(update_values    |
+				       update_gradients |
+				       update_q_points  |
+				       update_JxW_values));
+
+  const unsigned int dofs_per_cell = fe.dofs_per_cell;
+  const unsigned int n_q_points    = quadrature_formula.n_quadrature_points;
+
+  FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
+  Vector<double>     cell_rhs (dofs_per_cell);
+
+  vector<int>        local_dof_indices (dofs_per_cell);
+
+  vector<double>     coefficient_values (n_q_points);
+
+  DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
+					endc = dof_handler.end();
+  for (; cell!=endc; ++cell)
+    {
+      cell_matrix.clear ();
+      cell_rhs.clear ();
+
+      fe_values.reinit (cell);
+      const FullMatrix<double> 
+	& shape_values = fe_values.get_shape_values();
+      const vector<vector<Tensor<1,dim> > >
+	& shape_grads  = fe_values.get_shape_grads();
+      const vector<double>
+	& JxW_values   = fe_values.get_JxW_values();
+      const vector<Point<dim> >
+	& q_points     = fe_values.get_quadrature_points();
+
+      coefficient.value_list (q_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] *
+				   (shape_grads[i][q_point]    *
+				    shape_grads[j][q_point])   *
+				   JxW_values[q_point]);
+
+	    cell_rhs(i) += (shape_values (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);
+	};
+    };
+
+  map<int,double> boundary_values;
+  VectorTools::interpolate_boundary_values (dof_handler,
+					    0,
+					    ZeroFunction<dim>(),
+					    boundary_values);
+  MatrixTools<dim>::apply_boundary_values (boundary_values,
+					   system_matrix,
+					   solution,
+					   system_rhs);
+
+  hanging_node_constraints.condense (system_matrix);
+  hanging_node_constraints.condense (system_rhs);
+};
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::solve () 
+{
+  SolverControl           solver_control (1000, 1e-12);
+  PrimitiveVectorMemory<> vector_memory;
+  SolverCG<>              cg (solver_control, vector_memory);
+
+  PreconditionRelaxation<>
+    preconditioner(system_matrix,
+		   &SparseMatrix<double>::template precondition_SSOR<double>,
+		   1.2);
+
+  cg.solve (system_matrix, solution, system_rhs,
+	    preconditioner);
+
+  hanging_node_constraints.distribute (solution);
+};
+
+
+				 //...
+template <int dim>
+void LaplaceProblem<dim>::refine_grid ()
+{
+  Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+
+  KellyErrorEstimator<dim>::FunctionMap neumann_boundary;
+  KellyErrorEstimator<dim>::estimate (dof_handler,
+				      QGauss3<dim-1>(),
+				      neumann_boundary,
+				      solution,
+				      estimated_error_per_cell);
+
+  triangulation.refine_and_coarsen_fixed_number (estimated_error_per_cell,
+						 0.3, 0.03);
+  triangulation.execute_coarsening_and_refinement ();
+};
+
+template <int dim>
+void LaplaceProblem<dim>::output_results (const unsigned int cycle) const
+{
+				   // ...
+  string filename = "grid-";
+  filename += ('0' + cycle);
+  Assert (cycle < 10, ExcInternalError());
+  
+  filename += ".eps";
+  ofstream output (filename.c_str());
+
+  GridOut grid_out;
+  grid_out.write_eps (triangulation, output);
+};
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::run () 
+{
+  for (unsigned int cycle=0; cycle<8; ++cycle)
+    {
+      cout << "Cycle " << cycle << ':' << endl;
+
+      if (cycle == 0)
+	{
+					   //...
+	  GridGenerator::hyper_ball (triangulation);
+
+	  static const HyperBallBoundary<dim> boundary;
+	  triangulation.set_boundary (0, boundary);
+
+	  triangulation.refine_global (1);
+	}
+      else
+					 // ...
+	refine_grid ();
+
+      cout << "   Number of active cells: "
+	   << triangulation.n_active_cells()
+	   << endl;
+
+      setup_system ();
+      assemble_system ();
+      solve ();
+      output_results (cycle);
+    };
+
+  DataOut<dim>::EpsFlags eps_flags;
+  eps_flags.z_scaling = 4;
+  
+  DataOut<dim> data_out;
+  data_out.set_flags (eps_flags);
+
+  data_out.attach_dof_handler (dof_handler);
+  data_out.add_data_vector (solution, "solution");
+  data_out.build_patches ();
+  
+  ofstream output ("final-solution.eps");
+  data_out.write_eps (output);
+};
+
+    
+
+int main () 
+{
+  deallog.depth_console (0);
+
+  LaplaceProblem<2> laplace_problem_2d;
+  laplace_problem_2d.run ();
+  
+  return 0;
+};
diff --git a/deal.II/examples/step-6/Makefile b/deal.II/examples/step-6/Makefile
new file mode 100644
index 0000000000..dc2798f839
--- /dev/null
+++ b/deal.II/examples/step-6/Makefile
@@ -0,0 +1,94 @@
+# $Id$
+# Copyright W. Bangerth, University of Heidelberg, 1999
+
+# Template for makefiles for the examples subdirectory. In principle,
+# everything should be done automatically if you set the target file
+# here correctly. We get deduce it from the files in the present
+# directory:
+target   = $(basename $(shell echo step-*.cc))
+
+# All dependencies between files should be updated by the included
+# file Makefile.dep if necessary. Object files are compiled into
+# the archives ./Obj.a and ./Obj.g.a. By default, the debug version
+# is used to link. It you don't like that, change the following
+# variable to "off"
+debug-mode = on
+
+
+
+###############################################################################
+# Internals
+
+#deal include base path
+D = ../../../..
+
+include $D/common/Make.global_options
+
+# get lists of files we need
+
+
+# list of libraries needed to link with
+libs     = -ldeal_II_2d  -llac -lbase
+libs.g   = -ldeal_II_2d.g -llac.g -lbase.g
+
+
+# check whether we use debug mode or not
+ifeq ($(debug-mode),on)
+  libraries = $(target).go $(libs.g)
+  flags     = $(CXXFLAGS.g)
+else
+  libraries = $(target).go $(libs)
+  flags     = $(CXXFLAGS.o)
+endif
+
+
+
+# make rule for the target. $^ is the object file $(target).g?o
+$(target) : $(libraries)
+	@echo ============================ Linking $@
+	@$(CXX) $(flags) -o $@ $^
+
+# rule how to run the program
+run: $(target)
+	@echo ============================ Running $@
+	@./$(target)
+
+
+# rule to make object files
+%.go : %.cc
+	@echo ============================ Compiling with debugging information:   $<
+	@$(CXX) $(CXXFLAGS.g) -c $< -o $@
+%.o : %.cc
+	@echo ============================ Compiling with optimization:   $<
+	@$(CXX) $(CXXFLAGS) -c $< -o $@
+
+
+clean:
+	-rm -f *.o *.go *~ Makefile.dep $(target) *gmv *gnuplot *gpl *eps
+
+
+
+.PHONY: clean
+
+
+# Rule to generate the dependency file. This file is
+# automagically remade whenever needed, i.e. whenever
+# one of the cc-/h-files changed. Make detects whether
+# to remake this file upon inclusion at the bottom
+# of this file.
+#
+# Since the script prefixes the output names by lib/g?o, we have to
+# strip that again (the script was written for the main libraries and
+# large projects where object files are put into subdirs)
+Makefile.dep: $(target).cc Makefile \
+              $(shell echo $D/base/include/base/*.h \
+                           $D/lac/include/lac/*.h  \
+                           $D/deal.II/include/*/*.h)
+	@echo ============================ Remaking Makefile
+	@perl $D/common/scripts/make_dependencies.pl  $(INCLUDE) $(target).cc \
+		| perl -pi -e 's!lib/g?o/!!g;' \
+		> Makefile.dep
+
+
+include Makefile.dep
+
diff --git a/deal.II/examples/step-6/step-6.cc b/deal.II/examples/step-6/step-6.cc
new file mode 100644
index 0000000000..b87329d1dc
--- /dev/null
+++ b/deal.II/examples/step-6/step-6.cc
@@ -0,0 +1,349 @@
+/* $Id$ */
+/* Author: Wolfgang Bangerth, University of Heidelberg, 2000 */
+
+
+				 // still unfinished
+#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/vector_memory.h>
+#include <lac/precondition.h>
+#include <grid/tria.h>
+#include <dofs/dof_handler.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+#include <fe/fe_lib.lagrange.h>
+#include <fe/fe_values.h>
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <numerics/data_out.h>
+
+				 // out statt in
+#include <grid/grid_out.h>
+
+#include <grid/tria_boundary_lib.h>
+
+				 //...
+#include <dofs/dof_constraints.h>
+#include <numerics/error_estimator.h>
+
+#include <fstream>
+
+
+
+template <int dim>
+class LaplaceProblem 
+{
+  public:
+    LaplaceProblem ();
+    void run ();
+    
+  private:
+    void setup_system ();
+    void assemble_system ();
+    void solve ();
+    void refine_grid ();
+    void output_results (const unsigned int cycle) const;
+
+    Triangulation<dim>   triangulation;
+    FEQ1<dim>            fe;
+    DoFHandler<dim>      dof_handler;
+
+    ConstraintMatrix     hanging_node_constraints;
+
+    SparseMatrixStruct   sparsity_pattern;
+    SparseMatrix<double> system_matrix;
+
+    Vector<double>       solution;
+    Vector<double>       system_rhs;
+};
+
+
+
+template <int dim>
+class Coefficient : public Function<dim> 
+{
+  public:
+    virtual double value (const Point<dim>   &p,
+			  const unsigned int  component = 0) const;
+    
+    virtual void value_list (const vector<Point<dim> > &points,
+			     vector<double>            &values,
+			     const unsigned int         component = 0) const;
+};
+
+
+
+template <int dim>
+double Coefficient<dim>::value (const Point<dim> &p,
+				const unsigned int) const 
+{
+  if (p.square() < 0.5*0.5)
+    return 20;
+  else
+    return 1;
+};
+
+
+
+template <int dim>
+void Coefficient<dim>::value_list (const vector<Point<dim> > &points,
+				   vector<double>            &values,
+				   const unsigned int component) const 
+{
+  const unsigned int n_points = points.size();
+
+  Assert (values.size() == n_points, 
+	  ExcVectorHasWrongSize (values.size(), n_points));
+  
+  Assert (component == 0, 
+	  ExcWrongComponent (component, 1));
+  
+  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 () :
+		dof_handler (triangulation)
+{};
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::setup_system ()
+{
+  dof_handler.distribute_dofs (fe);
+
+				   //...
+  hanging_node_constraints.clear ();
+  DoFTools::make_hanging_node_constraints (dof_handler,
+					   hanging_node_constraints);
+  hanging_node_constraints.close ();
+
+  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);
+
+				   //...
+  hanging_node_constraints.condense (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 () 
+{  
+  const Coefficient<dim> coefficient;
+
+  QGauss3<dim>  quadrature_formula;
+
+  FEValues<dim> fe_values (fe, quadrature_formula, 
+			   UpdateFlags(update_values    |
+				       update_gradients |
+				       update_q_points  |
+				       update_JxW_values));
+
+  const unsigned int dofs_per_cell = fe.dofs_per_cell;
+  const unsigned int n_q_points    = quadrature_formula.n_quadrature_points;
+
+  FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
+  Vector<double>     cell_rhs (dofs_per_cell);
+
+  vector<int>        local_dof_indices (dofs_per_cell);
+
+  vector<double>     coefficient_values (n_q_points);
+
+  DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
+					endc = dof_handler.end();
+  for (; cell!=endc; ++cell)
+    {
+      cell_matrix.clear ();
+      cell_rhs.clear ();
+
+      fe_values.reinit (cell);
+      const FullMatrix<double> 
+	& shape_values = fe_values.get_shape_values();
+      const vector<vector<Tensor<1,dim> > >
+	& shape_grads  = fe_values.get_shape_grads();
+      const vector<double>
+	& JxW_values   = fe_values.get_JxW_values();
+      const vector<Point<dim> >
+	& q_points     = fe_values.get_quadrature_points();
+
+      coefficient.value_list (q_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] *
+				   (shape_grads[i][q_point]    *
+				    shape_grads[j][q_point])   *
+				   JxW_values[q_point]);
+
+	    cell_rhs(i) += (shape_values (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);
+	};
+    };
+
+  map<int,double> boundary_values;
+  VectorTools::interpolate_boundary_values (dof_handler,
+					    0,
+					    ZeroFunction<dim>(),
+					    boundary_values);
+  MatrixTools<dim>::apply_boundary_values (boundary_values,
+					   system_matrix,
+					   solution,
+					   system_rhs);
+
+  hanging_node_constraints.condense (system_matrix);
+  hanging_node_constraints.condense (system_rhs);
+};
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::solve () 
+{
+  SolverControl           solver_control (1000, 1e-12);
+  PrimitiveVectorMemory<> vector_memory;
+  SolverCG<>              cg (solver_control, vector_memory);
+
+  PreconditionRelaxation<>
+    preconditioner(system_matrix,
+		   &SparseMatrix<double>::template precondition_SSOR<double>,
+		   1.2);
+
+  cg.solve (system_matrix, solution, system_rhs,
+	    preconditioner);
+
+  hanging_node_constraints.distribute (solution);
+};
+
+
+				 //...
+template <int dim>
+void LaplaceProblem<dim>::refine_grid ()
+{
+  Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+
+  KellyErrorEstimator<dim>::FunctionMap neumann_boundary;
+  KellyErrorEstimator<dim>::estimate (dof_handler,
+				      QGauss3<dim-1>(),
+				      neumann_boundary,
+				      solution,
+				      estimated_error_per_cell);
+
+  triangulation.refine_and_coarsen_fixed_number (estimated_error_per_cell,
+						 0.3, 0.03);
+  triangulation.execute_coarsening_and_refinement ();
+};
+
+template <int dim>
+void LaplaceProblem<dim>::output_results (const unsigned int cycle) const
+{
+				   // ...
+  string filename = "grid-";
+  filename += ('0' + cycle);
+  Assert (cycle < 10, ExcInternalError());
+  
+  filename += ".eps";
+  ofstream output (filename.c_str());
+
+  GridOut grid_out;
+  grid_out.write_eps (triangulation, output);
+};
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::run () 
+{
+  for (unsigned int cycle=0; cycle<8; ++cycle)
+    {
+      cout << "Cycle " << cycle << ':' << endl;
+
+      if (cycle == 0)
+	{
+					   //...
+	  GridGenerator::hyper_ball (triangulation);
+
+	  static const HyperBallBoundary<dim> boundary;
+	  triangulation.set_boundary (0, boundary);
+
+	  triangulation.refine_global (1);
+	}
+      else
+					 // ...
+	refine_grid ();
+
+      cout << "   Number of active cells: "
+	   << triangulation.n_active_cells()
+	   << endl;
+
+      setup_system ();
+      assemble_system ();
+      solve ();
+      output_results (cycle);
+    };
+
+  DataOut<dim>::EpsFlags eps_flags;
+  eps_flags.z_scaling = 4;
+  
+  DataOut<dim> data_out;
+  data_out.set_flags (eps_flags);
+
+  data_out.attach_dof_handler (dof_handler);
+  data_out.add_data_vector (solution, "solution");
+  data_out.build_patches ();
+  
+  ofstream output ("final-solution.eps");
+  data_out.write_eps (output);
+};
+
+    
+
+int main () 
+{
+  deallog.depth_console (0);
+
+  LaplaceProblem<2> laplace_problem_2d;
+  laplace_problem_2d.run ();
+  
+  return 0;
+};
-- 
2.39.5