]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Bring up to standard with regard to multigrid on adaptive meshes.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 5 Feb 2010 04:07:25 +0000 (04:07 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 5 Feb 2010 04:07:25 +0000 (04:07 +0000)
git-svn-id: https://svn.dealii.org/trunk@20501 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-16/doc/intro.dox
deal.II/examples/step-16/step-16.cc

index 83f0c68a6538b6e7c52be48c0918d9397bd2bd1f..e86be0b0b6b57c81ffa69e57bd3edcdceaa146af 100644 (file)
@@ -1,3 +1,12 @@
+<br>
+
+<i>This program has evolved from a version originally written by Guido
+Kanschat in 2003. It has undergone significant revisions by B&auml;rbel
+Janssen, Guido Kanschat and Wolfgang Bangerth in 2009 and 2010 to demonstrate
+multigrid algorithms on adaptively refined meshes.
+</i>
+
+
 <a name="Intro"></a>
 <h1>Introduction</h1>
 
index 61f7d1f18da7ea843c1772b301482e0086f81531..2b562d310e35673e9a2a43ba3c35c726539cf1ff 100644 (file)
@@ -1,9 +1,11 @@
 /* $Id$ */
-/* Author: Guido Kanschat, University of Heidelberg, 2003 */
+/* Author: Guido Kanschat, University of Heidelberg, 2003  */
+/*         Baerbel Janssen, University of Heidelberg, 2010 */
+/*         Wolfgang Bangerth, Texas A&M University, 2010   */
 
 /*    $Id$       */
 /*                                                                */
-/*    Copyright (C) 2003, 2004, 2005, 2006, 2007, 2008, 2009 by the deal.II authors */
+/*    Copyright (C) 2003, 2004, 2006, 2007, 2008, 2009, 2010 by the deal.II authors                   */
 /*                                                                */
 /*    This file is subject to QPL and may not be  distributed     */
 /*    without copyright and license information. Please refer     */
 #include <base/quadrature_lib.h>
 #include <base/function.h>
 #include <base/logstream.h>
+#include <base/utilities.h>
+
+#include <lac/constraint_matrix.h>
 #include <lac/vector.h>
 #include <lac/full_matrix.h>
 #include <lac/sparse_matrix.h>
 #include <lac/solver_cg.h>
 #include <lac/precondition.h>
+
 #include <grid/tria.h>
 #include <grid/tria_accessor.h>
 #include <grid/tria_iterator.h>
-#include <grid/tria_boundary_lib.h>
 #include <grid/grid_generator.h>
-#include <dofs/dof_handler.h>
+
 #include <dofs/dof_accessor.h>
 #include <dofs/dof_tools.h>
+
 #include <fe/fe_q.h>
 #include <fe/fe_values.h>
+
 #include <numerics/vectors.h>
-#include <numerics/matrices.h>
 #include <numerics/data_out.h>
-                                // These are the new include files
-                                // required for multi-level methods.
-                                // First, the file defining the
-                                // multigrid method itself.
+
+//These are the same include files 
+//as in step-16 necessary for the
+//multi-level methods
 #include <multigrid/multigrid.h>
-                                // The DoFHandler is replaced by an
-                                // MGDoFHandler which is defined
-                                // here.
 #include <multigrid/mg_dof_handler.h>
 #include <multigrid/mg_dof_accessor.h>
-
-                                // Then, we need some pre-made
-                                // transfer routines between grids.
 #include <multigrid/mg_transfer.h>
-                                 // And a file in which equivalents to the
-                                 // DoFTools class are declared:
 #include <multigrid/mg_tools.h>
-
 #include <multigrid/mg_coarse.h>
 #include <multigrid/mg_smoother.h>
 #include <multigrid/mg_matrix.h>
                                 // previous programs:
 using namespace dealii;
 
-                                // This class is based on the same
-                                // class in step-5. Remark that we
-                                // replaced the DoFHandler by
-                                // MGDoFHandler. since this inherits
-                                // from DoFHandler, the new object
-                                // incorporates the old functionality
-                                // plus the new functions for degrees
-                                // of freedom on different
-                                // levels. Furthermore, we added
-                                // MGLevelObjects for sparsity
-                                // patterns and matrices.
+
+//This class is basically the same 
+//class as in step-16. The only 
+//difference is that here we solve Laplace's 
+//problem on an adaptively refined grid.
 template <int dim>
 class LaplaceProblem 
 {
   public:
-    LaplaceProblem ();
+    LaplaceProblem (const unsigned int deg);
     void run ();
-    
+
   private:
     void setup_system ();
     void assemble_system ();
-                                    // We add this function for
-                                    // assembling the multilevel
-                                    // matrices.
     void assemble_multigrid ();
     void solve ();
+    void refine_local ();
+    void test_get_coarse_cell ();
     void output_results (const unsigned int cycle) const;
 
     Triangulation<dim>   triangulation;
@@ -99,32 +89,34 @@ class LaplaceProblem
     SparsityPattern      sparsity_pattern;
     SparseMatrix<double> system_matrix;
 
-                                    // Here are the new objects for
-                                    // handling level matrices: sparsity
-                                    // patterns and matrices. We use number
-                                    // type float to save memory. It's only
-                                    // a preconditioner!
-    MGLevelObject<SparsityPattern>      mg_sparsity;
-    MGLevelObject<SparseMatrix<float> > mg_matrices;
-    
+                                    //This object holds the information f
+                                    //or the hanging nodes. 
+    ConstraintMatrix     constraints;
+
+    MGLevelObject<SparsityPattern> mg_sparsity;
+    MGLevelObject<SparseMatrix<double> > mg_matrices;
+
+                                    /* The matrices at the interface 
+                                     * between two refinement levels, 
+                                     * coupling coarse to fine.*/
+    MGLevelObject<SparseMatrix<double> > mg_interface_matrices_up;
+
     Vector<double>       solution;
     Vector<double>       system_rhs;
-};
 
+    const unsigned int degree;
+};
 
 
-                                // This function is as before.
 template <int dim>
