]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add a testcase for the multigrid components.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 28 Aug 1998 09:08:04 +0000 (09:08 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 28 Aug 1998 09:08:04 +0000 (09:08 +0000)
git-svn-id: https://svn.dealii.org/trunk@531 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/Attic/examples/Makefile
deal.II/deal.II/Attic/examples/multigrid/Makefile [new file with mode: 0644]
deal.II/deal.II/Attic/examples/multigrid/make_ps [new file with mode: 0644]
tests/big-tests/Makefile
tests/big-tests/multigrid/Makefile [new file with mode: 0644]
tests/big-tests/multigrid/convergence.cc [new file with mode: 0644]
tests/big-tests/multigrid/make_ps [new file with mode: 0644]

index 3c2149f907a042cdef498056cecc21841bf10861..1a75c07be4ce1f7a55cef6c9dbbbe3e3b149cccf 100644 (file)
@@ -3,7 +3,7 @@
 
 
 # list the directories we want to visit
-subdirs = grid/ dof/ poisson/ convergence/ error-estimation/
+subdirs = grid/ dof/ poisson/ convergence/ error-estimation/ multigrid/
 
 # define lists of targets: for each directory we produce a target name
 # for compilation, running and cleaning by appending the action to
diff --git a/deal.II/deal.II/Attic/examples/multigrid/Makefile b/deal.II/deal.II/Attic/examples/multigrid/Makefile
new file mode 100644 (file)
index 0000000..022319d
--- /dev/null
@@ -0,0 +1,118 @@
+# $Id$
+# Copyright W. Bangerth, University of Heidelberg, 1998
+
+# Template for makefiles for the examples subdirectory. In principle,
+# everything should be done automatically if you set the target file
+# here correctly:
+target   = convergence
+
+# 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 = off
+
+# If you want your program to be linked with extra object or library
+# files, specify them here:
+user-libs =
+
+# To run the program, use "make run"; to give parameters to the program,
+# give the parameters to the following variable:
+run-parameters  = 
+
+# To execute additional action apart from running the program, fill
+# in this list:
+additional-run-action = gnuplot make_ps
+
+# To specify which files are to be deleted by "make clean" (apart from
+# the usual ones: object files, executables, backups, etc), fill in the
+# following list
+delete-files = *gnuplot *.eps *ucd *history
+
+
+
+
+###############################################################################
+# Internals
+
+include ../../Make.global_options
+
+
+# get lists of files we need
+cc-files    = $(wildcard *.cc)
+o-files     = $(cc-files:.cc=.o)
+go-files    = $(cc-files:.cc=.go)
+h-files     = $(wildcard *.h)
+lib-h-files = $(wildcard ../../include/*/*.h)
+
+# list of libraries needed to link with
+libs     = ./Obj.a   ../../lib/libdeal_II.a   ../../../lac/lib/liblac.a
+libs.g   = ./Obj.g.a ../../lib/libdeal_II.g.a ../../../lac/lib/liblac.g.a
+
+
+# check whether we use debug mode or not
+ifeq ($(debug-mode),on)
+libraries = $(libs.g)
+flags     = $(CXXFLAGS.g)
+endif
+
+ifeq ($(debug-mode),off)
+libraries = $(libs)
+flags     = $(CXXFLAGS)
+endif
+
+
+
+# make rule for the target
+$(target) : $(libraries)
+       @echo ============================ Linking $@
+       @$(CXX) $(flags) -o $@ $(libraries) $(user-libs)
+
+# rule how to run the program
+run: $(target)
+       $(target) $(run-parameters)
+       $(additional-run-action)
+
+
+# rule to make object files
+%.go : %.cc
+       @echo ============================ Compiling with debugging information:   $<
+       @echo $(CXX) ... -c $< -o $@
+       @$(CXX) $(CXXFLAGS.g) -c $< -o $@
+%.o : %.cc
+       @echo ============================ Compiling with optimization:   $<
+       @echo $(CXX) ... -c $< -o $@
+       @$(CXX) $(CXXFLAGS) -c $< -o $@
+
+
+# rules which files the libraries depend upon
+Obj.a: ./Obj.a($(o-files))
+Obj.g.a: ./Obj.g.a($(go-files))
+
+
+clean:
+       rm -f *.o *.go *~ Makefile.dep Obj.a Obj.g.a $(target) $(delete-files)
+
+
+
+.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.
+#
+#use perl to generate rules for the .go files as well
+#as to make rules not for tria.o and the like, but
+#rather for libnumerics.a(tria.o)
+Makefile.dep: $(cc-files) $(h-files) $(lib-h-files)
+       @echo ============================ Remaking Makefile
+       @perl ../../Make_dep.pl ./Obj $(INCLUDE) $(cc-files) \
+               > Makefile.dep
+
+
+include Makefile.dep
+
diff --git a/deal.II/deal.II/Attic/examples/multigrid/make_ps b/deal.II/deal.II/Attic/examples/multigrid/make_ps
new file mode 100644 (file)
index 0000000..c2df845
--- /dev/null
@@ -0,0 +1,52 @@
+set term postscript eps
+set xlabel "Number of degrees of freedom"
+set data style linespoints
+set logscale xy
+
+
+
+set ylabel "Error"
+
+set output "criss-cross.eps"
+
+plot "criss-cross.history" using 1:2 title "L1 error","criss-cross.history" using 1:3 title "L2 error","criss-cross.history" using 1:4 title "Linfty error","criss-cross.history" using 1:5 title "H1 seminorm error","criss-cross.history" using 1:6 title "H1 error"
+
+
+
+set output "linear.eps"
+
+plot "linear.history" using 1:2 title "L1 error","linear.history" using 1:3 title "L2 error","linear.history" using 1:4 title "Linfty error","linear.history" using 1:5 title "H1 seminorm error","linear.history" using 1:6 title "H1 error"
+
+
+
+set output "quadratic.eps"
+
+plot "quadratic.history" using 1:2 title "L1 error","quadratic.history" using 1:3 title "L2 error","quadratic.history" using 1:4 title "Linfty error","quadratic.history" using 1:5 title "H1 seminorm error","quadratic.history" using 1:6 title "H1 error"
+
+
+
+set output "cubic.eps"
+
+plot "cubic.history" using 1:2 title "L1 error","cubic.history" using 1:3 title "L2 error","cubic.history" using 1:4 title "Linfty error","cubic.history" using 1:5 title "H1 seminorm error","cubic.history" using 1:6 title "H1 error"
+
+
+
+set output "quartic.eps"
+
+plot "quartic.history" using 1:2 title "L1 error","quartic.history" using 1:3 title "L2 error","quartic.history" using 1:4 title "Linfty error","quartic.history" using 1:5 title "H1 seminorm error","quartic.history" using 1:6 title "H1 error"
+
+
+
+set output "l2error.eps"
+set ylabel "L2-error"
+
+plot "criss-cross.history" using 1:3 title "Criss-cross elements", "linear.history" using 1:3 title "Linear elements", "quadratic.history" using 1:3 title "Quadratic elements", "cubic.history" using 1:3 title "Cubic elements", "quartic.history" using 1:3 title "Quartic elements"
+
+
+
+set output "h1error.eps"
+set ylabel "H1-error"
+
+plot "criss-cross.history" using 1:6 title "Criss-cross elements", "linear.history" using 1:6 title "Linear elements", "quadratic.history" using 1:6 title "Quadratic elements", "cubic.history" using 1:6 title "Cubic elements", "quartic.history" using 1:6 title "Quartic elements"
+
+
index 3c2149f907a042cdef498056cecc21841bf10861..1a75c07be4ce1f7a55cef6c9dbbbe3e3b149cccf 100644 (file)
@@ -3,7 +3,7 @@
 
 
 # list the directories we want to visit
-subdirs = grid/ dof/ poisson/ convergence/ error-estimation/
+subdirs = grid/ dof/ poisson/ convergence/ error-estimation/ multigrid/
 
 # define lists of targets: for each directory we produce a target name
 # for compilation, running and cleaning by appending the action to
diff --git a/tests/big-tests/multigrid/Makefile b/tests/big-tests/multigrid/Makefile
new file mode 100644 (file)
index 0000000..022319d
--- /dev/null
@@ -0,0 +1,118 @@
+# $Id$
+# Copyright W. Bangerth, University of Heidelberg, 1998
+
+# Template for makefiles for the examples subdirectory. In principle,
+# everything should be done automatically if you set the target file
+# here correctly:
+target   = convergence
+
+# 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 = off
+
+# If you want your program to be linked with extra object or library
+# files, specify them here:
+user-libs =
+
+# To run the program, use "make run"; to give parameters to the program,
+# give the parameters to the following variable:
+run-parameters  = 
+
+# To execute additional action apart from running the program, fill
+# in this list:
+additional-run-action = gnuplot make_ps
+
+# To specify which files are to be deleted by "make clean" (apart from
+# the usual ones: object files, executables, backups, etc), fill in the
+# following list
+delete-files = *gnuplot *.eps *ucd *history
+
+
+
+
+###############################################################################
+# Internals
+
+include ../../Make.global_options
+
+
+# get lists of files we need
+cc-files    = $(wildcard *.cc)
+o-files     = $(cc-files:.cc=.o)
+go-files    = $(cc-files:.cc=.go)
+h-files     = $(wildcard *.h)
+lib-h-files = $(wildcard ../../include/*/*.h)
+
+# list of libraries needed to link with
+libs     = ./Obj.a   ../../lib/libdeal_II.a   ../../../lac/lib/liblac.a
+libs.g   = ./Obj.g.a ../../lib/libdeal_II.g.a ../../../lac/lib/liblac.g.a
+
+
+# check whether we use debug mode or not
+ifeq ($(debug-mode),on)
+libraries = $(libs.g)
+flags     = $(CXXFLAGS.g)
+endif
+
+ifeq ($(debug-mode),off)
+libraries = $(libs)
+flags     = $(CXXFLAGS)
+endif
+
+
+
+# make rule for the target
+$(target) : $(libraries)
+       @echo ============================ Linking $@
+       @$(CXX) $(flags) -o $@ $(libraries) $(user-libs)
+
+# rule how to run the program
+run: $(target)
+       $(target) $(run-parameters)
+       $(additional-run-action)
+
+
+# rule to make object files
+%.go : %.cc
+       @echo ============================ Compiling with debugging information:   $<
+       @echo $(CXX) ... -c $< -o $@
+       @$(CXX) $(CXXFLAGS.g) -c $< -o $@
+%.o : %.cc
+       @echo ============================ Compiling with optimization:   $<
+       @echo $(CXX) ... -c $< -o $@
+       @$(CXX) $(CXXFLAGS) -c $< -o $@
+
+
+# rules which files the libraries depend upon
+Obj.a: ./Obj.a($(o-files))
+Obj.g.a: ./Obj.g.a($(go-files))
+
+
+clean:
+       rm -f *.o *.go *~ Makefile.dep Obj.a Obj.g.a $(target) $(delete-files)
+
+
+
+.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.
+#
+#use perl to generate rules for the .go files as well
+#as to make rules not for tria.o and the like, but
+#rather for libnumerics.a(tria.o)
+Makefile.dep: $(cc-files) $(h-files) $(lib-h-files)
+       @echo ============================ Remaking Makefile
+       @perl ../../Make_dep.pl ./Obj $(INCLUDE) $(cc-files) \
+               > Makefile.dep
+
+
+include Makefile.dep
+
diff --git a/tests/big-tests/multigrid/convergence.cc b/tests/big-tests/multigrid/convergence.cc
new file mode 100644 (file)
index 0000000..0aa8f56
--- /dev/null
@@ -0,0 +1,511 @@
+/* $Id$ */
+/* Copyright W. Bangerth, University of Heidelberg, 1998 */
+
+
+#include <grid/tria.h>
+#include <grid/dof.h>
+#include <grid/tria_accessor.h>
+#include <grid/dof_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_boundary.h>
+#include <grid/dof_constraints.h>
+#include <basic/function.h>
+#include <basic/data_io.h>
+#include <fe/fe_lib.lagrange.h>
+#include <fe/fe_lib.criss_cross.h>
+#include <fe/quadrature_lib.h>
+#include <numerics/base.h>
+#include <numerics/assembler.h>
+#include <numerics/vectors.h>
+
+#include <map>
+#include <fstream>
+#include <cmath>
+#include <string>
+#include <cstdlib>
+
+
+
+
+
+template <int dim>
+class PoissonEquation :  public Equation<dim> {
+  public:
+    PoissonEquation (const Function<dim> &rhs) :
+                   Equation<dim>(1),
+                   right_hand_side (rhs)  {};
+
+    virtual void assemble (dFMatrix            &cell_matrix,
+                          dVector             &rhs,
+                          const FEValues<dim> &fe_values,
+                          const Triangulation<dim>::cell_iterator &cell) const;
+    virtual void assemble (dFMatrix            &cell_matrix,
+                          const FEValues<dim> &fe_values,
+                          const Triangulation<dim>::cell_iterator &cell) const;
+    virtual void assemble (dVector             &rhs,
+                          const FEValues<dim> &fe_values,
+                          const Triangulation<dim>::cell_iterator &cell) const;
+  protected:
+    const Function<dim> &right_hand_side;
+};
+
+
+
+
+
+
+template <int dim>
+class PoissonProblem : public ProblemBase<dim> {
+  public:
+    PoissonProblem (unsigned int order);
+
+    void clear ();
+    void create_new ();
+    int run (unsigned int level);
+    void print_history (string filename) const;
+    
+  protected:
+    Triangulation<dim> *tria;
+    DoFHandler<dim>    *dof;
+    
+    Function<dim>      *rhs;
+    Function<dim>      *boundary_values;
+
+    vector<double> l1_error, l2_error, linfty_error, h1_seminorm_error, h1_error;
+    vector<int>    n_dofs;
+
+    unsigned int        order;
+};
+
+
+
+
+
+/**
+  Right hand side constructed such that the exact solution is
+  $sin(2 pi x) + sin(2 pi y)$
+  */
+template <int dim>
+class RHSPoly : public Function<dim> {
+  public:
+                                    /**
+                                     * Return the value of the function
+                                     * at the given point.
+                                     */
+    virtual double operator () (const Point<dim> &p) const;
+};
+
+
+
+template <int dim>
+class Solution : public Function<dim> {
+  public:
+                                    /**
+                                     * Return the value of the function
+                                     * at the given point.
+                                     */
+    virtual double operator () (const Point<dim> &p) const;
+                                    /**
+                                     * Return the gradient of the function
+                                     * at the given point.
+                                     */
+    virtual Point<dim> gradient (const Point<dim> &p) const;
+};
+
+
+
+
+double RHSPoly<2>::operator () (const Point<2> &p) const {
+  const double x = p(0),
+              y = p(1);
+  const double pi= 3.1415926536;
+  return 4*pi*pi*(sin(2*pi*x)+sin(2*pi*y));
+};
+
+
+
+double Solution<2>::operator () (const Point<2> &p) const {
+  const double x = p(0),
+              y = p(1);
+  const double pi= 3.1415926536;
+  return sin(2*pi*x)+sin(2*pi*y);
+};
+
+
+Point<2> Solution<2>::gradient (const Point<2> &p) const {
+  const double x = p(0),
+              y = p(1);
+  const double pi= 3.1415926536;
+  return Point<2> (2*pi*cos(2*pi*x),
+                  2*pi*cos(2*pi*y));
+};
+
+  
+
+
+
+
+void PoissonEquation<2>::assemble (dFMatrix            &cell_matrix,
+                                  dVector             &rhs,
+                                  const FEValues<2>   &fe_values,
+                                  const Triangulation<2>::cell_iterator &) const {
+  for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+    for (unsigned int i=0; i<fe_values.total_dofs; ++i) 
+      {
+       for (unsigned int j=0; j<fe_values.total_dofs; ++j)
+         cell_matrix(i,j) += (fe_values.shape_grad(i,point) *
+                              fe_values.shape_grad(j,point)) *
+                             fe_values.JxW(point);
+       rhs(i) += fe_values.shape_value(i,point) *
+                 right_hand_side(fe_values.quadrature_point(point)) *
+                 fe_values.JxW(point);
+      };
+};
+
+
+
+template <int dim>
+void PoissonEquation<dim>::assemble (dFMatrix            &,
+                                    const FEValues<dim> &,
+                                    const Triangulation<dim>::cell_iterator &) const {
+  Assert (false, ExcPureVirtualFunctionCalled());
+};
+
+
+
+template <int dim>
+void PoissonEquation<dim>::assemble (dVector             &,
+                                    const FEValues<dim> &,
+                                    const Triangulation<dim>::cell_iterator &) const {
+  Assert (false, ExcPureVirtualFunctionCalled());
+};
+
+
+
+
+
+
+
+
+
+template <int dim>
+PoissonProblem<dim>::PoissonProblem (unsigned int order) :
+               tria(0), dof(0), rhs(0),
+               boundary_values(0), order(order) {};
+
+
+
+
+template <int dim>
+void PoissonProblem<dim>::clear () {
+  if (tria != 0) {
+    delete tria;
+    tria = 0;
+  };
+  
+  if (dof != 0) {
+    delete dof;
+    dof = 0;
+  };
+
+  if (rhs != 0) 
+    {
+      delete rhs;
+      rhs = 0;
+    };
+
+  if (boundary_values != 0) 
+    {
+      delete boundary_values;
+      boundary_values = 0;
+    };
+
+  ProblemBase<dim>::clear ();
+};
+
+
+
+
+template <int dim>
+void PoissonProblem<dim>::create_new () {
+  clear ();
+  
+  tria = new Triangulation<dim>();
+  dof = new DoFHandler<dim> (tria);
+  set_tria_and_dof (tria, dof);
+};
+
+
+
+
+template <int dim>
+int PoissonProblem<dim>::run (const unsigned int level) {
+  create_new ();
+  
+  cout << "Refinement level = " << level
+       << ", using elements of type <";
+  switch (order)
+    {
+      case 0:
+           cout << "criss-cross";
+           break;
+      default:
+           cout << "Lagrange-" << order;
+           break;
+    };
+  cout << ">" << endl;
+  
+  cout << "    Making grid... ";
+  tria->create_hyper_ball ();
+  HyperBallBoundary<dim> boundary_description;
+  tria->set_boundary (&boundary_description);
+  tria->begin_active()->set_refine_flag();
+  (++(++(tria->begin_active())))->set_refine_flag();
+  tria->execute_refinement ();
+  tria->refine_global (level);
+  cout << tria->n_active_cells() << " active cells." << endl;
+
+  rhs             = new RHSPoly<dim>();
+  boundary_values = new Solution<dim> ();
+  
+
+  FiniteElement<dim>   *fe;
+  PoissonEquation<dim>  equation (*rhs);
+  Quadrature<dim>      *quadrature;
+  Quadrature<dim-1>    *boundary_quadrature;
+  switch (order) {
+    case 0:
+         fe         = new FECrissCross<dim>();
+         quadrature = new QCrissCross1<dim>();
+         boundary_quadrature = new QGauss2<dim-1>();
+         break;
+    case 1:
+         fe         = new FELinear<dim>();
+         quadrature = new QGauss3<dim>();
+         boundary_quadrature = new QGauss2<dim-1>();
+         break;
+    case 2:
+         fe         = new FEQuadraticSub<dim>();
+         quadrature = new QGauss4<dim>();
+         boundary_quadrature = new QGauss3<dim-1>();
+         break;
+    case 3:
+         fe         = new FECubicSub<dim>();
+         quadrature = new QGauss5<dim>();
+         boundary_quadrature = new QGauss4<dim-1>();
+         break;
+    case 4:
+         fe         = new FEQuarticSub<dim>();
+         quadrature = new QGauss6<dim>();
+         boundary_quadrature = new QGauss5<dim-1>();
+         break;
+    default:
+         return 100000;
+  };
+  
+  cout << "    Distributing dofs... "; 
+  dof->distribute_dofs (*fe);
+  cout << dof->n_dofs() << " degrees of freedom." << endl;
+  n_dofs.push_back (dof->n_dofs());
+
+  cout << "    Assembling matrices..." << endl;
+  UpdateFlags update_flags = UpdateFlags(update_q_points  | update_gradients |
+                                        update_JxW_values);
+  
+  ProblemBase<dim>::FunctionMap dirichlet_bc;
+  dirichlet_bc[0] = boundary_values;
+  assemble (equation, *quadrature, *fe, update_flags, dirichlet_bc);
+
+  cout << "    Solving..." << endl;
+  solve ();
+
+  Solution<dim> sol;
+  dVector       l1_error_per_cell, l2_error_per_cell, linfty_error_per_cell;
+  dVector       h1_seminorm_error_per_cell, h1_error_per_cell;
+  
+  cout << "    Calculating L1 error... ";
+  VectorTools<dim>::integrate_difference (*dof_handler,
+                                         solution, sol,
+                                         l1_error_per_cell,
+                                         *quadrature, *fe, L1_norm);
+  cout << l1_error_per_cell.l1_norm() << endl;
+  l1_error.push_back (l1_error_per_cell.l1_norm());
+
+  cout << "    Calculating L2 error... ";
+  VectorTools<dim>::integrate_difference (*dof_handler,
+                                         solution, sol,
+                                         l2_error_per_cell,
+                                         *quadrature, *fe, L2_norm);
+  cout << l2_error_per_cell.l2_norm() << endl;
+  l2_error.push_back (l2_error_per_cell.l2_norm());
+
+  cout << "    Calculating L-infinity error... ";
+  VectorTools<dim>::integrate_difference (*dof_handler,
+                                         solution, sol,
+                                         linfty_error_per_cell,
+                                         *quadrature, *fe, Linfty_norm);
+  cout << linfty_error_per_cell.linfty_norm() << endl;
+  linfty_error.push_back (linfty_error_per_cell.linfty_norm());
+  
+  cout << "    Calculating H1-seminorm error... ";
+  VectorTools<dim>::integrate_difference (*dof_handler,
+                                         solution, sol,
+                                         h1_seminorm_error_per_cell,
+                                         *quadrature, *fe, H1_seminorm);
+  cout << h1_seminorm_error_per_cell.l2_norm() << endl;
+  h1_seminorm_error.push_back (h1_seminorm_error_per_cell.l2_norm());
+
+  cout << "    Calculating H1 error... ";
+  VectorTools<dim>::integrate_difference (*dof_handler,
+                                         solution, sol,
+                                         h1_error_per_cell,
+                                         *quadrature, *fe, H1_norm);
+  cout << h1_error_per_cell.l2_norm() << endl;
+  h1_error.push_back (h1_error_per_cell.l2_norm());
+
+  if (dof->n_dofs()<=5000) 
+    {
+      dVector l1_error_per_dof, l2_error_per_dof, linfty_error_per_dof;
+      dVector h1_seminorm_error_per_dof, h1_error_per_dof;
+      dof->distribute_cell_to_dof_vector (l1_error_per_cell, l1_error_per_dof);
+      dof->distribute_cell_to_dof_vector (l2_error_per_cell, l2_error_per_dof);
+      dof->distribute_cell_to_dof_vector (linfty_error_per_cell,
+                                         linfty_error_per_dof);
+      dof->distribute_cell_to_dof_vector (h1_seminorm_error_per_cell,
+                                         h1_seminorm_error_per_dof);
+      dof->distribute_cell_to_dof_vector (h1_error_per_cell, h1_error_per_dof);
+
+//       dVector projected_solution;
+//       ConstraintMatrix constraints;
+//       constraints.close ();
+//       VectorTools<dim>::project (*dof, constraints, *fe,
+//                              StraightBoundary<dim>(), *quadrature, 
+//                              sol, projected_solution, false,
+//                              *boundary_quadrature);
+//       cout << "    Calculating L2 error of projected solution... ";
+//       VectorTools<dim>::integrate_difference (*dof_handler,
+//                                           projected_solution, sol,
+//                                           l2_error_per_cell,
+//                                           *quadrature, *fe, L2_norm);
+//       cout << l2_error_per_cell.l2_norm() << endl;
+
+
+      string filename;
+      filename = ('0'+order);
+      filename += ".";
+      filename += ('0'+level);
+      filename += ".ucd";
+      cout << "    Writing error plots to <" << filename << ">..." << endl;
+      
+      DataOut<dim> out;
+      ofstream o(filename.c_str());
+      fill_data (out);
+      out.add_data_vector (l1_error_per_dof, "L1-Error");
+      out.add_data_vector (l2_error_per_dof, "L2-Error");
+      out.add_data_vector (linfty_error_per_dof, "Linfty-Error");
+      out.add_data_vector (h1_seminorm_error_per_dof, "H1-seminorm-Error");
+      out.add_data_vector (h1_error_per_dof, "H1-Error");
+      out.write_ucd (o);
+      o.close ();
+    }
+  else
+    cout << "    Not writing error as grid." << endl;
+  
+  cout << endl;
+
+  delete fe;
+  delete quadrature;
+  delete boundary_quadrature;
+  
+  return dof->n_dofs();
+};
+
+
+template <int dim>
+void PoissonProblem<dim>::print_history (string filename) const {
+  ofstream out(filename.c_str());
+  out << "# n_dofs    l1_error l2_error linfty_error h1_seminorm_error h1_error"
+      << endl;
+  for (unsigned int i=0; i<n_dofs.size(); ++i)
+    out << n_dofs[i]
+       << "    "
+       << l1_error[i] << "  "
+       << l2_error[i] << "  "
+       << linfty_error[i] << "  "
+       << h1_seminorm_error[i] << "  "
+       << h1_error[i] << endl;
+
+  double average_l1=0,
+        average_l2=0,
+     average_linfty=0,
+    average_h1_semi=0,
+        average_h1=0;
+  for (unsigned int i=1; i<n_dofs.size(); ++i) 
+    {
+      average_l1 += l1_error[i]/l1_error[i-1];
+      average_l2 += l2_error[i]/l2_error[i-1];
+      average_linfty += linfty_error[i]/linfty_error[i-1];
+      average_h1_semi += h1_seminorm_error[i]/h1_seminorm_error[i-1];
+      average_h1 += h1_error[i]/h1_error[i-1];
+    };
+
+  average_l1 /= (l1_error.size()-1);
+  average_l2 /= (l1_error.size()-1);
+  average_linfty /= (l1_error.size()-1);
+  average_h1_semi /= (l1_error.size()-1);
+  average_h1 /= (l1_error.size()-1);
+
+  cout << "==========================================================\n";
+  cout << "Average error reduction rates for h->h/2:" << endl;
+  cout << "    L1 error         : " << 1./average_l1 << endl
+       << "    L2 error         : " << 1./average_l2 << endl
+       << "    Linfty error     : " << 1./average_linfty << endl
+       << "    H1 seminorm error: " << 1./average_h1_semi << endl
+       << "    H1 error         : " << 1./average_h1 << endl;
+  cout << "==========================================================\n";
+  cout << "==========================================================\n";
+};
+
+
+
+
+int main () {
+  for (unsigned int order=0; order<5; ++order) 
+    {
+      PoissonProblem<2> problem (order);
+      
+      unsigned int level=0;
+      unsigned int n_dofs;
+      do
+       n_dofs = problem.run (level++);
+      while (n_dofs<25000);
+
+      string filename;
+      switch (order) 
+       {
+         case 0:
+               filename = "criss_cross";
+               break;
+         case 1:
+               filename = "linear";
+               break;
+         case 2:
+               filename = "quadratic";
+               break;
+         case 3:
+               filename = "cubic";
+               break;
+         case 4:
+               filename = "quartic";
+               break;
+       };
+      filename += ".history";
+      
+      cout << endl << "Printing convergence history to <"
+          << filename << ">..." << endl;
+      problem.print_history (filename);
+      cout << endl << endl << endl;
+    };
+  
+  return 0;
+};
diff --git a/tests/big-tests/multigrid/make_ps b/tests/big-tests/multigrid/make_ps
new file mode 100644 (file)
index 0000000..c2df845
--- /dev/null
@@ -0,0 +1,52 @@
+set term postscript eps
+set xlabel "Number of degrees of freedom"
+set data style linespoints
+set logscale xy
+
+
+
+set ylabel "Error"
+
+set output "criss-cross.eps"
+
+plot "criss-cross.history" using 1:2 title "L1 error","criss-cross.history" using 1:3 title "L2 error","criss-cross.history" using 1:4 title "Linfty error","criss-cross.history" using 1:5 title "H1 seminorm error","criss-cross.history" using 1:6 title "H1 error"
+
+
+
+set output "linear.eps"
+
+plot "linear.history" using 1:2 title "L1 error","linear.history" using 1:3 title "L2 error","linear.history" using 1:4 title "Linfty error","linear.history" using 1:5 title "H1 seminorm error","linear.history" using 1:6 title "H1 error"
+
+
+
+set output "quadratic.eps"
+
+plot "quadratic.history" using 1:2 title "L1 error","quadratic.history" using 1:3 title "L2 error","quadratic.history" using 1:4 title "Linfty error","quadratic.history" using 1:5 title "H1 seminorm error","quadratic.history" using 1:6 title "H1 error"
+
+
+
+set output "cubic.eps"
+
+plot "cubic.history" using 1:2 title "L1 error","cubic.history" using 1:3 title "L2 error","cubic.history" using 1:4 title "Linfty error","cubic.history" using 1:5 title "H1 seminorm error","cubic.history" using 1:6 title "H1 error"
+
+
+
+set output "quartic.eps"
+
+plot "quartic.history" using 1:2 title "L1 error","quartic.history" using 1:3 title "L2 error","quartic.history" using 1:4 title "Linfty error","quartic.history" using 1:5 title "H1 seminorm error","quartic.history" using 1:6 title "H1 error"
+
+
+
+set output "l2error.eps"
+set ylabel "L2-error"
+
+plot "criss-cross.history" using 1:3 title "Criss-cross elements", "linear.history" using 1:3 title "Linear elements", "quadratic.history" using 1:3 title "Quadratic elements", "cubic.history" using 1:3 title "Cubic elements", "quartic.history" using 1:3 title "Quartic elements"
+
+
+
+set output "h1error.eps"
+set ylabel "H1-error"
+
+plot "criss-cross.history" using 1:6 title "Criss-cross elements", "linear.history" using 1:6 title "Linear elements", "quadratic.history" using 1:6 title "Quadratic elements", "cubic.history" using 1:6 title "Cubic elements", "quartic.history" using 1:6 title "Quartic elements"
+
+

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.