From ffb1cd645640fa1954201e971faac64072943e52 Mon Sep 17 00:00:00 2001
From: wolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Wed, 26 Jan 2000 23:15:27 +0000
Subject: [PATCH] Add step-7.

git-svn-id: https://svn.dealii.org/trunk@2274 0785d39b-7218-0410-832d-ea1e28bc413d
---
 .../examples/step-by-step/step-6/step-6.cc    |   1 +
 .../examples/step-by-step/step-7/Makefile     |  98 ++++
 .../examples/step-by-step/step-7/step-7.cc    | 486 ++++++++++++++++++
 .../step-7.data/intro.html                    |   0
 .../step-7.data/results.html                  |   0
 deal.II/examples/step-6/step-6.cc             |   1 +
 deal.II/examples/step-7/Makefile              |  98 ++++
 deal.II/examples/step-7/step-7.cc             | 486 ++++++++++++++++++
 8 files changed, 1170 insertions(+)
 create mode 100644 deal.II/deal.II/Attic/examples/step-by-step/step-7/Makefile
 create mode 100644 deal.II/deal.II/Attic/examples/step-by-step/step-7/step-7.cc
 create mode 100644 deal.II/doc/tutorial/chapter-2.step-by-step/step-7.data/intro.html
 create mode 100644 deal.II/doc/tutorial/chapter-2.step-by-step/step-7.data/results.html
 create mode 100644 deal.II/examples/step-7/Makefile
 create mode 100644 deal.II/examples/step-7/step-7.cc

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
index aa281775ed..2f33358c23 100644
--- 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
@@ -811,6 +811,7 @@ void LaplaceProblem<dim>::output_results (const unsigned int cycle) const
 };
 
 