-LaplaceProblem<dim>::LaplaceProblem () :
-                fe (1),
-               mg_dof_handler (triangulation)
+LaplaceProblem<dim>::LaplaceProblem (const unsigned int deg) :
+               triangulation (Triangulation<dim>::limit_level_difference_at_vertices),
+               fe (deg),
+               mg_dof_handler (triangulation),
+               degree(deg)
 {}
 
 
-
-                                // This is the function of step-5
-                                // augmented by the setup of the
-                                // multi-level matrices in the end.
 template <int dim>
 void LaplaceProblem<dim>::setup_system ()
 {
@@ -136,23 +128,33 @@ void LaplaceProblem<dim>::setup_system ()
                                   // multilevel structure
   deallog << "Number of degrees of freedom: "
          << mg_dof_handler.n_dofs();
+
   for (unsigned int l=0;l<triangulation.n_levels();++l)
     deallog << "   " << 'L' << l << ": "
            << mg_dof_handler.n_dofs(l);
   deallog  << std::endl;
-  
+
   sparsity_pattern.reinit (mg_dof_handler.n_dofs(),
                           mg_dof_handler.n_dofs(),
                           mg_dof_handler.max_couplings_between_dofs());
-  DoFTools::make_sparsity_pattern (static_cast<const DoFHandler<dim>&>(mg_dof_handler),
-                                   sparsity_pattern);
-  sparsity_pattern.compress();
-
-  system_matrix.reinit (sparsity_pattern);
+  DoFTools::make_sparsity_pattern (
+    static_cast<const DoFHandler<dim>&>(mg_dof_handler),
+    sparsity_pattern);
 
   solution.reinit (mg_dof_handler.n_dofs());
   system_rhs.reinit (mg_dof_handler.n_dofs());
 
+  constraints.clear ();
+  DoFTools::make_hanging_node_constraints (mg_dof_handler, constraints);
+  VectorTools::interpolate_boundary_values (mg_dof_handler,
+                                           0,
+                                           ZeroFunction<dim>(),
+                                           constraints);
+  constraints.close ();
+  constraints.condense (sparsity_pattern);
+  sparsity_pattern.compress();
+  system_matrix.reinit (sparsity_pattern);
+
                                   // The multi-level objects are
                                   // resized to hold matrices for
                                   // every level. The coarse level is
@@ -167,7 +169,11 @@ void LaplaceProblem<dim>::setup_system ()
                                   // SparsityPattern before it can be
                                   // destroyed.
   const unsigned int nlevels = triangulation.n_levels();
+
+  mg_interface_matrices_up.resize(0, nlevels-1);
+  mg_interface_matrices_up.clear ();
   mg_matrices.resize(0, nlevels-1);
+  mg_matrices.clear ();
   mg_sparsity.resize(0, nlevels-1);
 
                                   // Now, we have to build a matrix
@@ -183,8 +189,23 @@ void LaplaceProblem<dim>::setup_system ()
                                 mg_dof_handler.n_dofs(level),
                                 mg_dof_handler.max_couplings_between_dofs());
       MGTools::make_sparsity_pattern (mg_dof_handler, mg_sparsity[level], level);
