]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add step-5.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 29 Dec 1999 22:05:31 +0000 (22:05 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 29 Dec 1999 22:05:31 +0000 (22:05 +0000)
git-svn-id: https://svn.dealii.org/trunk@2126 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/Attic/examples/step-by-step/step-4/step-4.cc
deal.II/deal.II/Attic/examples/step-by-step/step-5/Makefile [new file with mode: 0644]
deal.II/deal.II/Attic/examples/step-by-step/step-5/step-5.cc [new file with mode: 0644]
deal.II/doc/tutorial/chapter-2.step-by-step/navbar.html
deal.II/doc/tutorial/chapter-2.step-by-step/step-4.results
deal.II/doc/tutorial/chapter-2.step-by-step/step-5.intro [new file with mode: 0644]
deal.II/doc/tutorial/chapter-2.step-by-step/step-5.results [new file with mode: 0644]
deal.II/doc/tutorial/chapter-2.step-by-step/toc.html
deal.II/examples/step-4/step-4.cc
deal.II/examples/step-5/Makefile [new file with mode: 0644]
deal.II/examples/step-5/step-5.cc [new file with mode: 0644]

index cf073c55093720d8d5f7994a5db750797567a992..5c95b88071e7c923f2a92906a14ef4bd1cc42470 100644 (file)
@@ -252,7 +252,7 @@ void LaplaceProblem<dim>::make_grid_and_dofs ()
   cout << "   Number of active cells: "
        << triangulation.n_active_cells()
        << endl
-       << "  Total number of cells: "
+       << "   Total number of cells: "
        << triangulation.n_cells()
        << endl;
 
diff --git a/deal.II/deal.II/Attic/examples/step-by-step/step-5/Makefile b/deal.II/deal.II/Attic/examples/step-by-step/step-5/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-5/step-5.cc b/deal.II/deal.II/Attic/examples/step-by-step/step-5/step-5.cc
new file mode 100644 (file)
index 0000000..16e2c5c
--- /dev/null
@@ -0,0 +1,504 @@
+/* $Id$ */
+
+                                // The first few (many?) include
+                                // files have already been used in
+                                // the previous example, so we will
+                                // not explain their meaning here
+                                // again.
+#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 <fe/fe_lib.lagrange.h>
+#include <dofs/dof_tools.h>
+#include <fe/fe_values.h>
+#include <base/quadrature_lib.h>
+#include <base/function.h>
+#include <numerics/vectors.h>
+#include <numerics/matrices.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 <numerics/data_out.h>
+#include <fstream>
+
+#include <base/logstream.h>
+
+
+template <int dim>
+class LaplaceProblem 
+{
+  public:
+    LaplaceProblem ();
+    void run ();
+    
+  private:
+    void make_grid_and_dofs (const unsigned int refinement);
+    void assemble_system ();
+    void solve ();
+    void output_results ();
+    void clear ();
+
+    Triangulation<dim>   triangulation;
+    FEQ1<dim>            fe;
+    DoFHandler<dim>      dof_handler;
+
+    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 10;
+  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] = 10;
+    else
+      values[i] = 1;
+};
+
+
+
+template <int dim>
+LaplaceProblem<dim>::LaplaceProblem () :
+               dof_handler (triangulation)
+{};
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::make_grid_and_dofs (const unsigned int refinement)
+{
+  GridGenerator::hyper_cube (triangulation, -1, 1);
+  triangulation.refine_global (refinement);
+  
+  cout << "   Number of active cells: "
+       << triangulation.n_active_cells()
+       << endl
+       << "   Total number of cells: "
+       << triangulation.n_cells()
+       << endl;
+
+  dof_handler.distribute_dofs (fe);
+
+  cout << "   Number of degrees of freedom: "
+       << dof_handler.n_dofs()
+       << endl;
+
+  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);
+  sparsity_pattern.compress();
+
+  system_matrix.reinit (sparsity_pattern);
+
+  solution.reinit (dof_handler.n_dofs());
+  system_rhs.reinit (dof_handler.n_dofs());
+};
+
+
+
+                                // As in the previous examples, this
+                                // function is not changed much with
+                                // regard to its functionality, but
+                                // there are still some optimizations
+                                // which we will show. For this, it
+                                // is important to note that if
+                                // efficient solvers are used (such
+                                // as the preconditions CG method),
+                                // assembling the matrix and right
+                                // hand side can take a comparable
+                                // time, and it is worth the effort
+                                // to use one or two optimizations at
+                                // some places.
+                                //
+                                // What we will show here is how we
+                                // can avoid calls to the
+                                // shape_value, shape_grad, and
+                                // quadrature_point functions of the
+                                // FEValues object, and in particular
+                                // optimize away most of the virtual
+                                // function calls of the Function
+                                // object. The way to do so will be
+                                // explained in the following, while
+                                // those parts of this function that
+                                // are not changed with respect to
+                                // the previous example are not
+                                // commented on.
+template <int dim>
+void LaplaceProblem<dim>::assemble_system () 
+{  
+                                  // This time, we will again use a
+                                  // constant right hand side
+                                  // function, but a variable
+                                  // coefficient. The following
+                                  // object will be used for this:
+  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 ();
+
+                                      // As before, we want the
+                                      // FEValues object to compute
+                                      // the quantities which we told
+                                      // him to compute in the
+                                      // constructor using the update
+                                      // flags.
+      fe_values.reinit (cell);
+                                      // Now, these quantities are
+                                      // stored in arrays in the
+                                      // FEValues object. Usually,
+                                      // the internals of how and
+                                      // where they are stored is not
+                                      // something that the outside
+                                      // world should know, but since
+                                      // this is a time critical
+                                      // function we decided to
+                                      // publicize these arrays a
+                                      // little bit, and provide
+                                      // facilities to export the
+                                      // address where this data is
+                                      // stored.
+                                      //
+                                      // For example, the values of
+                                      // shape function j at
+                                      // quadrature point q is stored
+                                      // in a matrix, of which we can
+                                      // get the address as follows
+                                      // (note that this is a
+                                      // reference to the matrix,
+                                      // symbolized by the ampersand,
+                                      // and that it must be a
+                                      // constant reference, since
+                                      // only read-only access is
+                                      // granted):
+      const FullMatrix<double> 
+       & shape_values = fe_values.get_shape_values();
+                                      // Instead of writing
+                                      // fe_values.shape_value(j,q)
+                                      // we can now write
+                                      // shape_values(j,q), i.e. the
+                                      // function call needed
+                                      // previously for each access
+                                      // has been otimized away.
+                                      //
+                                      // There are alike functions
+                                      // for almost all data elements
+                                      // in the FEValues class. The
+                                      // gradient are accessed as
+                                      // follows:
+      const vector<vector<Tensor<1,dim> > >
+       & shape_grads  = fe_values.get_shape_grads();
+                                      // The data type looks a bit
+                                      // unwieldy, since each entry
+                                      // in the matrix (j,q) now
+                                      // needs to be the gradient of
+                                      // the shape function, which is
+                                      // a vector.
+                                      //
+                                      // Similarly, access to the
+                                      // place where quadrature
+                                      // points and the determinants
+                                      // of the Jacobian matrices
+                                      // times the weights of the
+                                      // respective quadrature points
+                                      // are stored, can be obtained
+                                      // like this:
+      const vector<double>
+       & JxW_values   = fe_values.get_JxW_values();
+      const vector<Point<dim> >
+       & q_points     = fe_values.get_quadrature_points();
+                                      // Admittedly, the declarations
+                                      // above are not easily
+                                      // readable, but they can save
+                                      // many function calls in the
+                                      // inner loops and can thus
+                                      // make assemblage faster.
+                                      //
+                                      // An additional advantage is
+                                      // that the inner loops are
+                                      // simpler to read, since the
+                                      // fe_values object is no more
+                                      // explicitely needed to access
+                                      // the different fields (see
+                                      // below). Unfortunately,
+                                      // things became a bit
+                                      // inconsistent, since the
+                                      // shape values are accessed
+                                      // via the FullMatrix operator
+                                      // (), i.e. using parentheses,
+                                      // while all the other fields
+                                      // are accessed through vector
+                                      // operator [], i.e. using
+                                      // brackets. This is due to
+                                      // historical reasons and
+                                      // frequently leads to a bit of
+                                      // confusion, but since the
+                                      // places where this happens
+                                      // are few in well-written
+                                      // programs, this is not too
+                                      // big a problem.
+
+                                      // There is one more thing: in
+                                      // this example, we want to use
+                                      // a non-constant
+                                      // coefficient. In the previous
+                                      // example, we have called the
+                                      // ``value'' function of the
+                                      // right hand side object for
+                                      // each quadrature
+                                      // point. Unfortunately, that
+                                      // is a virtual function, so
+                                      // calling it is relatively
+                                      // expensive. Therefore, we use
+                                      // a function of the Function
+                                      // class which returns the
+                                      // values at all quadrature
+                                      // points at once; that
+                                      // function is still virtual,
+                                      // but it needs to be computed
+                                      // once per cell only, not once
+                                      // in the inner loop:
+      coefficient.value_list (q_points, coefficient_values);
+                                      // It should be noted that the
+                                      // creation of the
+                                      // coefficient_values object is
+                                      // done outside the loop over
+                                      // all cells to avoid memory
+                                      // allocation each time we
+                                      // visit a new cell. Contrary
+                                      // to this, the other variables
+                                      // above were created inside
+                                      // the loop, but they were only
+                                      // references to memory that
+                                      // has already been allocated
+                                      // (i.e. they are pointers to
+                                      // that memory) and therefore,
+                                      // no new memory needs to be
+                                      // allocated; in particular, by
+                                      // declaring the pointers as
+                                      // close to their use as
+                                      // possible, we give the
+                                      // compiler a better choice to
+                                      // optimize them away
+                                      // altogether, something which
+                                      // it definitely can't do with
+                                      // the coefficient_values
+                                      // object since it is too
+                                      // complicated, but mostly
+                                      // because it's address is
+                                      // passed to a virtual function
+                                      // which is not knows at
+                                      // compile time.
+      
+                                      // Using the various
+                                      // abbreviations, the loops
+                                      // then look like this (the
+                                      // parentheses around the
+                                      // product of the two gradients
+                                      // are needed to indicate the
+                                      // dot product; we have to
+                                      // overrule associativity of
+                                      // the operator* here, since
+                                      // the compiler would otherwise
+                                      // complain about an undefined
+                                      // product of double*gradient
+                                      // since it parses
+                                      // left-to-right):
+      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]);
+
+                                            // For the right hand
+                                            // side, a constant value
+                                            // is used again:
+           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);
+       };
+    };
+
+                                  // Again use zero boundary values:
+  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);
+};
+
+
+
+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);
+
+  cout << "   " << solver_control.last_step()
+       << " CG iterations needed to obtain convergence."
+       << endl;
+};
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::output_results () 
+{
+  DataOut<dim> data_out;
+
+  data_out.attach_dof_handler (dof_handler);
+  data_out.add_data_vector (solution, "solution");
+
+  data_out.build_patches ();
+
+  ofstream output (dim == 2 ?
+                  "solution-2d.gmv" :
+                  "solution-3d.gmv");
+                                  // ...
+  data_out.write_gnuplot (output);
+};
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::clear () 
+{
+  system_rhs.reinit (0);
+  solution.reinit (0);
+  system_matrix.reinit ();
+  sparsity_pattern.reinit (0, 0, 0);
+  dof_handler.clear ();
+  triangulation.clear ();
+};
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::run () 
+{
+  cout << "Solving problem in " << dim << " space dimensions." << endl;
+  
+  for (unsigned int refinement=0; refinement<7; ++refinement)
+    {
+      cout << "Refinement step: " << refinement << endl;
+      
+      make_grid_and_dofs(refinement);
+      assemble_system ();
+      solve ();
+      output_results ();
+
+      clear ();
+    };
+};
+
+    
+
+int main () 
+{
+  deallog.depth_console (0);
+
+  LaplaceProblem<2> laplace_problem_2d;
+  laplace_problem_2d.run ();
+  
+  return 0;
+};
index fbb14abf656a2bf70a5f1d3b38464967ec6bc543..0124224506b9bb869d10726fb5ebca60ee52dbc8 100644 (file)
          <a href="step-4.html" target="body">Step 4</a>
        </p>
       </li>
