]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add step-6.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Sun, 2 Jan 2000 20:34:42 +0000 (20:34 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Sun, 2 Jan 2000 20:34:42 +0000 (20:34 +0000)
git-svn-id: https://svn.dealii.org/trunk@2141 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/Attic/examples/step-by-step/step-6/Makefile [new file with mode: 0644]
deal.II/deal.II/Attic/examples/step-by-step/step-6/step-6.cc [new file with mode: 0644]
deal.II/examples/step-6/Makefile [new file with mode: 0644]
deal.II/examples/step-6/step-6.cc [new file with mode: 0644]

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 (file)
index 0000000..dc2798f
--- /dev/null
@@ -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 (file)
index 0000000..b87329d
--- /dev/null
@@ -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 (file)
index 0000000..dc2798f
--- /dev/null
@@ -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 (file)
index 0000000..b87329d
--- /dev/null
@@ -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;
+};

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.