+      CompressedSparsityPattern ci_sparsity;
+      if(level>0)
+       {
+         ci_sparsity.reinit(mg_dof_handler.n_dofs(level),
+                            mg_dof_handler.n_dofs(level));
+         MGTools::make_sparsity_pattern(mg_dof_handler, ci_sparsity, level);
+       }
+    }
+
+//And the same for the mg matrices 
+//for the interface. Note that there 
+//is no such interface on the coarsest level
+  for(unsigned int level=0; level<nlevels; ++level)
+    {
       mg_sparsity[level].compress();
       mg_matrices[level].reinit(mg_sparsity[level]);
+      mg_interface_matrices_up[level].reinit(mg_sparsity[level]);
     }
 }
 
@@ -194,11 +215,11 @@ void LaplaceProblem<dim>::setup_system ()
 template <int dim>
 void LaplaceProblem<dim>::assemble_system () 
 {  
-  QGauss<dim>  quadrature_formula(2);
+  QGauss<dim>  quadrature_formula(1+degree);
 
   FEValues<dim> fe_values (fe, quadrature_formula, 
                           update_values   | update_gradients |
-                           update_quadrature_points | update_JxW_values);
+                          update_quadrature_points | update_JxW_values);
 
   const unsigned int   dofs_per_cell = fe.dofs_per_cell;
   const unsigned int   n_q_points    = quadrature_formula.size();
@@ -208,27 +229,15 @@ void LaplaceProblem<dim>::assemble_system ()
 
   std::vector<unsigned int> local_dof_indices (dofs_per_cell);
 
-  typename DoFHandler<dim>::active_cell_iterator cell = mg_dof_handler.begin_active(),
-                                                endc = mg_dof_handler.end();
+  typename DoFHandler<dim>::active_cell_iterator
+    cell = mg_dof_handler.begin_active(),
+    endc = mg_dof_handler.end();
   for (; cell!=endc; ++cell)
     {
       cell_matrix = 0;
       cell_rhs = 0;
 
-                                      // As before, we want the FEValues
-                                      // object to compute the quantities
-                                      // which we told him to compute in
-                                      // the constructor using the update
-                                      // flags. Then, we loop over all
-                                      // quadrature points and the local
-                                      // matrix rows and columns for
-                                      // computing the element
-                                      // contribution. This is the same as
-                                      // in step-4. For the right hand
-                                      // side, we use a constant value of
-                                      // 1.
       fe_values.reinit (cell);
-
       for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
        for (unsigned int i=0; i<dofs_per_cell; ++i)
          {
@@ -239,34 +248,13 @@ void LaplaceProblem<dim>::assemble_system ()
 
            cell_rhs(i) += (fe_values.shape_value(i,q_point)
                            * 1.0 * fe_values.JxW(q_point));
-         };
+         }
 
       cell->get_dof_indices (local_dof_indices);
-      for (unsigned int i=0; i<dofs_per_cell; ++i)
-       {
-         for (unsigned int j=0; j<dofs_per_cell; ++j)
-           system_matrix.add (local_dof_indices[i],
-                              local_dof_indices[j],
-                              cell_matrix(i,j));
-         
-         system_rhs(local_dof_indices[i]) += cell_rhs(i);
-       };
-    };
-
-                                  // The Dirichlet boundary
-                                  // conditions on the finest level
-                                  // are handled as usual.
-  std::map<unsigned int,double> boundary_values;  
-
-  VectorTools::interpolate_boundary_values (mg_dof_handler,
-                                           0,
-                                           ZeroFunction<dim>(),
-                                           boundary_values);
-
-  MatrixTools::apply_boundary_values (boundary_values,
-                                     system_matrix,
-                                     solution,
-                                     system_rhs);
+      constraints.distribute_local_to_global (cell_matrix, cell_rhs,
+                                             local_dof_indices,
+                                             system_matrix, system_rhs);
+    }
 }
 
 
