]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add yet another example program.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 14 Sep 2001 09:28:50 +0000 (09:28 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 14 Sep 2001 09:28:50 +0000 (09:28 +0000)
git-svn-id: https://svn.dealii.org/trunk@5001 0785d39b-7218-0410-832d-ea1e28bc413d

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

diff --git a/deal.II/examples/step-11/.cvsignore b/deal.II/examples/step-11/.cvsignore
new file mode 100644 (file)
index 0000000..1f752be
--- /dev/null
@@ -0,0 +1,3 @@
+*.o *.go Makefile.dep *.gnuplot *.gmv *.eps
+step-11
+*.dat
diff --git a/deal.II/examples/step-11/Makefile b/deal.II/examples/step-11/Makefile
new file mode 100644 (file)
index 0000000..c02ef0f
--- /dev/null
@@ -0,0 +1,159 @@
+# $Id$
+
+
+# For the small projects Makefile, you basically need to fill in only
+# four fields.
+#
+# The first is the name of the application. It is assumed that the
+# application name is the same as the base file name of the single C++
+# file from which the application is generated.
+target = $(basename $(shell echo step-*.cc))
+
+# The second field determines whether you want to run your program in
+# debug or optimized mode. The latter is significantly faster, but no
+# run-time checking of parameters and internal states is performed, so
+# you should set this value to `on' while you develop your program,
+# and to `off' when running production computations.
+debug-mode = on
+
+
+# As third field, we need to give the path to the top-level deal.II
+# directory. You need to adjust this to your needs. Since this path is
+# probably the most often needed one in the Makefile internals, it is
+# designated by a single-character variable, since that can be
+# reference using $D only, i.e. without the parentheses that are
+# required for most other parameters, as e.g. in $(target).
+D = ../../
+
+
+# The last field specifies the names of data and other files that
+# shall be deleted when calling `make clean'. Object and backup files,
+# executables and the like are removed anyway. Here, we give a list of
+# files in the various output formats that deal.II supports.
+clean-up-files = *gmv *gnuplot *gpl *eps *pov
+
+
+
+
+#
+#
+# Usually, you will not need to change something beyond this point.
+#
+#
+# The next statement tell the `make' program where to find the
+# deal.II top level directory and to include the file with the global
+# settings
+include $D/common/Make.global_options
+
+
+# Since the whole project consists of only one file, we need not
+# consider difficult dependencies. We only have to declare the
+# libraries which we want to link to the object file, and there need
+# to be two sets of libraries: one for the debug mode version of the
+# application and one for the optimized mode. Here we have selected
+# the versions for 2d. Note that the order in which the libraries are
+# given here is important and that your applications won't link
+# properly if they are given in another order.
+#
+# You may need to augment the lists of libraries when compiling your
+# program for other dimensions, or when using third party libraries
+libs.g   = $(lib-deal2-2d.g) \
+          $(lib-lac.g)      \
+           $(lib-base.g)
+libs.o   = $(lib-deal2-2d.o) \
+          $(lib-lac.o)      \
+           $(lib-base.o)
+
+
+# We now use the variable defined above which switch between debug and
+# optimized mode to select the correct compiler flags and the set of
+# libraries to link with. Included in the list of libraries is the
+# name of the object file which we will produce from the single C++
+# file. Note that by default we use the extension .go for object files
+# compiled in debug mode and .o for object files in optimized mode.
+ifeq ($(debug-mode),on)
+  libraries = $(target).go $(libs.g)
+  flags     = $(CXXFLAGS.g)
+else
+  libraries = $(target).o $(libs.o)
+  flags     = $(CXXFLAGS.o)
+endif
+
+
+# Now comes the first production rule: how to link the single object
+# file produced from the single C++ file into the executable. Since
+# this is the first rule in the Makefile, it is the one `make' selects
+# if you call it without arguments.
+$(target) : $(libraries)
+       @echo ============================ Linking $@
+       @$(CXX) -o $@ $^ $(LIBS) $(LDFLAGS)
+
+
+# To make running the application somewhat independent of the actual
+# program name, we usually declare a rule `run' which simply runs the
+# program. You can then run it by typing `make run'. This is also
+# useful if you want to call the executable with arguments which do
+# not change frequently. You may then want to add them to the
+# following rule:
+run: $(target)
+       @echo ============================ Running $<
+       @./$(target)
+
+
+# As a last rule to the `make' program, we define what to do when
+# cleaning up a directory. This usually involves deleting object files
+# and other automatically created files such as the executable itself,
+# backup files, and data files. Since the latter are not usually quite
+# diverse, you needed to declare them at the top of this file.
+clean:
+       -rm -f *.o *.go *~ Makefile.dep $(target) $(clean-up-files)
+
+
+# Since we have not yet stated how to make an object file from a C++
+# file, we should do so now. Since the many flags passed to the
+# compiler are usually not of much interest, we suppress the actual
+# command line using the `at' sign in the first column of the rules
+# and write the string indicating what we do instead.
+%.go : %.cc
+       @echo ==============debug========= $(<F)
+       @$(CXX) $(CXXFLAGS.g) -c $< -o $@
+%.o : %.cc
+       @echo ==============optimized===== $(<F)
+       @$(CXX) $(CXXFLAGS.o) -c $< -o $@
+
+
+# The following statement tells make that the rules `run' and `clean'
+# are not expected to produce files of the same name as Makefile rules
+# usually do.
+.PHONY: run clean
+
+
+# Finally there is a rule which you normally need not care much about:
+# since the executable depends on some include files from the library,
+# besides the C++ application file of course, it is necessary to
+# re-generate the executable when one of the files it depends on has
+# changed. The following rule to created a dependency file
+# `Makefile.dep', which `make' uses to determine when to regenerate
+# the executable. This file is automagically remade whenever needed,
+# i.e. whenever one of the cc-(include-path-base)/baseh-files changed. Make detects whether
+# to remake this file upon inclusion at the bottom of this file.
+#
+# The dependency file is created using a perl script.  Since the
+# script prefixes the output names by `lib(include-path-base)/baseo' or `lib(include-path-base)/basego' (it was
+# written for the sublibraries' Makefile), we have to strip that again
+# since object files are placed in the present directory for this
+# application. All these things are made in the next rule:
+Makefile.dep: $(target).cc Makefile \
+              $(shell echo $(include-path-base)/base/*.h    \
+                           $(include-path-lac)/lac/*.h      \
+                           $(include-path-deal2)/*/*.h)
+       @echo ============================ Remaking Makefile
+       @perl $D/common/scripts/make_dependencies.pl  $(INCLUDE) $(target).cc \
+               | perl -pi -e 's!lib/g?o/!!g;' \
+               > Makefile.dep
+
+# To make the dependencies known to `make', we finally have to include
+# them:
+include Makefile.dep
+
+
diff --git a/deal.II/examples/step-11/step-11.cc b/deal.II/examples/step-11/step-11.cc
new file mode 100644 (file)
index 0000000..a4f38dd
--- /dev/null
@@ -0,0 +1,400 @@
+/* $Id$ */
+/* Author: Wolfgang Bangerth, University of Heidelberg, 2001 */
+
+                                // As usual, the program starts with
+                                // a rather long list of include
+                                // files which you are probably
+                                // already used to by now:
+#include <base/quadrature_lib.h>
+#include <base/function.h>
+#include <base/logstream.h>
+#include <lac/vector.h>
+#include <lac/sparse_matrix.h>
+#include <lac/solver_cg.h>
+#include <lac/vector_memory.h>
+#include <lac/precondition.h>
+#include <grid/tria.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_boundary_lib.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <dofs/dof_handler.h>
+#include <dofs/dof_constraints.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+#include <fe/fe_q.h>
+#include <fe/fe_values.h>
+#include <fe/mapping_q.h>
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+
+                                // Just this one is new: it declares
+                                // a class
+                                // ``CompressedSparsityPattern'',
+                                // which we will use and explain
+                                // further down below.
+#include <lac/compressed_sparsity_pattern.h>
+
+                                // We will make use of the std::find
+                                // algorithm of the C++ standard
+                                // library, so we have to include the
+                                // following file for its
+                                // declaration:
+#include <algorithm>
+
+
+template <int dim>
+double measure (const DoFHandler<dim> &dof_handler,
+               const Mapping<dim>    &mapping)
+{
+  QGauss4<dim> quadrature;
+  FEValues<dim> fe_values (mapping, dof_handler.get_fe(), quadrature,
+                          update_JxW_values);
+  
+  typename DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+  double measure = 0;
+  for (; cell!=endc; ++cell)
+    {
+      fe_values.reinit (cell);
+      for (unsigned int i=0; i<fe_values.n_quadrature_points; ++i)
+       measure += fe_values.JxW (i);
+    };
+  return measure;
+};
+
+
+template <int dim>
+double measure (const Triangulation<dim> &triangulation,
+               const Mapping<dim>       &mapping)
+{
+  FE_Q<dim> dummy_fe(1);
+  DoFHandler<dim> dof_handler (const_cast<Triangulation<dim>&>(triangulation));
+  dof_handler.distribute_dofs(dummy_fe);
+  return measure (dof_handler, mapping);
+};
+
+
+template <int dim>
+double surface (const DoFHandler<dim> &dof_handler,
+               const Mapping<dim>    &mapping)
+{
+  QGauss4<dim-1> quadrature;
+  FEFaceValues<dim> fe_values (mapping, dof_handler.get_fe(), quadrature,
+                              update_JxW_values);
+  
+  typename DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+  double surface = 0;
+  for (; cell!=endc; ++cell)
+    for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+      if (cell->face(face)->at_boundary())
+       {
+         fe_values.reinit (cell, face);
+         for (unsigned int i=0; i<fe_values.n_quadrature_points; ++i)
+           surface += fe_values.JxW (i);
+       };
+  return surface;
+};
+
+
+template <int dim>
+double surface (const Triangulation<dim> &triangulation,
+               const Mapping<dim>       &mapping)
+{
+  FE_Q<dim> dummy_fe(1);
+  DoFHandler<dim> dof_handler (const_cast<Triangulation<dim>&>(triangulation));
+  dof_handler.distribute_dofs(dummy_fe);
+  return surface (dof_handler, mapping);
+};
+
+
+template double surface (const Triangulation<2> &, const Mapping<2> &);
+template double measure (const Triangulation<2> &, const Mapping<2> &);
+
+
+
+
+                                // Then we declare a class which
+                                // represents the solution of a
+                                // Laplace problem. As this example
+                                // program is based on step-5, the
+                                // class looks rather the same, with
+                                // the sole structural difference
+                                // that we have merged the functions
+                                // ``assemble_system'' and ``solve'',
+                                // and the output function was
+                                // dropped since the solution
+                                // function is so boring that it is
+                                // not worth being viewed.
+                                //
+                                // The only other noteworthy change
+                                // is that the constructor takes a
+                                // value representing the polynomial
+                                // degree of the mapping to be used
+                                // later on, and that it has another
+                                // member variable representing
+                                // exactly this mapping. In general,
+                                // this variable will occur in real
+                                // applications at the same places
+                                // where the finite element is
+                                // declared or used.
+template <int dim>
+class LaplaceProblem 
+{
+  public:
+    LaplaceProblem (const unsigned int mapping_degree);
+    void run ();
+    
+  private:
+    void setup_system ();
+    void assemble_and_solve ();
+    void solve ();
+
+    Triangulation<dim>   triangulation;
+    FE_Q<dim>            fe;
+    DoFHandler<dim>      dof_handler;
+    MappingQ<dim>        mapping;
+
+    SparsityPattern      sparsity_pattern;
+    SparseMatrix<double> system_matrix;
+    ConstraintMatrix     mean_value_constraints;
+
+    Vector<double>       solution;
+    Vector<double>       system_rhs;
+};
+
+
+
+                                // Construct such an object, by
+                                // initializing the variables. Here,
+                                // we use linear finite elements (the
+                                // argument to the ``fe'' variable
+                                // denotes the polynomial degree),
+                                // and mappings of given order. Print
+                                // to screen what we are about to do.
+template <int dim>
+LaplaceProblem<dim>::LaplaceProblem (const unsigned int mapping_degree) :
+                fe (1),
+               dof_handler (triangulation),
+               mapping (mapping_degree)
+{
+  std::cout << "Using mapping with degree " << mapping_degree << ":"
+           << std::endl
+           << "============================"
+           << std::endl;
+};
+
+
+
+                                // The first task is to set up the
+                                // variables for this problem. This
+                                // includes generating a valid
+                                // ``DoFHandler'' object, as well as
+                                // the sparsity patterns for the
+                                // matrix, and the object
+                                // representing the constraints that
+                                // the mean value of the degrees of
+                                // freedom on the boundary be zero.
+template <int dim>
+void LaplaceProblem<dim>::setup_system ()
+{
+                                  // The first task is trivial:
+                                  // generate an enumeration of the
+                                  // degrees of freedom:
+  dof_handler.distribute_dofs (fe);
+
+                                  // Next task is to construct the
+                                  // object representing the
+                                  // constraint that the mean value
+                                  // of the degrees of freedom on the
+                                  // boundary shall be zero. For
+                                  // this, we first want a list of
+                                  // those nodes which are actually
+                                  // at the boundary. The
+                                  // ``DoFTools'' class has a
+                                  // function that returns an array
+                                  // of boolean values where ``true''
+                                  // indicates that the node is at
+                                  // the boundary. The second
+                                  // argument denotes a mask
+                                  // selecting which components of
+                                  // vector valued finite elements we
+                                  // want to be considered. Since we
+                                  // have a scalar finite element
+                                  // anyway, this mask consists of
+                                  // only one entry, and its value
+                                  // must be ``true''.
+  std::vector<bool> boundary_dofs (dof_handler.n_dofs(), false);
+  DoFTools::extract_boundary_dofs (dof_handler, std::vector<bool>(1,true),
+                                  boundary_dofs);
+  
+                                  // Let us first pick out the first
+                                  // boundary node from this list. We
+                                  // do that by searching for the
+                                  // first ``true'' value in the
+                                  // array (note that ``std::find''
+                                  // returns an iterator to this
+                                  // element), and computing its
+                                  // distance to the overall first
+                                  // element in the array to get its
+                                  // index:
+  const unsigned int first_boundary_dof
+    = std::distance (std::find (boundary_dofs.begin(),
+                               boundary_dofs.end(),
+                               true),
+                    boundary_dofs.begin());
+       
+  mean_value_constraints.clear ();
+  mean_value_constraints.add_line (first_boundary_dof);
+  for (unsigned int i=first_boundary_dof+1; i<dof_handler.n_dofs(); ++i)
+    if (boundary_dofs[i] == true)
+      mean_value_constraints.add_entry (first_boundary_dof,
+                                       i, -1);
+  mean_value_constraints.close ();
+
+  CompressedSparsityPattern csp (dof_handler.n_dofs(),
+                                dof_handler.n_dofs());
+  DoFTools::make_sparsity_pattern (dof_handler, csp);
+  mean_value_constraints.condense (csp);
+
+  sparsity_pattern.copy_from (csp);
+
+  system_matrix.reinit (sparsity_pattern);
+
+  solution.reinit (dof_handler.n_dofs());
+  system_rhs.reinit (dof_handler.n_dofs());
+};
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::assemble_and_solve () 
+{  
+  QGauss2<dim>  cell_quadrature;
+  QGauss2<dim-1> face_quadrature;
+  MatrixTools::create_laplace_matrix (mapping, dof_handler,
+                                     cell_quadrature,
+                                     system_matrix);
+  VectorTools::create_right_hand_side (mapping, dof_handler,
+                                      cell_quadrature,
+                                      ConstantFunction<dim>(-2),
+                                      system_rhs);
+  
+  Vector<double> tmp (system_rhs.size());
+  VectorTools::create_boundary_right_hand_side (mapping, dof_handler,
+                                               face_quadrature,
+                                               ConstantFunction<dim>(1),
+                                               tmp);
+  system_rhs += tmp;
+
+  mean_value_constraints.condense (system_matrix);
+  mean_value_constraints.condense (system_rhs);  
+
+  solve ();
+  mean_value_constraints.distribute (solution);
+  
+  Vector<float> difference_per_cell (triangulation.n_active_cells());
+  VectorTools::integrate_difference (mapping, dof_handler,
+                                    solution,
+                                    ZeroFunction<dim>(),
+                                    difference_per_cell,
+                                    QGauss3<dim>(),
+                                    H1_seminorm);
+  std::cout << "  " << triangulation.n_active_cells() << " cells:  "
+           << "  |u|_1="
+           << difference_per_cell.l2_norm()
+           << ", error="
+           << fabs(difference_per_cell.l2_norm()-sqrt(3.14159265358/2))
+           << std::endl;
+};
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::solve () 
+{
+  SolverControl           solver_control (1000, 1e-12);
+  PrimitiveVectorMemory<> vector_memory;
+  SolverCG<>              cg (solver_control, vector_memory);
+
+  PreconditionSSOR<> preconditioner;
+  preconditioner.initialize(system_matrix, 1.2);
+
+  cg.solve (system_matrix, solution, system_rhs,
+           preconditioner);
+};
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::run () 
+{
+  GridGenerator::hyper_ball (triangulation);
+  static const HyperBallBoundary<dim> boundary;
+  triangulation.set_boundary (0, boundary);
+  
+  for (unsigned int cycle=0; cycle<6; ++cycle, triangulation.refine_global(1))
+    {
+      setup_system ();
+      assemble_and_solve ();
+    };
+};
+
+    
+
+                                // Finally the main function. It's
+                                // structure is the same as that used
+                                // in several of the previous
+                                // examples, so probably needs no
+                                // more explanation.
+int main () 
+{
+  try
+    {
+      deallog.depth_console (0);
+      std::cout.precision(5);
+
+                                      // This is the main loop, doing
+                                      // the computations with
+                                      // mappings of linear through
+                                      // cubic mappings. Note that
+                                      // since we need the object of
+                                      // type ``LaplaceProblem<2>''
+                                      // only once, we do not even
+                                      // name it, but create an
+                                      // unnamed such object and call
+                                      // the ``run'' function of it,
+                                      // subsequent to which it is
+                                      // immediately destroyed again.
+      for (unsigned int mapping_degree=1; mapping_degree<=3; ++mapping_degree)
+       LaplaceProblem<2>(mapping_degree).run ();
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+               << exc.what() << std::endl
+               << "Aborting!" << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      return 1;
+    }
+  catch (...) 
+    {
+      std::cerr << std::endl << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+               << "Aborting!" << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      return 1;
+    };
+
+  return 0;
+};

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.