+
+      <li>
+       <p>
+         <a href="step-5.html" target="body">Step 5</a>
+       </p>
+      </li>
     </ol>
     <p>
       <a href="../index.html" target="_top">Back to the tutorial index</a>
index 3201c10650d3888808d27e0a5d23b14b7063299a..7514fac593112471e37462d58e5b4b8a1eb7660b 100644 (file)
@@ -1,7 +1,10 @@
 <a name="Results"></a>
 <h1>Results</h1>
 
-The output of the program looks as follows:
+The output of the program looks as follows (the number of iterations
+may vary by one or two, depending on your computer, since this is
+often dependent on the round-off accuracy of floating point
+operations, which differs between micro-processors):
 <pre>
 <code>
 Solving problem in 2 space dimensions.
diff --git a/deal.II/doc/tutorial/chapter-2.step-by-step/step-5.intro b/deal.II/doc/tutorial/chapter-2.step-by-step/step-5.intro
new file mode 100644 (file)
index 0000000..1518606
--- /dev/null
@@ -0,0 +1,3 @@
+<a name="Intro"></a>
+<h1>Introduction</h1>
+
diff --git a/deal.II/doc/tutorial/chapter-2.step-by-step/step-5.results b/deal.II/doc/tutorial/chapter-2.step-by-step/step-5.results
new file mode 100644 (file)
index 0000000..747bdbb
--- /dev/null
@@ -0,0 +1,3 @@
+<a name="Results"></a>
+<h1>Results</h1>
+
index 32b036a221306f2763b1622be7c25d4f595dfbf7..0bec85471b40b549908697c2f291fc34c7a31ce9 100644 (file)
@@ -65,6 +65,13 @@ At present, the following programs exist:
       3D, although the program is exactly the same. Non-constant right
       hand side function. Non-homogeneous boundary values.
       </dd>