@@ -275,7 +263,10 @@ void LaplaceProblem<dim>::assemble_system ()
                                 // the same as above. Only the loop
                                 // goes over all existing cells now
                                 // and the results must be entered
-                                // into the correct matrix.
+                                // into the correct matrix. Here comes 
+                                 // the difference to global refinement 
+                                 // into play. We have to fill the interface 
+                                 // matrices correctly.
 
                                 // Since we only do multi-level
                                 // preconditioning, no right-hand
@@ -283,26 +274,57 @@ void LaplaceProblem<dim>::assemble_system ()
 template <int dim>
 void LaplaceProblem<dim>::assemble_multigrid () 
 {  
-  QGauss<dim>  quadrature_formula(2);
+  QGauss<dim>  quadrature_formula(1+degree);
 
   FEValues<dim> fe_values (fe, quadrature_formula, 
                           update_values   | update_gradients |
-                           update_quadrature_points | update_JxW_values);
+                          update_quadrature_points | update_JxW_values);
 
-  const unsigned int   dofs_per_cell = fe.dofs_per_cell;
-  const unsigned int   n_q_points    = quadrature_formula.size();
+  const unsigned int   dofs_per_cell   = fe.dofs_per_cell;
+  const unsigned int   n_q_points      = quadrature_formula.size();
 
   FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
 
   std::vector<unsigned int> local_dof_indices (dofs_per_cell);
 
+  std::vector<std::vector<bool> > interface_dofs;
+  std::vector<std::vector<bool> > boundary_interface_dofs;
+  for (unsigned int level = 0; level<triangulation.n_levels(); ++level)
+    {
+      std::vector<bool> tmp (mg_dof_handler.n_dofs(level));
+      interface_dofs.push_back (tmp);
+      boundary_interface_dofs.push_back (tmp);
+    }
+  MGTools::extract_inner_interface_dofs (mg_dof_handler,
+                                        interface_dofs,
+                                        boundary_interface_dofs);
+
+  typename FunctionMap<dim>::type      dirichlet_boundary;
+  ZeroFunction<dim>                    homogeneous_dirichlet_bc (1);
+  dirichlet_boundary[0] = &homogeneous_dirichlet_bc;
+
+  std::vector<std::set<unsigned int> > boundary_indices(triangulation.n_levels());
+  MGTools::make_boundary_list (mg_dof_handler, dirichlet_boundary,
+                              boundary_indices);
+
+  std::vector<ConstraintMatrix> boundary_constraints (triangulation.n_levels());
+  std::vector<ConstraintMatrix> boundary_interface_constraints (triangulation.n_levels());
+  for (unsigned int level=0; level<triangulation.n_levels(); ++level)
+    {
+      boundary_constraints[level].add_lines (interface_dofs[level]);
+      boundary_constraints[level].add_lines (boundary_indices[level]);
+      boundary_constraints[level].close ();
+
+      boundary_interface_constraints[level]
+       .add_lines (boundary_interface_dofs[level]);
+      boundary_interface_constraints[level].close ();
+    }
+  
   typename MGDoFHandler<dim>::cell_iterator cell = mg_dof_handler.begin(),
                                            endc = mg_dof_handler.end();
+
   for (; cell!=endc; ++cell)
     {
-                                      // Remember the level of the
-                                      // current cell.
-      const unsigned int level = cell->level();
       cell_matrix = 0;
 
                                       // Compute the values specified
@@ -321,7 +343,6 @@ void LaplaceProblem<dim>::assemble_multigrid ()
                                  * fe_values.JxW(q_point);
          }
 
-
                                       // Oops! This is a tiny
                                       // difference easily
                                       // forgotten. The indices we
@@ -332,66 +353,26 @@ void LaplaceProblem<dim>::assemble_multigrid ()
                                       // 'mg' entered into the
                                       // function call.
       cell->get_mg_dof_indices (local_dof_indices);
+
+      const unsigned int level = cell->level();
+      boundary_constraints[level]
+       .distribute_local_to_global (cell_matrix,
+                                    local_dof_indices,
+                                    mg_matrices[level]);
+    
       for (unsigned int i=0; i<dofs_per_cell; ++i)