+
 template <int dim>
 void LaplaceProblem<dim>::run () 
 {
diff --git a/deal.II/deal.II/Attic/examples/step-by-step/step-7/Makefile b/deal.II/deal.II/Attic/examples/step-by-step/step-7/Makefile
new file mode 100644
index 0000000000..be82c25ebd
--- /dev/null
+++ b/deal.II/deal.II/Attic/examples/step-by-step/step-7/Makefile
@@ -0,0 +1,98 @@
+# $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.g   = $(lib-deal2-2d.g) \
+	   $(lib-lac.g)      \
+           $(lib-base.g)
+libs     = $(lib-deal2-2d.o) \
+	   $(lib-lac.o)      \
+           $(lib-base.o)
+
+
+# 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-7/step-7.cc b/deal.II/deal.II/Attic/examples/step-by-step/step-7/step-7.cc
new file mode 100644
index 0000000000..c14b80d9ff
--- /dev/null
+++ b/deal.II/deal.II/Attic/examples/step-by-step/step-7/step-7.cc
@@ -0,0 +1,486 @@
+/* $Id$ */
+/* Author: Wolfgang Bangerth, University of Heidelberg, 2000 */
+
+#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 <grid/tria_boundary_lib.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+#include <fe/fe_values.h>
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <fe/fe_lib.lagrange.h>
+#include <dofs/dof_constraints.h>
+#include <numerics/error_estimator.h>
+
+#include <numerics/dof_renumbering.h>
+#include <base/smartpointer.h>
+
+#include <fstream>
+
+
+template <int dim>
+class LaplaceProblem 
+{
+  public:
+    enum RefinementMode {
+	  global_refinement, adaptive_refinement
+    };
+    
+    LaplaceProblem (const FiniteElement<dim> &fe,
+		    const RefinementMode      refinement_mode);
+    ~LaplaceProblem ();
+    void run ();
+    
+  private:
+    void setup_system ();
+    void assemble_system ();
+    void solve ();
+    void refine_grid ();
+    void process_solution (const unsigned int cycle) const;
+
+    Triangulation<dim>                      triangulation;
+    DoFHandler<dim>                         dof_handler;
+				     //...
+    SmartPointer<const FiniteElement<dim> > fe;
+    ConstraintMatrix                        hanging_node_constraints;
+
+    SparsityPattern                         sparsity_pattern;
+    SparseMatrix<double>                    system_matrix;
+
+    Vector<double>                          solution;
+    Vector<double>                          system_rhs;
+
+    RefinementMode                          refinement_mode;
+};
+
+
+
+template <int dim>
+class SolutionBase 
+{
+  protected:
+    static const unsigned int n_source_centers = 3;    
+    static const Point<dim>   source_centers[n_source_centers];
+    static const double       width;
+};
+
+
+template <int dim>
+class Solution : public Function<dim>,
+		 protected SolutionBase<dim>
+{
+  public:
+    virtual double value (const Point<dim>   &p,
+			  const unsigned int  component = 0) const;
+    virtual Tensor<1,dim> gradient (const Point<dim>   &p,
+				    const unsigned int  component = 0) const;
+};
+
+
+
+template <int dim>
+class RightHandSide : public Function<dim>,
+		      protected SolutionBase<dim>
+{
+  public:
+    virtual double value (const Point<dim>   &p,
+			  const unsigned int  component = 0) const;
+};
+
+
+
+template <>
+const Point<1>
+SolutionBase<1>::source_centers[SolutionBase<1>::n_source_centers]
+= { Point<1>(-1.0 / 3.0), 
+    Point<1>(0.0), 
+    Point<1>(+1.0 / 3.0)   };
+
+template <>
+const Point<2>
+SolutionBase<2>::source_centers[SolutionBase<2>::n_source_centers]
+= { Point<2>(-0.5, +0.5), 
+    Point<2>(-0.5, -0.5), 
+    Point<2>(+0.5, -0.5)   };
+
+template <int dim>
+const double SolutionBase<dim>::width = 0.15;
+
+
+
+template <int dim>
+double Solution<dim>::value (const Point<dim>   &p,
+			     const unsigned int) const
+{
+  double return_value = 0;
+  for (unsigned int i=0; i<n_source_centers; ++i)
+    {
+      const Point<dim> shifted_point = p-source_centers[i];
+      
+      return_value += exp(-shifted_point.square() / (width*width));
+    };
+  
+  return return_value;
+};
+
+
+
+template <int dim>
+Tensor<1,dim> Solution<dim>::gradient (const Point<dim>   &p,
+				       const unsigned int) const
+{
+  Tensor<1,dim> return_value;
+  for (unsigned int i=0; i<n_source_centers; ++i)
+    {
+      const Point<dim> shifted_point = p-source_centers[i];
+      
+      return_value += (-2 / (width*width) *
+		       exp(-shifted_point.square() / (width*width)) *
+		       shifted_point);
+    };
+  
+  return return_value;
+};
+
+
+
+template <int dim>
+double RightHandSide<dim>::value (const Point<dim>   &p,
+				  const unsigned int) const
+{
+  double return_value = 0;
+  for (unsigned int i=0; i<n_source_centers; ++i)
+    {
+      const Point<dim> shifted_point = p-source_centers[i];
+      
+      return_value += ((2*dim - 4*shifted_point.square()/(width*width)) / (width*width) *
+		       exp(-shifted_point.square() / (width*width)));
+    };
+  
+  return return_value;
+};
+
+
+
+template <int dim>
+LaplaceProblem<dim>::LaplaceProblem (const FiniteElement<dim> &fe,
+				     const RefinementMode refinement_mode) :
+		dof_handler (triangulation),
+		fe (&fe),
+		refinement_mode (refinement_mode)
+{};
+
+
+
+template <int dim>
+LaplaceProblem<dim>::~LaplaceProblem () 
+{
+  dof_handler.clear ();
+};
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::setup_system ()
+{
+  dof_handler.distribute_dofs (*fe);
+				   // Renumber the degrees of freedom...
+  DoFRenumbering::Cuthill_McKee (dof_handler);
+
+  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 () 
+{  
+  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;
+
+  RightHandSide<dim> right_hand_side;
+  vector<double>     rhs_values (n_q_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);
+
+  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();
+
+      right_hand_side.value_list (q_points, rhs_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) += (shape_grads[i][q_point] *
+				   shape_grads[j][q_point] *
+				   JxW_values[q_point]);
+
+	    cell_rhs(i) += (shape_values (i,q_point) *
+			    rhs_values [q_point] *
+			    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);
+	};
+    };
+
+  hanging_node_constraints.condense (system_matrix);
+  hanging_node_constraints.condense (system_rhs);
+
+  map<int,double> boundary_values;
+  VectorTools::interpolate_boundary_values (dof_handler,
+					    0,
+					    Solution<dim>(),
+					    boundary_values);
+  MatrixTools<dim>::apply_boundary_values (boundary_values,
+					   system_matrix,
+					   solution,
+					   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 ()
+{
+  switch (refinement_mode) 
+    {
+      case global_refinement:
+      {
+	triangulation.refine_global (1);
+	break;
+      };
+       
+      case adaptive_refinement:
+      {
+	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 ();
+
+	break;
+      };
+    };
+};
+
+
+#include <numerics/data_out.h>
+
+template <int dim>
+void LaplaceProblem<dim>::process_solution (const unsigned int cycle) const
+{
+  Vector<float> difference_per_cell (triangulation.n_active_cells());
+  
+  VectorTools::integrate_difference (dof_handler,
+				     solution,
+				     Solution<dim>(),
+				     difference_per_cell,
+				     QGauss3<dim>(),
+				     L2_norm);
+  const double L2_error = difference_per_cell.l2_norm();
+
+  VectorTools::integrate_difference (dof_handler,
+				     solution,
+				     Solution<dim>(),
+				     difference_per_cell,
+				     QGauss3<dim>(),
+				     H1_seminorm);
+  const double H1_error = difference_per_cell.l2_norm();
+
+  VectorTools::integrate_difference (dof_handler,
+				     solution,
+				     Solution<dim>(),
+				     difference_per_cell,
+				     QGauss3<dim>(),
+				     Linfty_norm);
+  const double Linfty_error = difference_per_cell.linfty_norm();
+  
+  cout << "Cycle " << cycle << ':' 
+       << endl
+       << "   Number of active cells:       "
+       << triangulation.n_active_cells()
+       << endl
+       << "   Number of degrees of freedom: "
+       << dof_handler.n_dofs()
+       << endl;
+
+  cout << "   L2     error: " << L2_error      << endl
+       << "   H1     error: " << H1_error      << endl
+       << "   Linfty error: " << Linfty_error  << endl;
+};
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::run () 
+{
+  for (unsigned int cycle=0; cycle<12; ++cycle)
+    {
+      if (cycle == 0)
+	{
+	  GridGenerator::hyper_cube (triangulation, -1, 1);
+	  triangulation.refine_global (1);
+	}
+      else
+	refine_grid ();      
+
+      setup_system ();
+      
+      assemble_system ();
+      solve ();
+      process_solution (cycle);
+    };
+  
+  string filename;
+  switch (refinement_mode)
+    {
+      case global_refinement:
+	    filename = "solution-global";
+	    break;
+      case adaptive_refinement:
+	    filename = "solution-adaptive";
+	    break;
+      default:
+	    Assert (false, ExcInternalError());
+    };
+  filename += ".gmv";
+	    
+  ofstream output (filename.c_str());
+
+
+  DataOut<dim> data_out;
+  data_out.attach_dof_handler (dof_handler);
+  data_out.add_data_vector (solution, "solution");
+  data_out.build_patches ();
+  data_out.write_gmv (output);
+};
+
+
+
+int main () 
+{
+  try
+    {
+      deallog.depth_console (0);
+
+      FEQ1<2> fe;
+      LaplaceProblem<2> laplace_problem_2d (fe, LaplaceProblem<2>::adaptive_refinement);
+      laplace_problem_2d.run ();
+    }
+  catch (exception &exc)
+    {
+      cerr << endl << endl
+	   << "----------------------------------------------------"
+	   << endl;
+      cerr << "Exception on processing: " << endl
+	   << exc.what() << endl
+	   << "Aborting!" << endl
+	   << "----------------------------------------------------"
+	   << endl;
+      return 1;
+    }
+  catch (...) 
+    {
+      cerr << endl << endl
+	   << "----------------------------------------------------"
+	   << endl;
+      cerr << "Unknown exception!" << endl
+	   << "Aborting!" << endl
+	   << "----------------------------------------------------"
+	   << endl;
+      return 1;
+    };
+
+  return 0;
+};
diff --git a/deal.II/doc/tutorial/chapter-2.step-by-step/step-7.data/intro.html b/deal.II/doc/tutorial/chapter-2.step-by-step/step-7.data/intro.html
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/deal.II/doc/tutorial/chapter-2.step-by-step/step-7.data/results.html b/deal.II/doc/tutorial/chapter-2.step-by-step/step-7.data/results.html
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/deal.II/examples/step-6/step-6.cc b/deal.II/examples/step-6/step-6.cc
index aa281775ed..2f33358c23 100644
--- a/deal.II/examples/step-6/step-6.cc
+++ b/deal.II/examples/step-6/step-6.cc
@@ -811,6 +811,7 @@ void LaplaceProblem<dim>::output_results (const unsigned int cycle) const
 };
 
 
+
 template <int dim>
 void LaplaceProblem<dim>::run () 
 {
diff --git a/deal.II/examples/step-7/Makefile b/deal.II/examples/step-7/Makefile
new file mode 100644
index 0000000000..be82c25ebd
--- /dev/null
+++ b/deal.II/examples/step-7/Makefile
@@ -0,0 +1,98 @@
+# $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.g   = $(lib-deal2-2d.g) \
+	   $(lib-lac.g)      \
+           $(lib-base.g)
+libs     = $(lib-deal2-2d.o) \
+	   $(lib-lac.o)      \
+           $(lib-base.o)
+
+
+# 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-7/step-7.cc b/deal.II/examples/step-7/step-7.cc
new file mode 100644
index 0000000000..c14b80d9ff
--- /dev/null
+++ b/deal.II/examples/step-7/step-7.cc
@@ -0,0 +1,486 @@
+/* $Id$ */
+/* Author: Wolfgang Bangerth, University of Heidelberg, 2000 */
+
+#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 <grid/tria_boundary_lib.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+#include <fe/fe_values.h>
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <fe/fe_lib.lagrange.h>
+#include <dofs/dof_constraints.h>
+#include <numerics/error_estimator.h>
+
+#include <numerics/dof_renumbering.h>
+#include <base/smartpointer.h>
+
+#include <fstream>
+
+
+template <int dim>
+class LaplaceProblem 
+{
+  public:
+    enum RefinementMode {
+	  global_refinement, adaptive_refinement
+    };
+    
+    LaplaceProblem (const FiniteElement<dim> &fe,
+		    const RefinementMode      refinement_mode);
+    ~LaplaceProblem ();
+    void run ();
+    
+  private:
+    void setup_system ();
+    void assemble_system ();
+    void solve ();
+    void refine_grid ();
+    void process_solution (const unsigned int cycle) const;
+
+    Triangulation<dim>                      triangulation;
+    DoFHandler<dim>                         dof_handler;
+				     //...
+    SmartPointer<const FiniteElement<dim> > fe;
+    ConstraintMatrix                        hanging_node_constraints;
+
+    SparsityPattern                         sparsity_pattern;
+    SparseMatrix<double>                    system_matrix;
+
+    Vector<double>                          solution;
+    Vector<double>                          system_rhs;
+
+    RefinementMode                          refinement_mode;
+};
+
+
+
+template <int dim>
+class SolutionBase 
+{
+  protected:
+    static const unsigned int n_source_centers = 3;    
+    static const Point<dim>   source_centers[n_source_centers];
+    static const double       width;
+};
+
+
+template <int dim>
+class Solution : public Function<dim>,
+		 protected SolutionBase<dim>
+{
+  public:
+    virtual double value (const Point<dim>   &p,
+			  const unsigned int  component = 0) const;
+    virtual Tensor<1,dim> gradient (const Point<dim>   &p,
+				    const unsigned int  component = 0) const;
+};
+
+
+
+template <int dim>
+class RightHandSide : public Function<dim>,
+		      protected SolutionBase<dim>
+{
+  public:
+    virtual double value (const Point<dim>   &p,
+			  const unsigned int  component = 0) const;
+};
+
+
+
+template <>
+const Point<1>
+SolutionBase<1>::source_centers[SolutionBase<1>::n_source_centers]
+= { Point<1>(-1.0 / 3.0), 
+    Point<1>(0.0), 
+    Point<1>(+1.0 / 3.0)   };
+
+template <>
+const Point<2>
+SolutionBase<2>::source_centers[SolutionBase<2>::n_source_centers]
+= { Point<2>(-0.5, +0.5), 
+    Point<2>(-0.5, -0.5), 
+    Point<2>(+0.5, -0.5)   };
+
+template <int dim>
+const double SolutionBase<dim>::width = 0.15;
+
+
+
+template <int dim>
+double Solution<dim>::value (const Point<dim>   &p,
+			     const unsigned int) const
+{
+  double return_value = 0;
+  for (unsigned int i=0; i<n_source_centers; ++i)
+    {
+      const Point<dim> shifted_point = p-source_centers[i];
+      
+      return_value += exp(-shifted_point.square() / (width*width));
+    };
+  
+  return return_value;
+};
+
+
+
+template <int dim>
+Tensor<1,dim> Solution<dim>::gradient (const Point<dim>   &p,
+				       const unsigned int) const
+{
+  Tensor<1,dim> return_value;
+  for (unsigned int i=0; i<n_source_centers; ++i)
+    {
+      const Point<dim> shifted_point = p-source_centers[i];
+      
+      return_value += (-2 / (width*width) *
+		       exp(-shifted_point.square() / (width*width)) *
+		       shifted_point);
+    };
+  
+  return return_value;
+};
+
+
+
+template <int dim>
+double RightHandSide<dim>::value (const Point<dim>   &p,
+				  const unsigned int) const
+{
+  double return_value = 0;
+  for (unsigned int i=0; i<n_source_centers; ++i)
+    {
+      const Point<dim> shifted_point = p-source_centers[i];
+      
+      return_value += ((2*dim - 4*shifted_point.square()/(width*width)) / (width*width) *
+		       exp(-shifted_point.square() / (width*width)));
+    };
+  
+  return return_value;
+};
+
+
+
+template <int dim>
+LaplaceProblem<dim>::LaplaceProblem (const FiniteElement<dim> &fe,
+				     const RefinementMode refinement_mode) :
+		dof_handler (triangulation),
+		fe (&fe),
+		refinement_mode (refinement_mode)
+{};
+
+
+
+template <int dim>
+LaplaceProblem<dim>::~LaplaceProblem () 
+{
+  dof_handler.clear ();
+};
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::setup_system ()
+{
+  dof_handler.distribute_dofs (*fe);
+				   // Renumber the degrees of freedom...
+  DoFRenumbering::Cuthill_McKee (dof_handler);
+
+  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 () 
+{  
+  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;
+
+  RightHandSide<dim> right_hand_side;
+  vector<double>     rhs_values (n_q_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);
+
+  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();
+
+      right_hand_side.value_list (q_points, rhs_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) += (shape_grads[i][q_point] *
+				   shape_grads[j][q_point] *
+				   JxW_values[q_point]);
+
+	    cell_rhs(i) += (shape_values (i,q_point) *
+			    rhs_values [q_point] *
+			    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);
+	};
+    };
+
+  hanging_node_constraints.condense (system_matrix);
+  hanging_node_constraints.condense (system_rhs);
+
+  map<int,double> boundary_values;
+  VectorTools::interpolate_boundary_values (dof_handler,
+					    0,
+					    Solution<dim>(),
+					    boundary_values);
+  MatrixTools<dim>::apply_boundary_values (boundary_values,
+					   system_matrix,
+					   solution,
+					   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 ()
+{
+  switch (refinement_mode) 
+    {
+      case global_refinement:
+      {
+	triangulation.refine_global (1);
+	break;
+      };
+       
+      case adaptive_refinement:
+      {
+	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 ();
+
+	break;
+      };
+    };
+};
+
+
+#include <numerics/data_out.h>
+
+template <int dim>
+void LaplaceProblem<dim>::process_solution (const unsigned int cycle) const
+{
+  Vector<float> difference_per_cell (triangulation.n_active_cells());
+  
+  VectorTools::integrate_difference (dof_handler,
+				     solution,
+				     Solution<dim>(),
+				     difference_per_cell,
+				     QGauss3<dim>(),
+				     L2_norm);
+  const double L2_error = difference_per_cell.l2_norm();
+
+  VectorTools::integrate_difference (dof_handler,
+				     solution,
+				     Solution<dim>(),
+				     difference_per_cell,
+				     QGauss3<dim>(),
+				     H1_seminorm);
+  const double H1_error = difference_per_cell.l2_norm();
+
+  VectorTools::integrate_difference (dof_handler,
+				     solution,
+				     Solution<dim>(),
+				     difference_per_cell,
+				     QGauss3<dim>(),
+				     Linfty_norm);
+  const double Linfty_error = difference_per_cell.linfty_norm();
+  
+  cout << "Cycle " << cycle << ':' 
+       << endl
+       << "   Number of active cells:       "
+       << triangulation.n_active_cells()
+       << endl
+       << "   Number of degrees of freedom: "
+       << dof_handler.n_dofs()
+       << endl;
+
+  cout << "   L2     error: " << L2_error      << endl
+       << "   H1     error: " << H1_error      << endl
+       << "   Linfty error: " << Linfty_error  << endl;
+};
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::run () 
+{
+  for (unsigned int cycle=0; cycle<12; ++cycle)
+    {
+      if (cycle == 0)
+	{
+	  GridGenerator::hyper_cube (triangulation, -1, 1);
+	  triangulation.refine_global (1);
+	}
+      else
+	refine_grid ();      
+
+      setup_system ();
+      
+      assemble_system ();
+      solve ();
+      process_solution (cycle);
+    };
+  
+  string filename;
+  switch (refinement_mode)
+    {
+      case global_refinement:
+	    filename = "solution-global";
+	    break;
+      case adaptive_refinement:
+	    filename = "solution-adaptive";
+	    break;
+      default:
+	    Assert (false, ExcInternalError());
+    };
+  filename += ".gmv";
+	    
+  ofstream output (filename.c_str());
+
+
+  DataOut<dim> data_out;
+  data_out.attach_dof_handler (dof_handler);
+  data_out.add_data_vector (solution, "solution");
+  data_out.build_patches ();
+  data_out.write_gmv (output);
+};
+
+
+
+int main () 
+{
+  try
+    {
+      deallog.depth_console (0);
+
+      FEQ1<2> fe;
+      LaplaceProblem<2> laplace_problem_2d (fe, LaplaceProblem<2>::adaptive_refinement);
+      laplace_problem_2d.run ();
+    }
+  catch (exception &exc)
+    {
+      cerr << endl << endl
+	   << "----------------------------------------------------"
+	   << endl;
+      cerr << "Exception on processing: " << endl
+	   << exc.what() << endl
+	   << "Aborting!" << endl
+	   << "----------------------------------------------------"
+	   << endl;
+      return 1;
+    }
+  catch (...) 
+    {
+      cerr << endl << endl
+	   << "----------------------------------------------------"
+	   << endl;
+      cerr << "Unknown exception!" << endl
+	   << "Aborting!" << endl
+	   << "----------------------------------------------------"
+	   << endl;
+      return 1;
+    };
+
+  return 0;
+};
-- 
2.39.5