]> https://gitweb.dealii.org/ - dealii.git/commitdiff
New test
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 24 May 2007 01:28:41 +0000 (01:28 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 24 May 2007 01:28:41 +0000 (01:28 +0000)
git-svn-id: https://svn.dealii.org/trunk@14694 0785d39b-7218-0410-832d-ea1e28bc413d

tests/fail/hp-crash_18.cc [new file with mode: 0644]

diff --git a/tests/fail/hp-crash_18.cc b/tests/fail/hp-crash_18.cc
new file mode 100644 (file)
index 0000000..8282ce3
--- /dev/null
@@ -0,0 +1,735 @@
+//----------------------------  crash_18.cc  ---------------------------
+//    $Id: crash_18.cc 12732 2006-03-28 23:15:45Z wolf $
+//    Version: $Name$ 
+//
+//    Copyright (C) 2006, 2007 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  crash_18.cc  ---------------------------
+
+
+// a modified version of step-27 that crashes because of circular
+// constraints. like crash_17, but only check the second cycle. there,
+// we have dofs that are constrained against itself, and we don't
+// currently detect this...
+
+char logname[] = "hp-crash_18/output";
+
+
+#include "../tests.h"
+
+
+#include <base/quadrature_lib.h>
+#include <base/function.h>
+#include <base/logstream.h>
+#include <base/utilities.h>
+#include <base/timer.h>
+#include <lac/vector.h>
+#include <lac/full_matrix.h>
+#include <lac/sparse_matrix.h>
+#include <lac/compressed_sparsity_pattern.h>
+#include <lac/solver_cg.h>
+#include <lac/precondition.h>
+#include <grid/tria.h>
+#include <dofs/hp_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/hp_fe_values.h>
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <numerics/data_out.h>
+
+#include <fstream>
+#include <iostream>
+
+#include <fe/fe_q.h>
+#include <grid/grid_out.h>
+#include <dofs/dof_constraints.h>
+#include <grid/grid_refinement.h>
+#include <numerics/error_estimator.h>
+
+#include <complex>
+
+                                // Finally, this is as in previous
+                                // programs:
+using namespace dealii;
+
+template <int dim>
+class LaplaceProblem 
+{
+  public:
+    LaplaceProblem ();
+    ~LaplaceProblem ();
+
+    void run ();
+    
+  private:
+    void setup_system ();
+    void assemble_system ();
+    void solve ();
+    void create_coarse_grid ();
+    void refine_grid ();
+    void estimate_smoothness (Vector<float> &smoothness_indicators) const;
+    void output_results (const unsigned int cycle) const;
+
+    Triangulation<dim>   triangulation;
+
+    hp::DoFHandler<dim>      dof_handler;
+    hp::FECollection<dim>    fe_collection;
+    hp::QCollection<dim>     quadrature_collection;
+    hp::QCollection<dim-1>   face_quadrature_collection;
+
+    ConstraintMatrix     hanging_node_constraints;
+
+    SparsityPattern      sparsity_pattern;
+    SparseMatrix<double> system_matrix;
+
+    Vector<double>       solution;
+    Vector<double>       system_rhs;
+
+    Timer distr, condense, hang;
+};
+
+
+
+template <int dim>
+class RightHandSide : public Function<dim>
+{
+  public:
+    RightHandSide () : Function<dim> () {};
+    
+    virtual double value (const Point<dim>   &p,
+                         const unsigned int  component) const;
+};
+
+
+template <int dim>
+double
+RightHandSide<dim>::value (const Point<dim>   &p,
+                          const unsigned int  /*component*/) const
+{
+  double product = 1;
+  for (unsigned int d=0; d<dim; ++d)
+    product *= (p[d]+1);
+  return product;
+}
+
+
+
+
+template <int dim>
+LaplaceProblem<dim>::LaplaceProblem () :
+               dof_handler (triangulation)
+{
+  for (unsigned int degree=2; degree<5; ++degree)
+    {
+      fe_collection.push_back (FE_Q<dim>(degree));
+      quadrature_collection.push_back (QGauss<dim>(degree+2));
+      face_quadrature_collection.push_back (QGauss<dim-1>(degree+2));
+    }
+}
+
+
+template <int dim>
+LaplaceProblem<dim>::~LaplaceProblem () 
+{
+  dof_handler.clear ();
+}
+
+template <int dim>
+void LaplaceProblem<dim>::setup_system ()
+{
+  distr.reset();
+  distr.start();
+  dof_handler.distribute_dofs (fe_collection);
+  distr.stop();
+
+  solution.reinit (dof_handler.n_dofs());
+  system_rhs.reinit (dof_handler.n_dofs());
+
+  hanging_node_constraints.clear ();
+  hang.reset();
+  hang.start();
+  DoFTools::make_hanging_node_constraints (dof_handler,
+                                          hanging_node_constraints);
+
+                                  // check for self-referential constraints
+  
+  hanging_node_constraints.close ();
+  hang.stop();
+
+  if (dim < 3)
+    {
+      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);
+
+      condense.reset();
+      condense.start();
+      hanging_node_constraints.condense (sparsity_pattern);
+      condense.stop();
+      sparsity_pattern.compress();      
+    }
+  else
+    {
+      CompressedSparsityPattern csp (dof_handler.n_dofs(),
+                                    dof_handler.n_dofs());
+      DoFTools::make_sparsity_pattern (dof_handler, csp);
+
+      condense.reset();
+      condense.start();
+      hanging_node_constraints.condense (csp);
+      condense.stop();
+      
+      sparsity_pattern.copy_from (csp);
+    }  
+
+  system_matrix.reinit (sparsity_pattern);
+}
+
+template <int dim>
+void LaplaceProblem<dim>::assemble_system () 
+{
+  hp::FEValues<dim> hp_fe_values (fe_collection,
+                                 quadrature_collection, 
+                                 update_values    |  update_gradients |
+                                 update_q_points  |  update_JxW_values);
+
+  const RightHandSide<dim> rhs_function;
+  
+  typename hp::DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+  for (; cell!=endc; ++cell)
+    {
+      const unsigned int   dofs_per_cell = cell->get_fe().dofs_per_cell;
+      FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
+      Vector<double>       cell_rhs (dofs_per_cell);
+
+      std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+
+      cell_matrix = 0;
+      cell_rhs = 0;
+
+      hp_fe_values.reinit (cell);
+      
+      const FEValues<dim> &fe_values = hp_fe_values.get_present_fe_values ();
+
+      std::vector<double>  rhs_values (fe_values.n_quadrature_points);
+      rhs_function.value_list (fe_values.get_quadrature_points(),
+                              rhs_values);
+      
+      for (unsigned int q_point=0; q_point<fe_values.n_quadrature_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) += (fe_values.shape_grad(i,q_point) *
+                                  fe_values.shape_grad(j,q_point) *
+                                  fe_values.JxW(q_point));
+
+           cell_rhs(i) += (fe_values.shape_value(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);
+       }
+    }
+
+  condense.start();
+  hanging_node_constraints.condense (system_matrix);
+  hanging_node_constraints.condense (system_rhs);
+  condense.stop();
+  std::map<unsigned int,double> boundary_values;
+  VectorTools::interpolate_boundary_values (dof_handler,
+                                           0,
+                                           ZeroFunction<dim>(),
+                                           boundary_values);
+  MatrixTools::apply_boundary_values (boundary_values,
+                                     system_matrix,
+                                     solution,
+                                     system_rhs);
+}
+
+template <int dim>
+void LaplaceProblem<dim>::solve () 
+{
+  SolverControl           solver_control (system_rhs.size(),
+                                         1e-8*system_rhs.l2_norm());
+  SolverCG<>              cg (solver_control);
+
+  PreconditionSSOR<> preconditioner;
+  preconditioner.initialize(system_matrix, 1.2);
+
+  cg.solve (system_matrix, solution, system_rhs,
+           preconditioner);
+
+  hanging_node_constraints.distribute (solution);
+}
+
+
+unsigned int
+int_pow (const unsigned int x,
+        const unsigned int n)
+{
+  unsigned int p=1;
+  for (unsigned int i=0; i<n; ++i)
+    p *= x;
+  return p;
+}
+
+
+template <int dim>
+void
+LaplaceProblem<dim>::
+estimate_smoothness (Vector<float> &smoothness_indicators) const
+{
+  const unsigned int N = 4;
+
+                                  // form all the Fourier vectors
+                                  // that we want to
+                                  // consider. exclude k=0 to avoid
+                                  // problems with |k|^{-mu} and also
+                                  // logarithms of |k|
+  std::vector<Tensor<1,dim> > k_vectors;
+  std::vector<unsigned int>   k_vectors_magnitude;
+  switch (dim)
+    {
+      case 2:
+      {
+       for (unsigned int i=0; i<N; ++i)
+         for (unsigned int j=0; j<N; ++j)
+           if (!((i==0) && (j==0))
+               &&
+               (i*i + j*j < N*N))
+             {
+               k_vectors.push_back (Point<dim>(deal_II_numbers::PI * i,
+                                               deal_II_numbers::PI * j));
+               k_vectors_magnitude.push_back (i*i+j*j);
+             }
+       
+       break;
+      }
+
+      case 3:
+      {
+       for (unsigned int i=0; i<N; ++i)
+         for (unsigned int j=0; j<N; ++j)
+           for (unsigned int k=0; k<N; ++k)
+             if (!((i==0) && (j==0) && (k==0))
+                 &&
+                 (i*i + j*j + k*k < N*N))
+               {
+                 k_vectors.push_back (Point<dim>(deal_II_numbers::PI * i,
+                                                 deal_II_numbers::PI * j,
+                                                 deal_II_numbers::PI * k));
+                 k_vectors_magnitude.push_back (i*i+j*j+k*k);
+             }
+       
+       break;
+      }
+      
+      default:
+           Assert (false, ExcNotImplemented());
+    }
+
+  const unsigned n_fourier_modes = k_vectors.size();
+  std::vector<double> ln_k (n_fourier_modes);
+  for (unsigned int i=0; i<n_fourier_modes; ++i)
+    ln_k[i] = std::log (k_vectors[i].norm());
+  
+
+                                  // assemble the matrices that do
+                                  // the Fourier transforms for each
+                                  // of the finite elements we deal
+                                  // with. note that these matrices
+                                  // are complex-valued, so we can't
+                                  // use FullMatrix
+  QGauss<1>      base_quadrature (2);
+  QIterated<dim> quadrature (base_quadrature, N);
+  
+  std::vector<Table<2,std::complex<double> > >
+    fourier_transform_matrices (fe_collection.size());
+  for (unsigned int fe=0; fe<fe_collection.size(); ++fe)
+    {
+      fourier_transform_matrices[fe].reinit (n_fourier_modes,
+                                            fe_collection[fe].dofs_per_cell);
+
+      for (unsigned int k=0; k<n_fourier_modes; ++k)
+       for (unsigned int i=0; i<fe_collection[fe].dofs_per_cell; ++i)
+         {
+           std::complex<double> sum = 0;
+           for (unsigned int q=0; q<quadrature.n_quadrature_points; ++q)
+             {
+               const Point<dim> x_q = quadrature.point(q);
+               sum += std::exp(std::complex<double>(0,1) *
+                               (k_vectors[k] * x_q)) *
+                      fe_collection[fe].shape_value(i,x_q) *
+                      quadrature.weight(q);
+             }
+           fourier_transform_matrices[fe](k,i)
+             = sum / std::pow(2*deal_II_numbers::PI, 1.*dim/2);
+         }
+    }
+
+                                  // the next thing is to loop over
+                                  // all cells and do our work there,
+                                  // i.e. to locally do the Fourier
+                                  // transform and estimate the decay
+                                  // coefficient
+  std::vector<std::complex<double> > fourier_coefficients (n_fourier_modes);
+  Vector<double>                     local_dof_values;
+  
+  typename hp::DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+  for (unsigned int index=0; cell!=endc; ++cell, ++index)
+    {
+      local_dof_values.reinit (cell->get_fe().dofs_per_cell);
+      cell->get_dof_values (solution, local_dof_values);
+
+                                      // first compute the Fourier
+                                      // transform of the local
+                                      // solution
+      std::fill (fourier_coefficients.begin(), fourier_coefficients.end(), 0);
+      for (unsigned int f=0; f<n_fourier_modes; ++f)
+       for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
+         fourier_coefficients[f] += 
+           fourier_transform_matrices[cell->active_fe_index()](f,i)
+           *
+           local_dof_values(i);
+
+                                      // enter the Fourier
+                                      // coefficients into a map with
+                                      // the k-magnitudes, to make
+                                      // sure that we get only the
+                                      // largest magnitude for each
+                                      // value of |k|
+      std::map<unsigned int, double> k_to_max_U_map;
+      for (unsigned int f=0; f<n_fourier_modes; ++f)
+       if ((k_to_max_U_map.find (k_vectors_magnitude[f]) ==
+            k_to_max_U_map.end())
+           ||
+           (k_to_max_U_map[k_vectors_magnitude[f]] <
+            std::abs (fourier_coefficients[f])))
+         k_to_max_U_map[k_vectors_magnitude[f]]
+           = std::abs (fourier_coefficients[f]);
+      
+                                      // now we have to calculate the
+                                      // various contributions to the
+                                      // formula for mu. we'll only
+                                      // take those fourier
+                                      // coefficients with the
+                                      // largest value for a given
+                                      // |k|
+      double  sum_1           = 0,
+             sum_ln_k        = 0,
+             sum_ln_k_square = 0,
+             sum_ln_U        = 0,
+             sum_ln_U_ln_k   = 0;
+      for (unsigned int f=0; f<n_fourier_modes; ++f)
+       if (k_to_max_U_map[k_vectors_magnitude[f]] ==
+           std::abs (fourier_coefficients[f]))
+         {
+           sum_1 += 1;
+           sum_ln_k += ln_k[f];
+           sum_ln_k_square += ln_k[f]*ln_k[f];
+           sum_ln_U += std::log (std::abs (fourier_coefficients[f]));
+           sum_ln_U_ln_k += std::log (std::abs (fourier_coefficients[f])) * ln_k[f];
+         }
+
+      const double mu
+       = (1./(sum_1*sum_ln_k_square - sum_ln_k*sum_ln_k)
+          *
+          (sum_ln_k*sum_ln_U - sum_1*sum_ln_U_ln_k));
+      
+      smoothness_indicators(index) = mu - 1.*dim/2;
+    }
+}
+
+
+
+  
+template <int dim>
+void LaplaceProblem<dim>::refine_grid ()
+{
+  Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+  KellyErrorEstimator<dim>::estimate (dof_handler,
+                                     face_quadrature_collection,
+                                     typename FunctionMap<dim>::type(),
+                                     solution,
+                                     estimated_error_per_cell);
+
+  Vector<float> smoothness_indicators (triangulation.n_active_cells());
+  estimate_smoothness (smoothness_indicators);
+  
+  GridRefinement::refine_and_coarsen_fixed_number (triangulation,
+                                                  estimated_error_per_cell,
+                                                  0.3, 0.03);
+
+  float max_smoothness = 0,
+        min_smoothness = smoothness_indicators.linfty_norm();
+  {
+    typename hp::DoFHandler<dim>::active_cell_iterator
+      cell = dof_handler.begin_active(),
+      endc = dof_handler.end();
+    for (unsigned int index=0; cell!=endc; ++cell, ++index)
+      if (cell->refine_flag_set())
+       {
+         max_smoothness = std::max (max_smoothness,
+                                    smoothness_indicators(index));
+         min_smoothness = std::min (min_smoothness,
+                                    smoothness_indicators(index));
+       }
+  }
+  const float cutoff_smoothness = (max_smoothness + min_smoothness) / 2;
+  {
+    typename hp::DoFHandler<dim>::active_cell_iterator
+      cell = dof_handler.begin_active(),
+      endc = dof_handler.end();
+    for (unsigned int index=0; cell!=endc; ++cell, ++index)
+      if (cell->refine_flag_set()
+         &&
+         (smoothness_indicators(index) > cutoff_smoothness)
+         &&
+         !(cell->active_fe_index() == fe_collection.size() - 1))
+       {
+         cell->clear_refine_flag();
+         cell->set_active_fe_index (std::min (cell->active_fe_index() + 1,
+                                              fe_collection.size() - 1));
+       }
+  } 
+  
+  triangulation.execute_coarsening_and_refinement ();
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::output_results (const unsigned int cycle) const
+{
+  Assert (cycle < 10, ExcNotImplemented());
+  
+  {
+    const std::string filename = "grid-" +
+                                Utilities::int_to_string (cycle, 2) +
+                                ".eps";
+    std::ofstream output (filename.c_str());
+    
+    GridOut grid_out;
+    grid_out.write_eps (triangulation, output);
+  }
+  
+  {
+    Vector<float> smoothness_indicators (triangulation.n_active_cells());
+    estimate_smoothness (smoothness_indicators);
+
+    Vector<double> smoothness_field (dof_handler.n_dofs());
+    DoFTools::distribute_cell_to_dof_vector (dof_handler,
+                                            smoothness_indicators,
+                                            smoothness_field);
+
+    Vector<float> fe_indices (triangulation.n_active_cells());
+    {
+      typename hp::DoFHandler<dim>::active_cell_iterator
+       cell = dof_handler.begin_active(),
+       endc = dof_handler.end();
+      for (unsigned int index=0; cell!=endc; ++cell, ++index)
+       {
+         
+         fe_indices(index) = cell->active_fe_index();
+//       smoothness_indicators(index) *= std::sqrt(cell->diameter());
+       }
+    }
+    
+    const std::string filename = "solution-" +
+                                Utilities::int_to_string (cycle, 2) +
+                                ".vtk";
+    DataOut<dim,hp::DoFHandler<dim> > data_out;
+
+    data_out.attach_dof_handler (dof_handler);
+    data_out.add_data_vector (solution, "solution");
+    data_out.add_data_vector (smoothness_indicators, "smoothness1");
+    data_out.add_data_vector (smoothness_field, "smoothness2");
+    data_out.add_data_vector (fe_indices, "fe_index");
+    data_out.build_patches ();
+  
+    std::ofstream output (filename.c_str());
+    data_out.write_vtk (output);
+  }
+}
+
+
+template <>
+void LaplaceProblem<2>::create_coarse_grid ()
+{
+  const unsigned int dim = 2;
+  
+  static const Point<2> vertices_1[]
+    = {  Point<2> (-1.,   -1.),
+         Point<2> (-1./2, -1.),
+         Point<2> (0.,    -1.),
+         Point<2> (+1./2, -1.),
+         Point<2> (+1,    -1.),
+            
+         Point<2> (-1.,   -1./2.),
+         Point<2> (-1./2, -1./2.),
+         Point<2> (0.,    -1./2.),
+         Point<2> (+1./2, -1./2.),
+         Point<2> (+1,    -1./2.),
+            
+         Point<2> (-1.,   0.),
+         Point<2> (-1./2, 0.),
+         Point<2> (+1./2, 0.),
+         Point<2> (+1,    0.),
+            
+         Point<2> (-1.,   1./2.),
+         Point<2> (-1./2, 1./2.),
+         Point<2> (0.,    1./2.),
+         Point<2> (+1./2, 1./2.),
+         Point<2> (+1,    1./2.),
+            
+         Point<2> (-1.,   1.),
+         Point<2> (-1./2, 1.),
+         Point<2> (0.,    1.),                   
+         Point<2> (+1./2, 1.),
+         Point<2> (+1,    1.)    };
+  const unsigned int
+    n_vertices = sizeof(vertices_1) / sizeof(vertices_1[0]);
+  const std::vector<Point<dim> > vertices (&vertices_1[0],
+                                           &vertices_1[n_vertices]);
+  static const int cell_vertices[][GeometryInfo<dim>::vertices_per_cell]
+    = {{0, 1, 5, 6},
+       {1, 2, 6, 7},
+       {2, 3, 7, 8},
+       {3, 4, 8, 9},
+       {5, 6, 10, 11},
+       {8, 9, 12, 13},
+       {10, 11, 14, 15},
+       {12, 13, 17, 18},
+       {14, 15, 19, 20},
+       {15, 16, 20, 21},
+       {16, 17, 21, 22},
+       {17, 18, 22, 23}};
+  const unsigned int
+    n_cells = sizeof(cell_vertices) / sizeof(cell_vertices[0]);
+
+  std::vector<CellData<dim> > cells (n_cells, CellData<dim>());
+  for (unsigned int i=0; i<n_cells; ++i) 
+    {
+      for (unsigned int j=0;
+           j<GeometryInfo<dim>::vertices_per_cell;
+           ++j)
+        cells[i].vertices[j] = cell_vertices[i][j];
+      cells[i].material_id = 0;
+    }
+
+  triangulation.create_triangulation (vertices,
+                                    cells,
+                                    SubCellData());
+  triangulation.refine_global (3);
+}
+
+
+
+template <>
+void LaplaceProblem<3>::create_coarse_grid ()
+{
+  GridGenerator::hyper_cube (triangulation);
+  triangulation.refine_global (1);
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::run () 
+{
+  for (unsigned int cycle=0; cycle<2; ++cycle)
+    {
+      deallog << "Cycle " << cycle << ':' << std::endl;
+
+      if (cycle == 0)
+       create_coarse_grid ();
+      else
+       refine_grid ();
+      
+
+      deallog << "   Number of active cells:       "
+               << triangulation.n_active_cells()
+               << std::endl;
+
+      Timer all;
+      all.start();
+      setup_system ();
+
+      deallog << "   Number of degrees of freedom: "
+               << dof_handler.n_dofs()
+               << std::endl;
+      deallog << "   Number of constraints       : "
+               << hanging_node_constraints.n_constraints()
+               << std::endl;
+      
+      assemble_system ();
+      solve ();
+      all.stop();
+
+    }
+}
+
+int main () 
+{
+  std::ofstream logfile(logname);
+  logfile.precision (3);
+  
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+  try
+    {
+      LaplaceProblem<3> laplace_problem_2d;
+      laplace_problem_2d.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.