-       {
-         for (unsigned int j=0; j<dofs_per_cell; ++j)
-                                            // And now add everything
-                                            // to the matrix on the
-                                            // right level.
-           mg_matrices[level].add (local_dof_indices[i],
-                                    local_dof_indices[j],
-                                    cell_matrix(i,j));
-       }
+       for (unsigned int j=0; j<dofs_per_cell; ++j)
+         if( !(interface_dofs[level][local_dof_indices[i]]==true && 
+               interface_dofs[level][local_dof_indices[j]]==false))
+           cell_matrix(i,j) = 0;
+
+      boundary_interface_constraints[level]
+       .distribute_local_to_global (cell_matrix,
+                                    local_dof_indices,
+                                    mg_interface_matrices_up[level]);
     }
-
-                                  // Now we have to eliminate the
-                                  // boundary nodes on all
-                                  // levels. This is done exactly in
-                                  // the same fashion as on the
-                                  // finest level. Therefore, we need
-                                  // a function map first. On the
-                                  // other hand, since we use
-                                  // multigrid only as a
-                                  // preconditioner, we always use
-                                  // homogeneous boundary conditions
-                                  // and the ZeroFunction will be
-                                  // sufficient in all cases; using a
-                                  // different function here would
-                                  // not hurt on the other hand,
-                                  // since the values are ignored.
-
-  typename FunctionMap<dim>::type      dirichlet_boundary;
-  ZeroFunction<dim>                    homogeneous_dirichlet_bc (1);
-  dirichlet_boundary[0] = &homogeneous_dirichlet_bc;
-  
-                                  // Next we generate the set of
-                                  // dof indices to be eliminated on
-                                  // each level.
-  std::vector<std::set<unsigned int> > boundary_indices(triangulation.n_levels());
-  MGTools::make_boundary_list (mg_dof_handler,
-                              dirichlet_boundary,
-                              boundary_indices);
-
-                                  // And finally we eliminate these
-                                  // degrees of freedom from every
-                                  // matrix in the multilevel hierarchy.
- const unsigned int nlevels = triangulation.n_levels();
- for (unsigned int level=0;level<nlevels;++level)
-   {
-     MGTools::apply_boundary_values(boundary_indices[level],
-                                   mg_matrices[level],
-                                   true);
-   }
 }
 
-
-
-                                // The solution process again looks
-                                // mostly like in the previous
-                                // examples. However, we will now use
-                                // a preconditioned conjugate
-                                // gradient algorithm. It is not very
-                                // difficult to make this change:
 template <int dim>
 void LaplaceProblem<dim>::solve () 
 {
@@ -406,7 +387,7 @@ void LaplaceProblem<dim>::solve ()
                                   // the transfer of functions
                                   // between different grid
                                   // levels.
-  MGTransferPrebuilt<Vector<double> > mg_transfer;
+  MGTransferPrebuilt<Vector<double> > mg_transfer(constraints);
   mg_transfer.build_matrices(mg_dof_handler);
 
                                   // Next, we need a coarse grid
@@ -415,9 +396,9 @@ void LaplaceProblem<dim>::solve ()
                                   // direct solver, even if its
                                   // implementation is not very
                                   // clever.
-  FullMatrix<float> coarse_matrix;
+  FullMatrix<double> coarse_matrix;
   coarse_matrix.copy_from (mg_matrices[0]);
-  MGCoarseGridHouseholder<float, Vector<double> > mg_coarse;
+  MGCoarseGridHouseholder<double, Vector<double> > mg_coarse;
   mg_coarse.initialize(coarse_matrix);
 
                                   // The final ingredient for the
@@ -427,20 +408,22 @@ void LaplaceProblem<dim>::solve ()
                                   // here. Names are getting quite
                                   // long here, so we help with
                                   // typedefs.
-  typedef PreconditionSOR<SparseMatrix<float> > RELAXATION;
-  MGSmootherRelaxation<SparseMatrix<float>, RELAXATION, Vector<double> >
+  typedef PreconditionSOR<SparseMatrix<double> > RELAXATION;
+//  typedef PreconditionJacobi<SparseMatrix<double> > RELAXATION;
+//  typedef SparseILU<double> RELAXATION;
+  MGSmootherRelaxation<SparseMatrix<double>, RELAXATION, Vector<double> >
     mg_smoother(vector_memory);
 
                                   // Initialize the smoother with our
-                                  // level matrices and the (required)
+                                  // level matrices and the required,
                                   // additional data for the
-                                  // relaxation method with default
+                                  // relaxaton method with default
                                   // values.
-  RELAXATION::AdditionalData smoother_data;
+  RELAXATION::AdditionalData smoother_data;//(0, 9,false);
   mg_smoother.initialize(mg_matrices, smoother_data);
-  
+
                                   // Do two smoothing steps per level
-  mg_smoother.set_steps(2);
+  mg_smoother.set_steps(1);
                                   // Since the SOR method is not
                                   // symmetric, but we use conjugate
                                   // gradient iteration below, here
@@ -449,12 +432,19 @@ void LaplaceProblem<dim>::solve ()
                                   // symmetric operator even for
                                   // nonsymmetric smoothers.
   mg_smoother.set_symmetric(true);
+  mg_smoother.set_variable(true);
+
 
                                   // We must wrap our matrices in an
                                   // object having the required
                                   // multiplication functions.
-  MGMatrix<SparseMatrix<float>, Vector<double> >
+  MGMatrix<SparseMatrix<double>, Vector<double> >
     mg_matrix(&mg_matrices);
+                                  //do the same for the interface matrices
+  MGMatrix<SparseMatrix<double>, Vector<double> >
+    mg_interface_up(&mg_interface_matrices_up);
+  MGMatrix<SparseMatrix<double>, Vector<double> >
+    mg_interface_down(&mg_interface_matrices_up);
                                   // Now, we are ready to set up the
                                   // V-cycle operator and the
                                   // multilevel preconditioner.
@@ -464,66 +454,177 @@ void LaplaceProblem<dim>::solve ()
                                mg_transfer,
                                mg_smoother,
                                mg_smoother);
