]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Linear transport in non-conservative form.
authorhartmann <hartmann@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 9 Nov 2001 13:14:19 +0000 (13:14 +0000)
committerhartmann <hartmann@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 9 Nov 2001 13:14:19 +0000 (13:14 +0000)
git-svn-id: https://svn.dealii.org/trunk@5187 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-12/step-12.cc [new file with mode: 0644]

diff --git a/deal.II/examples/step-12/step-12.cc b/deal.II/examples/step-12/step-12.cc
new file mode 100644 (file)
index 0000000..ff3e5ee
--- /dev/null
@@ -0,0 +1,716 @@
+/* $Id$ */
+/* Author: Ralf Hartmann, University of Heidelberg, 2000 */
+
+                                // The first few files have already
+                                // been covered in previous examples
+                                // and will thus not be further
+                                // commented on.
+#include <base/quadrature_lib.h>
+#include <base/function.h>
+#include <lac/vector.h>
+#include <lac/sparse_matrix.h>
+#include <lac/vector_memory.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_values.h>
+#include <numerics/data_out.h>
+#include <grid/grid_out.h>
+#include <grid/grid_refinement.h>
+#include <numerics/error_estimator.h>
+
+#include <fe/mapping_q1.h>
+#include <fe/fe_dgq.h>
+#include <lac/precondition_block.h>
+#include <lac/solver_richardson.h>
+
+
+#include <fstream>
+
+
+template <int dim>
+class Beta
+{
+  public:
+    Beta () {};
+
+    void value_list (const std::vector<Point<dim> > &points,
+                    std::vector<Point<dim> > &values) const;
+};
+
+
+template <int dim>
+class RHS:  public Function<dim>
+{
+  public:
+    RHS() {};
+    
+    virtual void value_list (const std::vector<Point<dim> > &points,
+                            std::vector<double> &values,
+                            const unsigned int component=0) const;
+};
+
+
+template <int dim>
+class BoundaryFunction:  public Function<dim>
+{
+  public:
+    BoundaryFunction() {};
+    
+    virtual void value_list (const std::vector<Point<dim> > &points,
+                            std::vector<double> &values,
+                            const unsigned int component=0) const;
+};
+
+template <int dim>
+class DGAssembler
+{
+  public:
+    DGAssembler() {};
+
+    void assemble_cell_term(const FEValuesBase<dim>& fe_v,
+                           FullMatrix<double> &cell_matrix,
+                           Vector<double> &cell_vector);
+    
+    void assemble_face_term(const FEFaceValuesBase<dim>& fe_v,
+                           const FEFaceValuesBase<dim>& fe_v_neighbor,
+                           FullMatrix<double> &cell_matrix,
+                           FullMatrix<double> &cell_inflow_matrix,
+                           Vector<double> &cell_vector);
+    
+  private:
+    Beta<dim> beta_function;
+    RHS<dim> rhs_function;
+    BoundaryFunction<dim> boundary_function;
+};
+
+                                // The main class is again almost
+                                // unchanged. Two additions, however,
+                                // are made: we have added the
+                                // ``refine'' function, which is used
+                                // to adaptively refine the grid
+                                // (instead of the global refinement
+                                // in the previous examples), and a
+                                // variable which will hold the
+                                // constraints associated to the
+                                // hanging nodes.
+template <int dim>
+class TransportProblem
+{
+  public:
+    TransportProblem ();
+    ~TransportProblem ();
+
+    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;
+    MappingQ1<dim>       mapping;
+    
+                                    // We need a finite element
+                                    // again. This time, we will want
+                                    // to use quadratic polynomials
+                                    // (but this is only specified in
+                                    // the constructor):
+    FE_DGQ<dim>          fe;
+    DoFHandler<dim>      dof_handler;
+
+    SparsityPattern      sparsity_pattern;
+    SparseMatrix<double> system_matrix;
+    
+    Vector<double>       solution;
+    Vector<double>       right_hand_side;
+
+    DGAssembler<dim>     dg_assembler;
+};
+
+
+template <>
+void Beta<2>::value_list(const std::vector<Point<2> > &points,
+                        std::vector<Point<2> > &values) const
+{
+  Assert(values.size()==points.size(), ExcDimensionMismatch(values.size(),points.size()));
+  for (unsigned int i=0; i<points.size(); ++i)
+    {
+      const Point<2> &p=points[i];
+      Point<2> &beta=values[i];
+
+      beta(0)=-p(1);
+      beta(1)=p(0);
+      beta/=sqrt(beta*beta);
+    }
+}
+
+
+
+
+template <int dim>
+void RHS<dim>::value_list(const std::vector<Point<dim> > &,
+                         std::vector<double> &values,
+                         const unsigned int) const
+{
+  for (unsigned int i=0; i<values.size(); ++i)
+    values[i]=0;
+}
+
+
+
+
+template <int dim>
+void BoundaryFunction<dim>::value_list(const std::vector<Point<dim> > &points,
+                                      std::vector<double> &values,
+                                      const unsigned int) const
+{
+  Assert(values.size()==points.size(), ExcDimensionMismatch(values.size(),points.size()));
+  for (unsigned int i=0; i<values.size(); ++i)
+    {
+      if (points[i](0)<0.5)
+       values[i]=1.;
+      else
+       values[i]=0.;
+    }
+}
+
+
+
+
+template <int dim>
+void DGAssembler<dim>::assemble_cell_term(const FEValuesBase<dim>& fe_v,
+                                         FullMatrix<double> &cell_matrix,
+                                         Vector<double> &cell_vector)
+{
+  const vector<vector<Tensor<1,2> > > &grad_v = fe_v.get_shape_grads ();
+  const FullMatrix<double> &v = fe_v.get_shape_values ();
+  const vector<double> &JxW = fe_v.get_JxW_values ();
+
+  vector<Point<dim> > beta (fe_v.n_quadrature_points);
+  vector<double> rhs (fe_v.n_quadrature_points);
+  
+  beta_function.value_list (fe_v.get_quadrature_points(), beta);
+  rhs_function.value_list (fe_v.get_quadrature_points(), rhs);
+  
+  for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
+    for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i) 
+      {
+       for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
+         cell_matrix(i,j) += beta[point]*grad_v[j][point]*
+                             v(i,point) *
+                             JxW[point];
+       
+       cell_vector(i) += rhs[point] *v(i,point) * JxW[point];
+      }
+}
+
+
+template <int dim>
+void DGAssembler<dim>::assemble_face_term(const FEFaceValuesBase<dim>& fe_v,
+                                         const FEFaceValuesBase<dim>& fe_v_neighbor,      
+                                         FullMatrix<double> &cell_matrix,
+                                         FullMatrix<double> &cell_inflow_matrix,
+                                         Vector<double> &cell_vector)
+{
+  DoFHandler<dim>::face_iterator face=fe_v.get_face();
+  
+  const FullMatrix<double> &v = fe_v.get_shape_values ();
+  const FullMatrix<double> &v_neighbor = fe_v_neighbor.get_shape_values ();  
+  const vector<double> &JxW = fe_v.get_JxW_values ();
+  const vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
+
+  vector<Point<dim> > beta (fe_v.n_quadrature_points);
+  vector<double> g(fe_v.n_quadrature_points);
+  
+  beta_function.value_list (fe_v.get_quadrature_points(), beta);
+
+  if (face->at_boundary())
+    boundary_function.value_list (fe_v.get_quadrature_points(), g);
+
+  for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
+    {
+      double beta_n=beta[point] * normals[point];
+      if (beta_n<0)
+       {
+         for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
+           {
+             for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
+               cell_matrix(i,j) -= beta_n *
+                                   v(j,point) *
+                                   v(i,point) *
+                                   JxW[point];
+
+             if (!face->at_boundary())
+               for (unsigned int k=0; k<fe_v_neighbor.dofs_per_cell; ++k)
+                 cell_inflow_matrix(i,k) += beta_n *
+                                            v_neighbor(k,point) *
+                                            v(i,point) *
+                                            JxW[point];
+             else
+               cell_vector(i) -= beta_n *
+                                 g[point] *
+                                 v(i,point) *
+                                 JxW[point];
+           }
+       }
+    }
+}
+
+  
+template <int dim>
+TransportProblem<dim>::TransportProblem () :
+                fe (1),
+               dof_handler (triangulation)
+{}
+
+
+template <int dim>
+TransportProblem<dim>::~TransportProblem () 
+{
+  dof_handler.clear ();
+};
+
+
+
+template <int dim>
+void TransportProblem<dim>::setup_system ()
+{
+                                  // To distribute degrees of
+                                  // freedom, the ``dof_handler''
+                                  // variable takes only the finite
+                                  // element object. In this case, it
+                                  // will distribute four degrees of
+                                  // freedom per cell.
+  dof_handler.distribute_dofs (fe);
+
+  sparsity_pattern.reinit (dof_handler.n_dofs(),
+                          dof_handler.n_dofs(),
+                          dof_handler.max_couplings_between_dofs());
+  DoFTools::make_flux_sparsity_pattern (dof_handler, sparsity_pattern);
+  sparsity_pattern.compress();
+  
+  system_matrix.reinit (sparsity_pattern);
+
+  solution.reinit (dof_handler.n_dofs());
+  right_hand_side.reinit (dof_handler.n_dofs());
+};
+
+
+
+template <int dim>
+void TransportProblem<dim>::assemble_system () 
+{
+                                  // See Cockburn paper for the proper quadrature.
+  QGauss2<dim>  quadrature;
+  QGauss2<dim-1>  face_quadrature;
+  
+  const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
+  vector<unsigned int> dofs (dofs_per_cell);
+  vector<unsigned int> dofs_neighbor (dofs_per_cell);
+
+  UpdateFlags update_flags = UpdateFlags(update_values
+                                        | update_gradients
+                                        | update_q_points
+                                        | update_JxW_values);
+  
+  UpdateFlags face_update_flags = UpdateFlags(update_values
+                                             | update_q_points
+                                             | update_JxW_values
+                                             | update_normal_vectors);
+   
+
+  FEValues<dim> fe_v (
+    mapping, fe, quadrature, update_flags);
+  FEFaceValues<dim> fe_v_face (
+    mapping, fe, face_quadrature, face_update_flags);
+  FESubfaceValues<dim> fe_v_subface (
+    mapping, fe, face_quadrature, face_update_flags);
+  FEFaceValues<dim> fe_v_face_neighbor (
+    mapping, fe, face_quadrature, UpdateFlags(update_values | update_default));
+  FESubfaceValues<dim> fe_v_subface_neighbor (
+    mapping, fe, face_quadrature, UpdateFlags(update_values | update_default));
+
+                                  // includes the u and v terms
+  FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
+                                  // includes u_hat and v terms
+  FullMatrix<double> cell_inflow_matrix (dofs_per_cell, dofs_per_cell);
+
+  Vector<double>  cell_vector (dofs_per_cell);
+
+  DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
+                                       endc = dof_handler.end();
+  DoFHandler<dim>::face_iterator face;
+  DoFHandler<dim>::cell_iterator neighbor;
+  DoFHandler<dim>::cell_iterator neighbor_child;
+
+  for (;cell!=endc; ++cell) 
+    {
+                                      // re-init fe values for this cell
+      fe_v.reinit (cell);
+
+      cell_matrix.clear ();
+      cell_vector.clear ();
+
+      dg_assembler.assemble_cell_term(fe_v,
+                                     cell_matrix,
+                                     cell_vector);
+      
+      cell->get_dof_indices (dofs);
+
+      for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
+       {
+         face = cell->face(face_no);
+         
+         cell_inflow_matrix.clear();
+
+         if (face->at_boundary())
+           {
+             fe_v_face.reinit (cell, face_no);
+
+             dg_assembler.assemble_face_term(fe_v_face,
+                                             fe_v_face,
+                                             cell_matrix,
+                                             cell_inflow_matrix,
+                                             cell_vector);
+           }
+         else // if (!face->at_boundary())
+           {
+             Assert (cell->neighbor(face_no).state() == valid, ExcInternalError());
+             neighbor = cell->neighbor(face_no);
+             
+             if (face->has_children())  // i.e. neighbor is one level more refined than cell
+               {
+                                                  // store which number #cell# has in the
+                                                  // list of neighbors of #neighbor#
+                 const unsigned int neighbor2=cell->neighbor_of_neighbor(face_no);
+                 
+                 
+                                                  // loop over all subfaces
+                 for (unsigned int subface_no=0; subface_no<GeometryInfo<dim>::subfaces_per_face;
+                      ++subface_no)
+                   {
+                                                      // get an iterator pointing to the
+                                                      // cell behind the present subface
+                     neighbor_child = neighbor->child(GeometryInfo<dim>::
+                                                      child_cell_on_face(neighbor2,subface_no));
+                     Assert (neighbor_child->face(neighbor2) == face->child(subface_no),
+                             ExcInternalError());
+                     Assert (!neighbor_child->has_children(), ExcInternalError());
+
+                     fe_v_subface.reinit (cell, face_no, subface_no);
+                     fe_v_face_neighbor.reinit (neighbor_child, neighbor2);
+       
+                     cell_inflow_matrix.clear();
+           
+                     dg_assembler.assemble_face_term(fe_v_subface,
+                                                     fe_v_face_neighbor,
+                                                     cell_matrix,
+                                                     cell_inflow_matrix,
+                                                     cell_vector);
+                 
+                                                      // get indices of dofs of neighbor_child cell
+                     neighbor_child->get_dof_indices (dofs_neighbor);
+                                                               
+                                                      // distribute cell matrix
+                     for (unsigned int i=0; i<dofs_per_cell; ++i)
+                       for (unsigned int k=0; k<dofs_per_cell; ++k)
+                         system_matrix.add(dofs[i], dofs_neighbor[k],
+                                           cell_inflow_matrix(i,k));
+                   }
+               }
+             else // if (!face->has_children())
+               {
+                 if (neighbor->level() == cell->level()) 
+                   {
+                                                      // store which number #cell# has in the
+                                                      // list of neighbors of #neighbor#
+                     const unsigned int neighbor2=cell->neighbor_of_neighbor(face_no);
+
+                     fe_v_face.reinit (cell, face_no);
+                     fe_v_face_neighbor.reinit (neighbor, neighbor2);
+                     
+                     dg_assembler.assemble_face_term(fe_v_face,
+                                                     fe_v_face_neighbor,
+                                                     cell_matrix,
+                                                     cell_inflow_matrix,
+                                                     cell_vector);
+                   }
+                 else // if (neighbor->level() < cell->level()) i.e. neighbor is one level coarser than cell
+                   {
+                     Assert(neighbor->level() < cell->level(), ExcInternalError());
+
+                     const std::pair<unsigned int, unsigned int> faceno_subfaceno=
+                       cell->neighbor_of_coarser_neighbor(face_no);
+                     const unsigned int neighbor_face_no=faceno_subfaceno.first,
+                                     neighbor_subface_no=faceno_subfaceno.second;
+
+                     Assert (neighbor->neighbor(neighbor_face_no)
+                             ->child(GeometryInfo<dim>::child_cell_on_face(
+                               face_no,neighbor_subface_no)) == cell, ExcInternalError());
+                       
+                                                      // now 'neighbor_face_no' stores the number
+                                                      // of a face in the list of faces of 'neighbor'.
+                                                      // This face has got a subface that is 
+                                                      // between 'cell' and 'neighbor'.
+                                                      // 'neighbor_subface_no' stores the number
+                                                      // of this subface in the list of subfaces of this
+                                                      // face 'neighbor->face(neighbor_face_no)'
+                                                      // that is between 'cell' and 'neighbor'
+                     fe_v_face.reinit (cell, face_no);
+                     fe_v_subface_neighbor.reinit (neighbor, neighbor_face_no,
+                                                   neighbor_subface_no);
+                     
+                     dg_assembler.assemble_face_term(fe_v_face,
+                                                     fe_v_subface_neighbor,
+                                                     cell_matrix,
+                                                     cell_inflow_matrix,
+                                                     cell_vector);
+                   } // else // if (neighbor->level() < cell->level())
+
+                                                  // get indices of dofs of neighbor_child cell
+                 neighbor->get_dof_indices (dofs_neighbor);
+                                                               
+                                                  // distribute cell_inflow_matrix
+                 for (unsigned int i=0; i<dofs_per_cell; ++i)
+                   for (unsigned int k=0; k<dofs_per_cell; ++k)
+                     system_matrix.add(dofs[i], dofs_neighbor[k],
+                                       cell_inflow_matrix(i,k));
+               } // else // if (!face->has_children())
+           }  // else // if (!face->at_boundary())
+       } //for (face_no...)
+      
+                                      // distribute cell matrix
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+       for (unsigned int j=0; j<dofs_per_cell; ++j)
+         system_matrix.add(dofs[i], dofs[j], cell_matrix(i,j));
+      
+                                      // distribute cell vector
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+       right_hand_side(dofs[i]) += cell_vector(i);
+    }  // for (cell...)
+};
+
+
+
+template <int dim>
+void TransportProblem<dim>::solve () 
+{  
+  SolverControl           solver_control (1000, 1e-12);
+  PrimitiveVectorMemory<> vector_memory;
+  SolverRichardson<>      solver (solver_control, vector_memory);
+
+  PreconditionBlockSSOR<double> preconditioner;
+  preconditioner.initialize(system_matrix, fe.dofs_per_cell);
+  preconditioner.invert_diagblocks();
+  
+  solver.solve (system_matrix, solution, right_hand_side,
+               preconditioner);
+};
+
+
+template <int dim>
+void TransportProblem<dim>::refine_grid ()
+{
+  Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+
+  FunctionMap<dim>::type neumann_boundary;
+
+  KellyErrorEstimator<dim>::estimate (dof_handler,
+                                     QGauss3<dim-1>(),
+                                     neumann_boundary,
+                                     solution,
+                                     estimated_error_per_cell);
+
+  GridRefinement::refine_and_coarsen_fixed_number (triangulation,
+                                                  estimated_error_per_cell,
+                                                  0.3, 0.03);
+
+  triangulation.execute_coarsening_and_refinement ();
+};
+
+
+
+template <int dim>
+void TransportProblem<dim>::output_results (const unsigned int cycle) const
+{
+                                  // We want to write the grid in
+                                  // each cycle. Here is another way
+                                  // to quickly produce a filename
+                                  // based on the cycle number. It
+                                  // assumes that the numbers `0'
+                                  // through `9' are represented
+                                  // consecutively in the character
+                                  // set (which is the case in all
+                                  // known character sets). However,
+                                  // this will only work if the cycle
+                                  // number is less than ten, which
+                                  // we check by an assertion.
+  std::string filename = "grid-";
+  filename += ('0' + cycle);
+  Assert (cycle < 10, ExcInternalError());
+  
+  filename += ".eps";
+  std::ofstream eps_output (filename.c_str());
+
+                                  // Using this filename, we write
+                                  // each grid as a postscript file.
+  GridOut grid_out;
+  grid_out.write_eps (triangulation, eps_output);
+
+                                  // output of the solution
+  filename = "sol-";
+  filename += ('0' + cycle);
+  Assert (cycle < 10, ExcInternalError());
+  
+  filename += ".gnuplot";
+  std::ofstream gnuplot_output (filename.c_str());
+  
+  DataOut<dim> data_out;
+  data_out.attach_dof_handler (dof_handler);
+  data_out.add_data_vector (solution, "u");
+
+  data_out.build_patches ();
+  
+  data_out.write_gnuplot(gnuplot_output);
+};
+
+
+
+template <int dim>
+void TransportProblem<dim>::run () 
+{
+  for (unsigned int cycle=0; cycle<3; ++cycle)
+    {
+      std::cout << "Cycle " << cycle << ':' << std::endl;
+
+      if (cycle == 0)
+       {
+         GridGenerator::hyper_cube (triangulation);
+
+         triangulation.refine_global (3);
+       }
+      else
+                                        // In case this is not the
+                                        // first cycle, we want to
+                                        // refine the grid. Unlike
+                                        // the global refinement
+                                        // employed in the last
+                                        // example, we now use the
+                                        // adaptive procedure
+                                        // described in the function
+                                        // which we now call:
+       {
+         refine_grid ();
+       };
+      
+
+      std::cout << "   Number of active cells:       "
+               << triangulation.n_active_cells()
+               << std::endl;
+
+      setup_system ();
+
+      std::cout << "   Number of degrees of freedom: "
+               << dof_handler.n_dofs()
+               << std::endl;
+      
+      assemble_system ();
+      solve ();
+      output_results (cycle);
+    }
+}
+
+int main () 
+{
+
+                                  // The general idea behind the
+                                  // layout of this function is as
+                                  // follows: let's try to run the
+                                  // program as we did before...
+  try
+    {
+      TransportProblem<2> Transport_problem_2d;
+      Transport_problem_2d.run ();
+    }
+                                  // ...and if this should fail, try
+                                  // to gather as much information as
+                                  // possible. Specifically, if the
+                                  // exception that was thrown is an
+                                  // object of a class that is
+                                  // derived from the C++ standard
+                                  // class ``exception'', then we can
+                                  // use the ``what'' member function
+                                  // to get a string which describes
+                                  // the reason why the exception was
+                                  // thrown. 
+                                  //
+                                  // The deal.II exception classes
+                                  // are all derived from the
+                                  // standard class, and in
+                                  // particular, the ``exc.what()''
+                                  // function will return
+                                  // approximately the same string as
+                                  // would be generated if the
+                                  // exception was thrown using the
+                                  // ``Assert'' macro. You have seen
+                                  // the output of such an exception
+                                  // in the previous example, and you
+                                  // then know that it contains the
+                                  // file and line number of where
+                                  // the exception occured, and some
+                                  // other information. This is also
+                                  // what would be printed in the
+                                  // following.
+  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;
+                                      // We can't do much more than
+                                      // printing as much information
+                                      // as we can get to, so abort
+                                      // with error:
+      return 1;
+    }
+                                  // If the exception that was thrown
+                                  // somewhere was not an object of a
+                                  // class derived from the standard
+                                  // ``exception'' class, then we
+                                  // can't do anything at all. We
+                                  // then simply print an error
+                                  // message and exit.
+  catch (...) 
+    {
+      std::cerr << std::endl << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+               << "Aborting!" << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      return 1;
+    };
+
+                                  // If we got to this point, there
+                                  // was no exception which
+                                  // propagated up to the main
+                                  // function (maybe there were some,
+                                  // but they were caught somewhere
+                                  // in the program or the
+                                  // library). Therefore, the program
+                                  // performed as was expected and we
+                                  // can return without error.
+  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.