]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add step-4.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 22 Dec 1999 21:28:11 +0000 (21:28 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 22 Dec 1999 21:28:11 +0000 (21:28 +0000)
git-svn-id: https://svn.dealii.org/trunk@2112 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/Attic/examples/step-by-step/step-4/Makefile [new file with mode: 0644]
deal.II/deal.II/Attic/examples/step-by-step/step-4/step-4.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.data/grid-3d.jpg [new file with mode: 0644]
deal.II/doc/tutorial/chapter-2.step-by-step/step-4.data/solution-3d.jpg [new file with mode: 0644]
deal.II/doc/tutorial/chapter-2.step-by-step/step-4.intro [new file with mode: 0644]
deal.II/doc/tutorial/chapter-2.step-by-step/step-4.results [new file with mode: 0644]
deal.II/doc/tutorial/chapter-2.step-by-step/toc.html
deal.II/examples/step-4/Makefile [new file with mode: 0644]
deal.II/examples/step-4/step-4.cc [new file with mode: 0644]

diff --git a/deal.II/deal.II/Attic/examples/step-by-step/step-4/Makefile b/deal.II/deal.II/Attic/examples/step-by-step/step-4/Makefile
new file mode 100644 (file)
index 0000000..9212069
--- /dev/null
@@ -0,0 +1,120 @@
+# $Id$
+# Copyright W. Bangerth, University of Heidelberg, 1998
+
+# Template for makefiles for the examples subdirectory. In principle,
+# everything should be done automatically if you set the target file
+# here correctly:
+target   = step-4
+
+# 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
+
+# If you want your program to be linked with extra object or library
+# files, specify them here:
+user-libs =
+
+# To run the program, use "make run"; to give parameters to the program,
+# give the parameters to the following variable:
+run-parameters  = 
+
+# To execute additional action apart from running the program, fill
+# in this list:
+additional-run-action =
+
+# To specify which files are to be deleted by "make clean" (apart from
+# the usual ones: object files, executables, backups, etc), fill in the
+# following list
+delete-files = *gpl *inp *history
+
+
+
+
+###############################################################################
+# Internals
+
+#deal include base path
+D = ../../../..
+
+include $D/common/Make.global_options
+
+# get lists of files we need
+cc-files    = $(filter-out *%, $(shell echo *.cc))
+o-files     = $(cc-files:.cc=.o)
+go-files    = $(cc-files:.cc=.go)
+h-files     = $(filter-out *%, $(shell echo *.h))
+lib-h-files = $(filter-out *%, $(shell echo ../../include/*/*.h))
+
+# list of libraries needed to link with
+libs     = ./Obj.a   -ldeal_II_2d  -llac -lbase
+libs.g   = ./Obj.g.a -ldeal_II_3d.g -llac.g -lbase.g
+
+
+# check whether we use debug mode or not
+ifeq ($(debug-mode),on)
+libraries = $(libs.g)
+flags     = $(CXXFLAGS.g)
+endif
+
+ifeq ($(debug-mode),off)
+libraries = $(libs)
+flags     = $(CXXFLAGS)
+endif
+
+
+
+# make rule for the target
+$(target) : $(libraries)
+       @echo ============================ Linking $@
+       @$(CXX) $(flags) -o $@ $^ $(user-libs)
+
+# rule how to run the program
+run: $(target)
+       $(target) $(run-parameters)
+       $(additional-run-action)
+
+
+# rule to make object files
+%.go : %.cc
+       @echo ============================ Compiling with debugging information:   $<
+       @echo $(CXX) ... -c $< -o $@
+       @$(CXX) $(CXXFLAGS.g) -c $< -o $@
+%.o : %.cc
+       @echo ============================ Compiling with optimization:   $<
+       @echo $(CXX) ... -c $< -o $@
+       @$(CXX) $(CXXFLAGS) -c $< -o $@
+
+
+# rules which files the libraries depend upon
+Obj.a: ./Obj.a($(o-files))
+Obj.g.a: ./Obj.g.a($(go-files))
+
+
+clean:
+       -rm -f *.o *.go *~ Makefile.dep Obj.a Obj.g.a $(target) $(delete-files)
+
+
+
+.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.
+#
+#use perl to generate rules for the .go files as well
+#as to make rules not for tria.o and the like, but
+#rather for libnumerics.a(tria.o)
+Makefile.dep: $(cc-files) $(h-files) $(lib-h-files)
+       @echo ============================ Remaking Makefile
+       @perl $D/common/scripts/Make_dep.pl ./Obj $(INCLUDE) $(cc-files) \
+               > Makefile.dep
+
+
+include Makefile.dep
+
diff --git a/deal.II/deal.II/Attic/examples/step-by-step/step-4/step-4.cc b/deal.II/deal.II/Attic/examples/step-by-step/step-4/step-4.cc
new file mode 100644 (file)
index 0000000..901bfa6
--- /dev/null
@@ -0,0 +1,546 @@
+/* $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>
+
+                                // This is new, however: in the
+                                // previous example we got some
+                                // unwanted output from the linear
+                                // solvers. If we want to suppress
+                                // it, we have to include this file
+                                // and add a line somewhere to the
+                                // program; in this program, it was
+                                // added to the main function.
+#include <base/logstream.h>
+
+
+
+
+
+       
+
+
+
+template <int dim>
+class LaplaceProblem 
+{
+  public:
+    LaplaceProblem ();
+    void run ();
+    
+  private:
+    void make_grid_and_dofs ();
+    void assemble_system ();
+    void solve ();
+    void output_results ();
+
+    Triangulation<dim>   triangulation;
+    FEQ1<dim>            fe;
+    DoFHandler<dim>      dof_handler;
+
+    SparseMatrixStruct   sparsity_pattern;
+    SparseMatrix<double> system_matrix;
+
+    Vector<double>       solution;
+    Vector<double>       system_rhs;
+};
+
+
+                                // In the following, we declare two
+                                // more classes, which will represent
+                                // the functions of the
+                                // dim-dimensional space denoting the
+                                // right hand side and the
+                                // non-homogeneous Dirichlet boundary
+                                // values.
+                                //
+                                // Each of these classes is derived
+                                // from a common, abstract base class
+                                // Function, which declares the
+                                // common interface which all
+                                // functions have to follow. In
+                                // particular, concrete classes have
+                                // to overload the `value' function,
+                                // which takes a point in
+                                // dim-dimensional space as
+                                // parameters and shall return the
+                                // value at that point as a `double'
+                                // variable.
+template <int dim>
+class RightHandSide : public Function<dim> 
+{
+  public:
+    virtual double value (const Point<dim>   &p,
+                         const unsigned int  component) const;
+};
+
+
+
+template <int dim>
+class BoundaryValues : public Function<dim> 
+{
+  public:
+    virtual double value (const Point<dim>   &p,
+                         const unsigned int  component) const;
+};
+
+
+
+
+                                // We wanted the right hand side
+                                // function to be 4*(x**4+y**4) in
+                                // 2D, or 4*(x**4+y**4+z**4) in
+                                // 3D. Unfortunately, this is not as
+                                // elegantly feasible dimension
+                                // independently as much of the rest
+                                // of this program, so we have to do
+                                // it using a small
+                                // loop. Fortunately, the compiler
+                                // knows the size of the loop at
+                                // compile time, i.e. the number of
+                                // times the body will be executed,
+                                // so it can optimize away the
+                                // overhead needed for the loop and
+                                // the result will be as fast as if
+                                // we had used the formulas above
+                                // right away.
+                                //
+                                // Note that the different
+                                // coordinates of the point are
+                                // accessed using the () operator.
+template <int dim>
+double RightHandSide<dim>::value (const Point<dim> &p,
+                                 const unsigned int) const 
+{
+  double return_value = 0;
+  for (unsigned int i=0; i<dim; ++i)
+    return_value += 4*pow(p(i), 4);
+
+  return return_value;
+};
+
+
+                                // The boundary values were to be
+                                // chosen to be x*x+y*y in 2D, and
+                                // x*x+y*y+z*z in 3D. This happens to
+                                // be equal to the square of the
+                                // vector from the origin to the
+                                // point at which we would like to
+                                // evaluate the function,
+                                // irrespective of the dimension. So
+                                // that is what we return:
+template <int dim>
+double BoundaryValues<dim>::value (const Point<dim> &p,
+                                  const unsigned int) const 
+{
+  return p.square();
+};
+
+
+
+
+                                // This is the constructor of the
+                                // LaplaceProblem class. It
+                                // associates the DoFHandler to the
+                                // triangulation just as in the
+                                // previous example.
+template <int dim>
+LaplaceProblem<dim>::LaplaceProblem () :
+               dof_handler (triangulation)
+{};
+
+
+
+                                // Grid creation is something
+                                // inherently dimension
+                                // dependent. However, as long as the
+                                // domains are sufficiently similar
+                                // in 2D or 3D, the library can
+                                // abstract for you. In our case, we
+                                // would like to again solve on the
+                                // square [-1,1]x[-1,1] in 2D, or on
+                                // the cube [-1,1]x[-1,1]x[-1,1] in
+                                // 3D; both can be termed
+                                // ``hyper_cube'', so we may use the
+                                // same function in whatever
+                                // dimension we are. Of course, the
+                                // functions that create a hypercube
+                                // in two and three dimensions are
+                                // very much different, but that is
+                                // something you need not care
+                                // about. Let the library handle the
+                                // difficult things.
+                                //
+                                // Likewise, associating a degree of
+                                // freedom with each vertex is
+                                // something which certainly looks
+                                // different in 2D and 3D, but that
+                                // does not need to bother you. This
+                                // function therefore looks exactly
+                                // like in the previous example,
+                                // although it performs actions that
+                                // in their details are quite
+                                // different. The only significant
+                                // difference is the number of cells
+                                // resulting, which is much higher in
+                                // three than in two space
+                                // dimensions!
+template <int dim>
+void LaplaceProblem<dim>::make_grid_and_dofs ()
+{
+  GridGenerator::hyper_cube (triangulation, -1, 1);
+  triangulation.refine_global (4);
+  
+  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());
+};
+
+
+
+                                // Unlike in the previous example, we
+                                // would now like to use a
+                                // non-constant right hand side
+                                // function and non-zero boundary
+                                // values. Both are tasks that are
+                                // readily achieved with a only a few
+                                // new lines of code in the
+                                // assemblage of the matrix and right
+                                // hand side.
+                                //
+                                // More interesting, though, is they
+                                // way we assemble matrix and right
+                                // hand side vector dimension
+                                // independently: there is simply no
+                                // difference to the pure
+                                // two-dimensional case. Since the
+                                // important objects used in this
+                                // function (quadrature formula,
+                                // FEValues) depend on the dimension
+                                // by way of a template parameter as
+                                // well, they can take care of
+                                // setting up properly everything for
+                                // the dimension for which this
+                                // function is compiled. By declaring
+                                // all classes which might depend on
+                                // the dimension using a template
+                                // parameter, the library can make
+                                // nearly all work for you and you
+                                // don't have to care about most
+                                // things.
+template <int dim>
+void LaplaceProblem<dim>::assemble_system () 
+{  
+  QGauss3<dim>  quadrature_formula;
+
+                                  // We wanted to have a non-constant
+                                  // right hand side, so we use an
+                                  // object of the class declared
+                                  // above to generate the necessary
+                                  // data. Since this right hand side
+                                  // object is only used in this
+                                  // function, we only declare it
+                                  // here, rather than as a member
+                                  // variable of the LaplaceProblem
+                                  // class, or somewhere else.
+  const RightHandSide<dim> right_hand_side;
+
+                                  // Compared to the previous
+                                  // example, in order to evaluate
+                                  // the non-constant right hand side
+                                  // function we now also need the
+                                  // quadrature points on the cell we
+                                  // are presently on (previously,
+                                  // they were only needed on the
+                                  // unit cell, in order to compute
+                                  // the values and gradients of the
+                                  // shape function, which are
+                                  // defined on the unit cell
+                                  // however). We can tell the
+                                  // FEValues object to do for us by
+                                  // giving it the update_q_points
+                                  // flag:
+  FEValues<dim> fe_values (fe, quadrature_formula, 
+                          UpdateFlags(update_values    |
+                                      update_gradients |
+                                      update_q_points  |
+                                      update_JxW_values));
+
+                                  // Note that the following numbers
+                                  // depend on the dimension which we
+                                  // are presently using. However,
+                                  // the FE and Quadrature classes do
+                                  // all the necessary work for you
+                                  // and you don't have to care about
+                                  // the dimension dependent parts:
+  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);
+
+                                  // Note here, that a cell is a
+                                  // quadrilateral in two space
+                                  // dimensions, but a hexahedron in
+                                  // 3D. In fact, the
+                                  // active_cell_iterator data type
+                                  // is something different,
+                                  // depending on the dimension we
+                                  // are in, but to the outside world
+                                  // they look alike and you will
+                                  // probably never see a difference
+                                  // although they are totally
+                                  // unrelated.
+  DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
+                                       endc = dof_handler.end();
+  for (; cell!=endc; ++cell)
+    {
+      fe_values.reinit (cell);
+      cell_matrix.clear ();
+      cell_rhs.clear ();
+
+                                      // Now we have to assemble the
+                                      // local matrix and right hand
+                                      // side. This is done exactly
+                                      // like in the previous
+                                      // example, but now we revert
+                                      // the order of the loops
+                                      // (which we can safely do
+                                      // since they are independent
+                                      // of each other) and merge the
+                                      // loops for the local matrix
+                                      // and the local vector as far
+                                      // as possible; this makes
+                                      // things a bit faster.
+      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) += (fe_values.shape_grad (i, q_point) *
+                                  fe_values.shape_grad (j, q_point) *
+                                  fe_values.JxW (q_point));
+
+                                            // Here is about the only
+                                            // difference to the
+                                            // previous example:
+                                            // instead of using a
+                                            // constant right hand
+                                            // side, we use the
+                                            // respective object and
+                                            // evaluate it at the
+                                            // quadrature points.
+           cell_rhs(i) += (fe_values.shape_value (i, q_point) *
+                           right_hand_side.value (fe_values.quadrature_point (q_point)) *
+                           fe_values.JxW (q_point));
+         };
+
+
+                                      // The transfer into the global
+                                      // matrix and right hand side
+                                      // is done exactly as before,
+                                      // but here we have again
+                                      // merged some loops for
+                                      // efficiency:
+      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);
+       };
+    };
+
+  
+                                  // We wanted to have
+                                  // non-homogeneous boundary values
+                                  // in this example, contrary to the
+                                  // one before. This is a simple
+                                  // task, we only have to replace
+                                  // the ZeroFunction used there by
+                                  // an object of the class which
+                                  // describes the boundary values we
+                                  // would like to use (i.e. the
+                                  // BoundaryValues class declared
+                                  // above):
+  map<int,double> boundary_values;
+  VectorTools::interpolate_boundary_values (dof_handler,
+                                           0,
+                                           BoundaryValues<dim>(),
+                                           boundary_values);
+  MatrixTools<dim>::apply_boundary_values (boundary_values,
+                                          system_matrix,
+                                          solution,
+                                          system_rhs);
+};
+
+
+                                // Solving the linear system of
+                                // equation is something that looks
+                                // almost identical in most
+                                // programs. In particular, it is
+                                // dimension independent, so this
+                                // function is mostly copied from the
+                                // previous example.
+template <int dim>
+void LaplaceProblem<dim>::solve () 
+{
+  SolverControl           solver_control (1000, 1e-12);
+  PrimitiveVectorMemory<> vector_memory;
+  SolverCG<>              cg (solver_control, vector_memory);
+  cg.solve (system_matrix, solution, system_rhs,
+           PreconditionIdentity());
+
+                                  // We have made one addition,
+                                  // though: since we suppress output
+                                  // from the linear solvers, we have
+                                  // to print the number of
+                                  // iterations by hand.
+  cout << "  " << solver_control.last_step()
+       << " CG iterations needed to obtain convergence."
+       << endl;
+};
+
+
+
+                                // This function also does what the
+                                // respective one did in the previous
+                                // example. No changes here for
+                                // dimension independentce either.
+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 ();
+
+                                  // Only difference to the previous
+                                  // example: write output in GMV
+                                  // format, rather than for gnuplot.
+  ofstream output ("solution.gmv");
+  data_out.write_gmv (output);
+};
+
+
+
+                                // This is the function which has the
+                                // top-level control over
+                                // everything. Apart from one line of
+                                // additional output, it is the same
+                                // as for the previous example.
+template <int dim>
+void LaplaceProblem<dim>::run () 
+{
+  cout << "Solving problem in " << dim << " space dimensions." << endl;
+  
+  make_grid_and_dofs();
+  assemble_system ();
+  solve ();
+  output_results ();
+};
+
+    
+
+                                // And this is the main function. It
+                                // also looks mostly like in the
+                                // previous example:
+int main () 
+{
+                                  // In the previous example, we had
+                                  // the output from the linear
+                                  // solvers about the starting
+                                  // residual and the number of the
+                                  // iteration where convergence was
+                                  // detected. This can be suppressed
+                                  // like this:
+  deallog.depth_console (0);
+                                  // The rationale here is the
+                                  // following: the deallog
+                                  // (i.e. deal-log, not de-allog)
+                                  // variable represents a stream to
+                                  // which some parts of the library
+                                  // write output. It redirects this
+                                  // output to the console and if
+                                  // required to a file. The output
+                                  // is nested in a way that each
+                                  // function can use a prefix string
+                                  // (separated by colons) for each
+                                  // line of output; if it calls
+                                  // another function, that may also
+                                  // use its prefix which is then
+                                  // printed after the one of the
+                                  // calling function. Since output
+                                  // from functions which are nested
+                                  // deep below is usually not as
+                                  // important as top-level output,
+                                  // you can give the deallog
+                                  // variable a maximal depth of
+                                  // nested output for output to
+                                  // console and file. The depth zero
+                                  // which we gave here means that no
+                                  // output is written.
+
+                                  // After having done this
+                                  // administrative stuff, we can go
+                                  // on just as before: define one of
+                                  // these top-level objects and
+                                  // transfer control to it:
+//  LaplaceProblem<2> laplace_problem_2d;
+//  laplace_problem_2d.run ();
+
+    LaplaceProblem<3> laplace_problem_3d;
+    laplace_problem_3d.run ();
+  return 0;
+};
index 7307d29ed049b7100d7138e2dddf11dfa6b98212..fbb14abf656a2bf70a5f1d3b38464967ec6bc543 100644 (file)
          <a href="step-3.html" target="body">Step 3</a>
        </p>
       </li>
+
+      <li>
+       <p>
+         <a href="step-4.html" target="body">Step 4</a>
+       </p>
+      </li>
     </ol>
     <p>
       <a href="../index.html" target="_top">Back to the tutorial index</a>
diff --git a/deal.II/doc/tutorial/chapter-2.step-by-step/step-4.data/grid-3d.jpg b/deal.II/doc/tutorial/chapter-2.step-by-step/step-4.data/grid-3d.jpg
new file mode 100644 (file)
index 0000000..da8b54f
Binary files /dev/null and b/deal.II/doc/tutorial/chapter-2.step-by-step/step-4.data/grid-3d.jpg differ
diff --git a/deal.II/doc/tutorial/chapter-2.step-by-step/step-4.data/solution-3d.jpg b/deal.II/doc/tutorial/chapter-2.step-by-step/step-4.data/solution-3d.jpg
new file mode 100644 (file)
index 0000000..d7db56d
Binary files /dev/null and b/deal.II/doc/tutorial/chapter-2.step-by-step/step-4.data/solution-3d.jpg differ
diff --git a/deal.II/doc/tutorial/chapter-2.step-by-step/step-4.intro b/deal.II/doc/tutorial/chapter-2.step-by-step/step-4.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-4.results b/deal.II/doc/tutorial/chapter-2.step-by-step/step-4.results
new file mode 100644 (file)
index 0000000..4c30e00
--- /dev/null
@@ -0,0 +1,23 @@
+<a name="Results"></a>
+<h1>Results</h1>
+
+The output of the program looks as follows:
+<pre>
+<code>
+...
+</code>
+</pre>
+
+
+<p>
+<TABLE WIDTH="100%">
+<tr>
+<td>
+<IMG SRC="step-4.data/solution-3d.jpg" ALT="solution-3d" HEIGHT="300">
+</td>
+<td>
+<IMG SRC="step-4.data/grid-3d.jpg" ALT="grid-3d" HEIGHT="300">
+</td>
+</tr>
+</table>
+<p>
index 2e05de8a0d2496746a5531f3753dd21005bec708..42b75b9e5cd24f565902d7a7587c33eead7adda8 100644 (file)
@@ -57,6 +57,14 @@ At present, the following programs exist:
       problem. Object-orientation. Assembling matrices and
       vectors. Boundary values.
       </dd>
+
+      <dt><a href="step-4.html">Step 4</a></dt> 
+      <dd><strong>What's new:</strong> This example is programmed in a
+      way that it is independent of the dimension for which we want to
+      solve Laplace's equation; we will solve the equation in 2D and
+      3D, although the program is exactly the same. Non-constant right
+      hand side function. Non-homogeneous boundary values.
+      </dd>
 </dl>
 
 <!-- Page Foot -->
diff --git a/deal.II/examples/step-4/Makefile b/deal.II/examples/step-4/Makefile
new file mode 100644 (file)
index 0000000..9212069
--- /dev/null
@@ -0,0 +1,120 @@
+# $Id$
+# Copyright W. Bangerth, University of Heidelberg, 1998
+
+# Template for makefiles for the examples subdirectory. In principle,
+# everything should be done automatically if you set the target file
+# here correctly:
+target   = step-4
+
+# 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
+
+# If you want your program to be linked with extra object or library
+# files, specify them here:
+user-libs =
+
+# To run the program, use "make run"; to give parameters to the program,
+# give the parameters to the following variable:
+run-parameters  = 
+
+# To execute additional action apart from running the program, fill
+# in this list:
+additional-run-action =
+
+# To specify which files are to be deleted by "make clean" (apart from
+# the usual ones: object files, executables, backups, etc), fill in the
+# following list
+delete-files = *gpl *inp *history
+
+
+
+
+###############################################################################
+# Internals
+
+#deal include base path
+D = ../../../..
+
+include $D/common/Make.global_options
+
+# get lists of files we need
+cc-files    = $(filter-out *%, $(shell echo *.cc))
+o-files     = $(cc-files:.cc=.o)
+go-files    = $(cc-files:.cc=.go)
+h-files     = $(filter-out *%, $(shell echo *.h))
+lib-h-files = $(filter-out *%, $(shell echo ../../include/*/*.h))
+
+# list of libraries needed to link with
+libs     = ./Obj.a   -ldeal_II_2d  -llac -lbase
+libs.g   = ./Obj.g.a -ldeal_II_3d.g -llac.g -lbase.g
+
+
+# check whether we use debug mode or not
+ifeq ($(debug-mode),on)
+libraries = $(libs.g)
+flags     = $(CXXFLAGS.g)
+endif
+
+ifeq ($(debug-mode),off)
+libraries = $(libs)
+flags     = $(CXXFLAGS)
+endif
+
+
+
+# make rule for the target
+$(target) : $(libraries)
+       @echo ============================ Linking $@
+       @$(CXX) $(flags) -o $@ $^ $(user-libs)
+
+# rule how to run the program
+run: $(target)
+       $(target) $(run-parameters)
+       $(additional-run-action)
+
+
+# rule to make object files
+%.go : %.cc
+       @echo ============================ Compiling with debugging information:   $<
+       @echo $(CXX) ... -c $< -o $@
+       @$(CXX) $(CXXFLAGS.g) -c $< -o $@
+%.o : %.cc
+       @echo ============================ Compiling with optimization:   $<
+       @echo $(CXX) ... -c $< -o $@
+       @$(CXX) $(CXXFLAGS) -c $< -o $@
+
+
+# rules which files the libraries depend upon
+Obj.a: ./Obj.a($(o-files))
+Obj.g.a: ./Obj.g.a($(go-files))
+
+
+clean:
+       -rm -f *.o *.go *~ Makefile.dep Obj.a Obj.g.a $(target) $(delete-files)
+
+
+
+.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.
+#
+#use perl to generate rules for the .go files as well
+#as to make rules not for tria.o and the like, but
+#rather for libnumerics.a(tria.o)
+Makefile.dep: $(cc-files) $(h-files) $(lib-h-files)
+       @echo ============================ Remaking Makefile
+       @perl $D/common/scripts/Make_dep.pl ./Obj $(INCLUDE) $(cc-files) \
+               > Makefile.dep
+
+
+include Makefile.dep
+
diff --git a/deal.II/examples/step-4/step-4.cc b/deal.II/examples/step-4/step-4.cc
new file mode 100644 (file)
index 0000000..901bfa6
--- /dev/null
@@ -0,0 +1,546 @@
+/* $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>
+
+                                // This is new, however: in the
+                                // previous example we got some
+                                // unwanted output from the linear
+                                // solvers. If we want to suppress
+                                // it, we have to include this file
+                                // and add a line somewhere to the
+                                // program; in this program, it was
+                                // added to the main function.
+#include <base/logstream.h>
+
+
+
+
+
+       
+
+
+
+template <int dim>
+class LaplaceProblem 
+{
+  public:
+    LaplaceProblem ();
+    void run ();
+    
+  private:
+    void make_grid_and_dofs ();
+    void assemble_system ();
+    void solve ();
+    void output_results ();
+
+    Triangulation<dim>   triangulation;
+    FEQ1<dim>            fe;
+    DoFHandler<dim>      dof_handler;
+
+    SparseMatrixStruct   sparsity_pattern;
+    SparseMatrix<double> system_matrix;
+
+    Vector<double>       solution;
+    Vector<double>       system_rhs;
+};
+
+
+                                // In the following, we declare two
+                                // more classes, which will represent
+                                // the functions of the
+                                // dim-dimensional space denoting the
+                                // right hand side and the
+                                // non-homogeneous Dirichlet boundary
+                                // values.
+                                //
+                                // Each of these classes is derived
+                                // from a common, abstract base class
+                                // Function, which declares the
+                                // common interface which all
+                                // functions have to follow. In
+                                // particular, concrete classes have
+                                // to overload the `value' function,
+                                // which takes a point in
+                                // dim-dimensional space as
+                                // parameters and shall return the
+                                // value at that point as a `double'
+                                // variable.
+template <int dim>
+class RightHandSide : public Function<dim> 
+{
+  public:
+    virtual double value (const Point<dim>   &p,
+                         const unsigned int  component) const;
+};
+
+
+
+template <int dim>
+class BoundaryValues : public Function<dim> 
+{
+  public:
+    virtual double value (const Point<dim>   &p,
+                         const unsigned int  component) const;
+};
+
+
+
+
+                                // We wanted the right hand side
+                                // function to be 4*(x**4+y**4) in
+                                // 2D, or 4*(x**4+y**4+z**4) in
+                                // 3D. Unfortunately, this is not as
+                                // elegantly feasible dimension
+                                // independently as much of the rest
+                                // of this program, so we have to do
+                                // it using a small
+                                // loop. Fortunately, the compiler
+                                // knows the size of the loop at
+                                // compile time, i.e. the number of
+                                // times the body will be executed,
+                                // so it can optimize away the
+                                // overhead needed for the loop and
+                                // the result will be as fast as if
+                                // we had used the formulas above
+                                // right away.
+                                //
+                                // Note that the different
+                                // coordinates of the point are
+                                // accessed using the () operator.
+template <int dim>
+double RightHandSide<dim>::value (const Point<dim> &p,
+                                 const unsigned int) const 
+{
+  double return_value = 0;
+  for (unsigned int i=0; i<dim; ++i)
+    return_value += 4*pow(p(i), 4);
+
+  return return_value;
+};
+
+
+                                // The boundary values were to be
+                                // chosen to be x*x+y*y in 2D, and
+                                // x*x+y*y+z*z in 3D. This happens to
+                                // be equal to the square of the
+                                // vector from the origin to the
+                                // point at which we would like to
+                                // evaluate the function,
+                                // irrespective of the dimension. So
+                                // that is what we return:
+template <int dim>
+double BoundaryValues<dim>::value (const Point<dim> &p,
+                                  const unsigned int) const 
+{
+  return p.square();
+};
+
+
+
+
+                                // This is the constructor of the
+                                // LaplaceProblem class. It
+                                // associates the DoFHandler to the
+                                // triangulation just as in the
+                                // previous example.
+template <int dim>
+LaplaceProblem<dim>::LaplaceProblem () :
+               dof_handler (triangulation)
+{};
+
+
+
+                                // Grid creation is something
+                                // inherently dimension
+                                // dependent. However, as long as the
+                                // domains are sufficiently similar
+                                // in 2D or 3D, the library can
+                                // abstract for you. In our case, we
+                                // would like to again solve on the
+                                // square [-1,1]x[-1,1] in 2D, or on
+                                // the cube [-1,1]x[-1,1]x[-1,1] in
+                                // 3D; both can be termed
+                                // ``hyper_cube'', so we may use the
+                                // same function in whatever
+                                // dimension we are. Of course, the
+                                // functions that create a hypercube
+                                // in two and three dimensions are
+                                // very much different, but that is
+                                // something you need not care
+                                // about. Let the library handle the
+                                // difficult things.
+                                //
+                                // Likewise, associating a degree of
+                                // freedom with each vertex is
+                                // something which certainly looks
+                                // different in 2D and 3D, but that
+                                // does not need to bother you. This
+                                // function therefore looks exactly
+                                // like in the previous example,
+                                // although it performs actions that
+                                // in their details are quite
+                                // different. The only significant
+                                // difference is the number of cells
+                                // resulting, which is much higher in
+                                // three than in two space
+                                // dimensions!
+template <int dim>
+void LaplaceProblem<dim>::make_grid_and_dofs ()
+{
+  GridGenerator::hyper_cube (triangulation, -1, 1);
+  triangulation.refine_global (4);
+  
+  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());
+};
+
+
+
+                                // Unlike in the previous example, we
+                                // would now like to use a
+                                // non-constant right hand side
+                                // function and non-zero boundary
+                                // values. Both are tasks that are
+                                // readily achieved with a only a few
+                                // new lines of code in the
+                                // assemblage of the matrix and right
+                                // hand side.
+                                //
+                                // More interesting, though, is they
+                                // way we assemble matrix and right
+                                // hand side vector dimension
+                                // independently: there is simply no
+                                // difference to the pure
+                                // two-dimensional case. Since the
+                                // important objects used in this
+                                // function (quadrature formula,
+                                // FEValues) depend on the dimension
+                                // by way of a template parameter as
+                                // well, they can take care of
+                                // setting up properly everything for
+                                // the dimension for which this
+                                // function is compiled. By declaring
+                                // all classes which might depend on
+                                // the dimension using a template
+                                // parameter, the library can make
+                                // nearly all work for you and you
+                                // don't have to care about most
+                                // things.
+template <int dim>
+void LaplaceProblem<dim>::assemble_system () 
+{  
+  QGauss3<dim>  quadrature_formula;
+
+                                  // We wanted to have a non-constant
+                                  // right hand side, so we use an
+                                  // object of the class declared
+                                  // above to generate the necessary
+                                  // data. Since this right hand side
+                                  // object is only used in this
+                                  // function, we only declare it
+                                  // here, rather than as a member
+                                  // variable of the LaplaceProblem
+                                  // class, or somewhere else.
+  const RightHandSide<dim> right_hand_side;
+
+                                  // Compared to the previous
+                                  // example, in order to evaluate
+                                  // the non-constant right hand side
+                                  // function we now also need the
+                                  // quadrature points on the cell we
+                                  // are presently on (previously,
+                                  // they were only needed on the
+                                  // unit cell, in order to compute
+                                  // the values and gradients of the
+                                  // shape function, which are
+                                  // defined on the unit cell
+                                  // however). We can tell the
+                                  // FEValues object to do for us by
+                                  // giving it the update_q_points
+                                  // flag:
+  FEValues<dim> fe_values (fe, quadrature_formula, 
+                          UpdateFlags(update_values    |
+                                      update_gradients |
+                                      update_q_points  |
+                                      update_JxW_values));
+
+                                  // Note that the following numbers
+                                  // depend on the dimension which we
+                                  // are presently using. However,
+                                  // the FE and Quadrature classes do
+                                  // all the necessary work for you
+                                  // and you don't have to care about
+                                  // the dimension dependent parts:
+  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);
+
+                                  // Note here, that a cell is a
+                                  // quadrilateral in two space
+                                  // dimensions, but a hexahedron in
+                                  // 3D. In fact, the
+                                  // active_cell_iterator data type
+                                  // is something different,
+                                  // depending on the dimension we
+                                  // are in, but to the outside world
+                                  // they look alike and you will
+                                  // probably never see a difference
+                                  // although they are totally
+                                  // unrelated.
+  DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
+                                       endc = dof_handler.end();
+  for (; cell!=endc; ++cell)
+    {
+      fe_values.reinit (cell);
+      cell_matrix.clear ();
+      cell_rhs.clear ();
+
+                                      // Now we have to assemble the
+                                      // local matrix and right hand
+                                      // side. This is done exactly
+                                      // like in the previous
+                                      // example, but now we revert
+                                      // the order of the loops
+                                      // (which we can safely do
+                                      // since they are independent
+                                      // of each other) and merge the
+                                      // loops for the local matrix
+                                      // and the local vector as far
+                                      // as possible; this makes
+                                      // things a bit faster.
+      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) += (fe_values.shape_grad (i, q_point) *
+                                  fe_values.shape_grad (j, q_point) *
+                                  fe_values.JxW (q_point));
+
+                                            // Here is about the only
+                                            // difference to the
+                                            // previous example:
+                                            // instead of using a
+                                            // constant right hand
+                                            // side, we use the
+                                            // respective object and
+                                            // evaluate it at the
+                                            // quadrature points.
+           cell_rhs(i) += (fe_values.shape_value (i, q_point) *
+                           right_hand_side.value (fe_values.quadrature_point (q_point)) *
+                           fe_values.JxW (q_point));
+         };
+
+
+                                      // The transfer into the global
+                                      // matrix and right hand side
+                                      // is done exactly as before,
+                                      // but here we have again
+                                      // merged some loops for
+                                      // efficiency:
+      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);
+       };
+    };
+
+  
+                                  // We wanted to have
+                                  // non-homogeneous boundary values
+                                  // in this example, contrary to the
+                                  // one before. This is a simple
+                                  // task, we only have to replace
+                                  // the ZeroFunction used there by
+                                  // an object of the class which
+                                  // describes the boundary values we
+                                  // would like to use (i.e. the
+                                  // BoundaryValues class declared
+                                  // above):
+  map<int,double> boundary_values;
+  VectorTools::interpolate_boundary_values (dof_handler,
+                                           0,
+                                           BoundaryValues<dim>(),
+                                           boundary_values);
+  MatrixTools<dim>::apply_boundary_values (boundary_values,
+                                          system_matrix,
+                                          solution,
+                                          system_rhs);
+};
+
+
+                                // Solving the linear system of
+                                // equation is something that looks
+                                // almost identical in most
+                                // programs. In particular, it is
+                                // dimension independent, so this
+                                // function is mostly copied from the
+                                // previous example.
+template <int dim>
+void LaplaceProblem<dim>::solve () 
+{
+  SolverControl           solver_control (1000, 1e-12);
+  PrimitiveVectorMemory<> vector_memory;
+  SolverCG<>              cg (solver_control, vector_memory);
+  cg.solve (system_matrix, solution, system_rhs,
+           PreconditionIdentity());
+
+                                  // We have made one addition,
+                                  // though: since we suppress output
+                                  // from the linear solvers, we have
+                                  // to print the number of
+                                  // iterations by hand.
+  cout << "  " << solver_control.last_step()
+       << " CG iterations needed to obtain convergence."
+       << endl;
+};
+
+
+
+                                // This function also does what the
+                                // respective one did in the previous
+                                // example. No changes here for
+                                // dimension independentce either.
+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 ();
+
+                                  // Only difference to the previous
+                                  // example: write output in GMV
+                                  // format, rather than for gnuplot.
+  ofstream output ("solution.gmv");
+  data_out.write_gmv (output);
+};
+
+
+
+                                // This is the function which has the
+                                // top-level control over
+                                // everything. Apart from one line of
+                                // additional output, it is the same
+                                // as for the previous example.
+template <int dim>
+void LaplaceProblem<dim>::run () 
+{
+  cout << "Solving problem in " << dim << " space dimensions." << endl;
+  
+  make_grid_and_dofs();
+  assemble_system ();
+  solve ();
+  output_results ();
+};
+
+    
+
+                                // And this is the main function. It
+                                // also looks mostly like in the
+                                // previous example:
+int main () 
+{
+                                  // In the previous example, we had
+                                  // the output from the linear
+                                  // solvers about the starting
+                                  // residual and the number of the
+                                  // iteration where convergence was
+                                  // detected. This can be suppressed
+                                  // like this:
+  deallog.depth_console (0);
+                                  // The rationale here is the
+                                  // following: the deallog
+                                  // (i.e. deal-log, not de-allog)
+                                  // variable represents a stream to
+                                  // which some parts of the library
+                                  // write output. It redirects this
+                                  // output to the console and if
+                                  // required to a file. The output
+                                  // is nested in a way that each
+                                  // function can use a prefix string
+                                  // (separated by colons) for each
+                                  // line of output; if it calls
+                                  // another function, that may also
+                                  // use its prefix which is then
+                                  // printed after the one of the
+                                  // calling function. Since output
+                                  // from functions which are nested
+                                  // deep below is usually not as
+                                  // important as top-level output,
+                                  // you can give the deallog
+                                  // variable a maximal depth of
+                                  // nested output for output to
+                                  // console and file. The depth zero
+                                  // which we gave here means that no
+                                  // output is written.
+
+                                  // After having done this
+                                  // administrative stuff, we can go
+                                  // on just as before: define one of
+                                  // these top-level objects and
+                                  // transfer control to it:
+//  LaplaceProblem<2> laplace_problem_2d;
+//  laplace_problem_2d.run ();
+
+    LaplaceProblem<3> laplace_problem_3d;
+    laplace_problem_3d.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.