+  mg.set_edge_matrices(mg_interface_down, mg_interface_up);
+
   PreconditionMG<dim, Vector<double>,
     MGTransferPrebuilt<Vector<double> > >
-    preconditioner(mg_dof_handler, mg, mg_transfer);
-  
-                                  // Finally, create the solver
-                                  // object and solve the system
-  SolverControl           solver_control (1000, 1e-12);
-  SolverCG<>              cg (solver_control);
+  preconditioner(mg_dof_handler, mg, mg_transfer);
 
-  
-  cg.solve (system_matrix, solution, system_rhs,
-           preconditioner);
+                                // Finally, create the solver
+                                // object and solve the system
+ReductionControl solver_control (100, 1.e-20, 1.e-10, true, true);
+SolverCG<>    cg (solver_control);
+
+solution = 0;
+
+cg.solve (system_matrix, solution, system_rhs,
+         preconditioner);
+constraints.distribute (solution);
+
+std::cout << "   " << solver_control.last_step()
+<< " CG iterations needed to obtain convergence."
+<< std::endl;
 }
 
+template <int dim>
+void LaplaceProblem<dim>::refine_local ()
+{
+  bool cell_refined = false;
+  for (typename Triangulation<dim>::active_cell_iterator
+        cell = triangulation.begin_active();
+       cell != triangulation.end(); ++cell)
+    {
+      if(false)
+       {
+         for (unsigned int vertex=0;
+              vertex < GeometryInfo<dim>::vertices_per_cell;
+              ++vertex)
+           {
+             if(cell->vertex(vertex)[dim-1]==0 && cell->vertex(vertex)[0]==0 && cell->vertex(vertex)[dim-2]==0)
+               {
+                 cell->set_refine_flag ();
+                 cell_refined = true;
+                 break;
+               }
+           }
+       }
+      else if(true) //Kreis
+       {
+         for (unsigned int vertex=0;
+              vertex < GeometryInfo<dim>::vertices_per_cell;
+              ++vertex)
+           {
+             const Point<dim> p = cell->vertex(vertex);
+             const Point<dim> origin = (dim == 2 ?
+                                        Point<dim>(0,0) :
+                                        Point<dim>(0,0,0));
+             const double dist = p.distance(origin);
+             if(dist<0.25/M_PI)
+               {
+                 cell->set_refine_flag ();
+                 cell_refined = true;
+                 break;
+               }
+           }
+       }
+      else if(false) //linke Diagonale
+       {
+         for (unsigned int vertex=0;
+              vertex < GeometryInfo<dim>::vertices_per_cell;
+              ++vertex)
+           {
+             const Point<dim> p = cell->vertex(vertex);
+             if(p[0]==p[1])
+               {
+                 cell->set_refine_flag ();
+                 cell_refined = true;
+                 break;
+               }
+           }
+       }
+      else if(false) //inneres Quadrat
+       {
+         for (unsigned int vertex=0;
+              vertex < GeometryInfo<dim>::vertices_per_cell;
+              ++vertex)
+           {
+             const Point<dim> p = cell->vertex(vertex);
+             const double dist = std::max(std::fabs(p[0]),std::fabs(p[1]));
+             if(dist<0.5)
+               {
+                 cell->set_refine_flag ();
+                 cell_refined = true;
+                 break;
+               }
+           }
+       }
+      else if(false) //Raute
+       {
+         for (unsigned int vertex=0;
+              vertex < GeometryInfo<dim>::vertices_per_cell;
+              ++vertex)
+           {
+             const Point<dim> p = cell->vertex(vertex);
+             const double dist = std::fabs(p[0])+std::fabs(p[1]);
+             if(dist<0.5)
+               {
+                 cell->set_refine_flag ();
+                 cell_refined = true;
+                 break;
+               }
+           }
+       }
+      else //erster Quadrant 
+       {
+         const Point<dim> p = cell->center();
+         bool positive = p(0) > 0;
+         if (dim>1 && p(1) <= 0)
+           positive = false;
+         if (dim>2 && p(2) <= 0)
+           positive = false;
+         if (positive)
+           {
+             cell->set_refine_flag();
+             cell_refined = true;
+           }
+       }
+    }
+                                  //Wenn nichts verfeinert wurde bisher, global verfeinern!
+  if(!cell_refined)
+    for (typename Triangulation<dim>::active_cell_iterator
+          cell = triangulation.begin_active();
+        cell != triangulation.end(); ++cell)
+      cell->set_refine_flag();
+
 
+  triangulation.execute_coarsening_and_refinement ();
+}
 