+
+      <dt><a href="step-5.html">Step 5</a></dt> 
+      <dd><strong>What's new:</strong> Computations on successively
+      refined grids. Some optimizations.  Non-constant coefficient in
+      the elliptic operator. Preconditioning of the CG solver for the
+      linear system of equations.
+      </dd>
 </dl>
 
 <!-- Page Foot -->
index cf073c55093720d8d5f7994a5db750797567a992..5c95b88071e7c923f2a92906a14ef4bd1cc42470 100644 (file)
@@ -252,7 +252,7 @@ void LaplaceProblem<dim>::make_grid_and_dofs ()
   cout << "   Number of active cells: "
        << triangulation.n_active_cells()
        << endl
-       << "  Total number of cells: "
+       << "   Total number of cells: "
        << triangulation.n_cells()
        << endl;
 
diff --git a/deal.II/examples/step-5/Makefile b/deal.II/examples/step-5/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-5/step-5.cc b/deal.II/examples/step-5/step-5.cc
new file mode 100644 (file)
index 0000000..16e2c5c
--- /dev/null
@@ -0,0 +1,504 @@
+/* $Id$ */
+
+                                // The first few (many?) include
+                                // files have already been used in
+                                // the previous example, so we will
+                                // not explain their meaning here
+                                // again.
+#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 <fe/fe_lib.lagrange.h>
+#include <dofs/dof_tools.h>
+#include <fe/fe_values.h>
+#include <base/quadrature_lib.h>
+#include <base/function.h>
+#include <numerics/vectors.h>
+#include <numerics/matrices.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 <numerics/data_out.h>
+#include <fstream>
+
+#include <base/logstream.h>
+
+
+template <int dim>
+class LaplaceProblem 
+{
+  public:
+    LaplaceProblem ();
+    void run ();
+    
+  private:
+    void make_grid_and_dofs (const unsigned int refinement);
+    void assemble_system ();
+    void solve ();
+    void output_results ();
+    void clear ();
+
+    Triangulation<dim>   triangulation;
+    FEQ1<dim>            fe;
+    DoFHandler<dim>      dof_handler;
+
+    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 10;
+  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] = 10;
+    else
+      values[i] = 1;
+};
+
+
+
+template <int dim>
+LaplaceProblem<dim>::LaplaceProblem () :
+               dof_handler (triangulation)
+{};
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::make_grid_and_dofs (const unsigned int refinement)
+{
+  GridGenerator::hyper_cube (triangulation, -1, 1);
+  triangulation.refine_global (refinement);
+  
+  cout << "   Number of active cells: "
+       << triangulation.n_active_cells()
+       << endl
+       << "   Total number of cells: "
+       << triangulation.n_cells()
+       << endl;
+
+  dof_handler.distribute_dofs (fe);
+
+  cout << "   Number of degrees of freedom: "
+       << dof_handler.n_dofs()
+       << endl;
+
+  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);
+  sparsity_pattern.compress();
+
+  system_matrix.reinit (sparsity_pattern);
+
+  solution.reinit (dof_handler.n_dofs());
+  system_rhs.reinit (dof_handler.n_dofs());
+};
+
+
+
+                                // As in the previous examples, this
+                                // function is not changed much with
+                                // regard to its functionality, but
+                                // there are still some optimizations
+                                // which we will show. For this, it
+                                // is important to note that if
+                                // efficient solvers are used (such
+                                // as the preconditions CG method),
+                                // assembling the matrix and right
+                                // hand side can take a comparable
+                                // time, and it is worth the effort
+                                // to use one or two optimizations at
+                                // some places.
+                                //
+                                // What we will show here is how we
+                                // can avoid calls to the
+                                // shape_value, shape_grad, and
+                                // quadrature_point functions of the
+                                // FEValues object, and in particular
+                                // optimize away most of the virtual
+                                // function calls of the Function
+                                // object. The way to do so will be
+                                // explained in the following, while
+                                // those parts of this function that
+                                // are not changed with respect to
+                                // the previous example are not
+                                // commented on.
+template <int dim>
+void LaplaceProblem<dim>::assemble_system () 
+{  
+                                  // This time, we will again use a
+                                  // constant right hand side
+                                  // function, but a variable
+                                  // coefficient. The following
+                                  // object will be used for this:
+  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 ();
+
+                                      // As before, we want the
+                                      // FEValues object to compute
+                                      // the quantities which we told
+                                      // him to compute in the
+                                      // constructor using the update
+                                      // flags.
+      fe_values.reinit (cell);
+                                      // Now, these quantities are
+                                      // stored in arrays in the
+                                      // FEValues object. Usually,
+                                      // the internals of how and
+                                      // where they are stored is not
+                                      // something that the outside
+                                      // world should know, but since
+                                      // this is a time critical
+                                      // function we decided to
+                                      // publicize these arrays a
+                                      // little bit, and provide
+                                      // facilities to export the
+                                      // address where this data is
+                                      // stored.
+                                      //
+                                      // For example, the values of
+                                      // shape function j at
+                                      // quadrature point q is stored
+                                      // in a matrix, of which we can
+                                      // get the address as follows
+                                      // (note that this is a
+                                      // reference to the matrix,
+                                      // symbolized by the ampersand,
+                                      // and that it must be a
+                                      // constant reference, since
+                                      // only read-only access is
+                                      // granted):
+      const FullMatrix<double> 
+       & shape_values = fe_values.get_shape_values();
+                                      // Instead of writing
+                                      // fe_values.shape_value(j,q)
+                                      // we can now write
+                                      // shape_values(j,q), i.e. the
+                                      // function call needed
+                                      // previously for each access
+                                      // has been otimized away.
+                                      //
+                                      // There are alike functions
+                                      // for almost all data elements
+                                      // in the FEValues class. The
+                                      // gradient are accessed as
+                                      // follows:
+      const vector<vector<Tensor<1,dim> > >
+       & shape_grads  = fe_values.get_shape_grads();
+                                      // The data type looks a bit
+                                      // unwieldy, since each entry
+                                      // in the matrix (j,q) now
+                                      // needs to be the gradient of
+                                      // the shape function, which is
+                                      // a vector.
+                                      //
+                                      // Similarly, access to the
+                                      // place where quadrature
+                                      // points and the determinants
+                                      // of the Jacobian matrices
+                                      // times the weights of the
+                                      // respective quadrature points
+                                      // are stored, can be obtained
+                                      // like this:
+      const vector<double>
+       & JxW_values   = fe_values.get_JxW_values();
+      const vector<Point<dim> >
+       & q_points     = fe_values.get_quadrature_points();
+                                      // Admittedly, the declarations
+                                      // above are not easily
+                                      // readable, but they can save
+                                      // many function calls in the
+                                      // inner loops and can thus
+                                      // make assemblage faster.
+                                      //
+                                      // An additional advantage is
+                                      // that the inner loops are
+                                      // simpler to read, since the
+                                      // fe_values object is no more
+                                      // explicitely needed to access
+                                      // the different fields (see
+                                      // below). Unfortunately,
+                                      // things became a bit
+                                      // inconsistent, since the
+                                      // shape values are accessed
+                                      // via the FullMatrix operator
+                                      // (), i.e. using parentheses,
+                                      // while all the other fields
+                                      // are accessed through vector
+                                      // operator [], i.e. using
+                                      // brackets. This is due to
+                                      // historical reasons and
+                                      // frequently leads to a bit of
+                                      // confusion, but since the
+                                      // places where this happens
+                                      // are few in well-written
+                                      // programs, this is not too
+                                      // big a problem.
+
+                                      // There is one more thing: in
+                                      // this example, we want to use
+                                      // a non-constant
+                                      // coefficient. In the previous
+                                      // example, we have called the
+                                      // ``value'' function of the
+                                      // right hand side object for
+                                      // each quadrature
+                                      // point. Unfortunately, that
+                                      // is a virtual function, so
+                                      // calling it is relatively
+                                      // expensive. Therefore, we use
+                                      // a function of the Function
+                                      // class which returns the
+                                      // values at all quadrature
+                                      // points at once; that
+                                      // function is still virtual,
+                                      // but it needs to be computed
+                                      // once per cell only, not once
+                                      // in the inner loop:
+      coefficient.value_list (q_points, coefficient_values);
+                                      // It should be noted that the
+                                      // creation of the
+                                      // coefficient_values object is
+                                      // done outside the loop over
+                                      // all cells to avoid memory
+                                      // allocation each time we
+                                      // visit a new cell. Contrary
+                                      // to this, the other variables
+                                      // above were created inside
+                                      // the loop, but they were only
+                                      // references to memory that
+                                      // has already been allocated
+                                      // (i.e. they are pointers to
+                                      // that memory) and therefore,
+                                      // no new memory needs to be
+                                      // allocated; in particular, by
+                                      // declaring the pointers as
+                                      // close to their use as
+                                      // possible, we give the
+                                      // compiler a better choice to
+                                      // optimize them away
+                                      // altogether, something which
+                                      // it definitely can't do with
+                                      // the coefficient_values
+                                      // object since it is too
+                                      // complicated, but mostly
+                                      // because it's address is
+                                      // passed to a virtual function
+                                      // which is not knows at
+                                      // compile time.
+      
+                                      // Using the various
+                                      // abbreviations, the loops
+                                      // then look like this (the
+                                      // parentheses around the
+                                      // product of the two gradients
+                                      // are needed to indicate the
+                                      // dot product; we have to
+                                      // overrule associativity of
+                                      // the operator* here, since
+                                      // the compiler would otherwise
+                                      // complain about an undefined
+                                      // product of double*gradient
+                                      // since it parses
+                                      // left-to-right):
+      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]);
+
+                                            // For the right hand
+                                            // side, a constant value
+                                            // is used again:
+           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);
+       };
+    };
+
+                                  // Again use zero boundary values:
+  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);
+};
+
+
+
+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);
+
+  cout << "   " << solver_control.last_step()
+       << " CG iterations needed to obtain convergence."
+       << endl;
+};
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::output_results () 
+{
+  DataOut<dim> data_out;
+
+  data_out.attach_dof_handler (dof_handler);
+  data_out.add_data_vector (solution, "solution");
+
+  data_out.build_patches ();
+
+  ofstream output (dim == 2 ?
+                  "solution-2d.gmv" :
+                  "solution-3d.gmv");
+                                  // ...
+  data_out.write_gnuplot (output);
+};
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::clear () 
+{
+  system_rhs.reinit (0);
+  solution.reinit (0);
+  system_matrix.reinit ();
+  sparsity_pattern.reinit (0, 0, 0);
+  dof_handler.clear ();
+  triangulation.clear ();
+};
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::run () 
+{
+  cout << "Solving problem in " << dim << " space dimensions." << endl;
+  
+  for (unsigned int refinement=0; refinement<7; ++refinement)
+    {
+      cout << "Refinement step: " << refinement << endl;
+      
+      make_grid_and_dofs(refinement);
+      assemble_system ();
+      solve ();
+      output_results ();
+
+      clear ();
+    };
+};
+
+    
+
+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.