-                                // Here is the data output, which is
-                                // a simplified version of step-5. We
-                                // do a standard gnuplot output for
-                                // each grid produced in the
-                                // refinement process.
 template <int dim>
 void LaplaceProblem<dim>::output_results (const unsigned int cycle) const
 {
-                                  // Construct and initialize a DataOut object
   DataOut<dim> data_out;
 
   data_out.attach_dof_handler (mg_dof_handler);
   data_out.add_data_vector (solution, "solution");
   data_out.build_patches ();
 
-                                  // The following block generates
-                                  // the file name and opens the
-                                  // file:
   std::ostringstream filename;
   filename << "solution-"
           << cycle
-          << ".gnuplot";
+          << ".vtk";
 
   std::ofstream output (filename.str().c_str());
-  data_out.write_gnuplot (output);
+  data_out.write_vtk (output);
 }
 
-
-
 template <int dim>
 void LaplaceProblem<dim>::run () 
 {
-  for (unsigned int cycle=0; cycle<6; ++cycle)
+  for (unsigned int cycle=0; cycle<9; ++cycle)
     {
-      deallog << "Cycle " << cycle << std::endl;
-
       if (cycle == 0)
        {
-                                          // Generate a simple hyperball grid.
-         GridGenerator::hyper_ball(triangulation);
-         static const HyperBallBoundary<dim> boundary;
-         triangulation.set_boundary (0, boundary);
+         GridGenerator::hyper_cube(triangulation, -1, 1);
+         triangulation.refine_global (1);
        }
-      triangulation.refine_global (1);
+      refine_local ();
+
+      std::cout << "Cycle " << cycle
+               << " with " << triangulation.n_active_cells()
+               << " cells."
+               << std::endl;
+
       setup_system ();
       assemble_system ();
       assemble_multigrid ();
@@ -532,16 +633,41 @@ void LaplaceProblem<dim>::run ()
     };
 }
 
-    
 
-                                // The main function looks mostly
-                                // like the one in the previous
-                                // example, so we won't comment on it
-                                // further.
+
 int main () 
 {
-  LaplaceProblem<2> laplace_problem_2d;
-  laplace_problem_2d.run ();
-  
+  try
+    {
+      deallog.depth_console (0);
+
+      LaplaceProblem<2> laplace_problem(1);
+      laplace_problem.run ();
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+               << exc.what() << std::endl
+               << "Aborting!" << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+
+      return 1;
+    }
+  catch (...) 
+    {
+      std::cerr << std::endl << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+               << "Aborting!" << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      return 1;
+    }
+
   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.