]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add step-9, which is by no means complete, however.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 10 May 2000 13:45:15 +0000 (13:45 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 10 May 2000 13:45:15 +0000 (13:45 +0000)
git-svn-id: https://svn.dealii.org/trunk@2832 0785d39b-7218-0410-832d-ea1e28bc413d

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

diff --git a/deal.II/examples/step-9/Makefile b/deal.II/examples/step-9/Makefile
new file mode 100644 (file)
index 0000000..e24b954
--- /dev/null
@@ -0,0 +1,167 @@
+# $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
+
+
+# If in multithread mode, add the ACE library to the libraries which
+# we need to link with:
+ifneq ($(with-multithreading),no)
+  libraries += $(lib-ACE)
+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) $(flags) -o $@ $^
+
+
+# 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-9/step-9.cc b/deal.II/examples/step-9/step-9.cc
new file mode 100644 (file)
index 0000000..7de8c2e
--- /dev/null
@@ -0,0 +1,533 @@
+/* $Id$ */
+/* Author: Wolfgang Bangerth, University of Heidelberg, 2000 */
+
+#include <base/quadrature_lib.h>
+#include <base/function.h>
+#include <base/logstream.h>
+#include <lac/vector.h>
+#include <lac/full_matrix.h>
+#include <lac/sparse_matrix.h>
+#include <lac/solver_bicgstab.h>
+#include <lac/vector_memory.h>
+#include <lac/precondition.h>
+#include <grid/tria.h>
+#include <dofs/dof_handler.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_boundary_lib.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+#include <fe/fe_values.h>
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <numerics/data_out.h>
+#include <fe/fe_lib.lagrange.h>
+#include <grid/grid_out.h>
+
+#include <dofs/dof_constraints.h>
+
+#include <numerics/error_estimator.h>
+
+#include <fstream>
+
+
+template <int dim>
+class AdvectionProblem 
+{
+  public:
+    AdvectionProblem ();
+    ~AdvectionProblem ();
+    void run ();
+    
+  private:
+    void setup_system ();
+    void assemble_system ();
+    void solve ();
+    void refine_grid ();
+    void output_results (const unsigned int cycle) const;
+
+    Triangulation<dim>   triangulation;
+    DoFHandler<dim>      dof_handler;
+
+    FEQ1<dim>            fe;
+
+    ConstraintMatrix     hanging_node_constraints;
+
+    SparsityPattern      sparsity_pattern;
+    SparseMatrix<double> system_matrix;
+
+    Vector<double>       solution;
+    Vector<double>       system_rhs;
+};
+
+
+
+template <int dim>
+class AdvectionField 
+{
+  public:
+    Point<dim> value (const Point<dim> &p) const;
+    
+    void value_list (const vector<Point<dim> > &points,
+                    vector<Point<dim> >       &values) const;
+
+                                    /**
+                                     * Exception
+                                     */
+    DeclException2 (ExcVectorHasWrongSize,
+                   int, int,
+                   << "The vector has size " << arg1 << " but should have "
+                   << arg2 << " elements.");
+};
+
+
+
+template <int dim>
+Point<dim>
+AdvectionField<dim>::value (const Point<dim> &p) const 
+{
+  Point<dim> value;
+  value[0] = 2;
+  for (unsigned int i=1; i<dim; ++i)
+    value[i] = 1+0.8*sin(8*M_PI*p[0]);
+
+  return value;
+};
+
+
+
+template <int dim>
+void
+AdvectionField<dim>::value_list (const vector<Point<dim> > &points,
+                                vector<Point<dim> >       &values) const 
+{
+  Assert (values.size() == points.size(),
+         ExcVectorHasWrongSize (values.size(), points.size()));
+  
+  for (unsigned int i=0; i<points.size(); ++i)
+    values[i] = AdvectionField<dim>::value (points[i]);
+};
+
+
+
+
+template <int dim>
+class RightHandSide : 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;
+    
+  private:
+    static const Point<dim> center_point;
+};
+
+
+template <>
+const Point<1> RightHandSide<1>::center_point = Point<1> (-0.75);
+
+template <>
+const Point<2> RightHandSide<2>::center_point = Point<2> (-0.75, -0.75);
+
+template <>
+const Point<3> RightHandSide<3>::center_point = Point<3> (-0.75, -0.75, -0.75);
+
+
+
+template <int dim>
+double
+RightHandSide<dim>::value (const Point<dim>   &p,
+                          const unsigned int  component) const 
+{
+  Assert (component == 0, ExcWrongComponent (component, 1));
+  const double diameter = 0.1;
+  return ( (p-center_point).square() < diameter*diameter ?
+          .1/pow(diameter,dim) :
+          0);
+};
+
+
+
+template <int dim>
+void
+RightHandSide<dim>::value_list (const vector<Point<dim> > &points,
+                               vector<double>            &values,
+                               const unsigned int         component) const 
+{
+  Assert (values.size() == points.size(),
+         ExcVectorHasWrongSize (values.size(), points.size()));
+  
+  for (unsigned int i=0; i<points.size(); ++i)
+    values[i] = RightHandSide<dim>::value (points[i], component);
+};
+
+
+
+template <int dim>
+class BoundaryValues : 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
+BoundaryValues<dim>::value (const Point<dim>   &p,
+                           const unsigned int  component) const 
+{
+  Assert (component == 0, ExcWrongComponent (component, 1));
+
+  const double sine_term = sin(16*M_PI*sqrt(p.square()));
+  const double weight    = exp(-5*p.square()) / exp(-5.);
+  return sine_term * weight;
+};
+
+
+
+template <int dim>
+void
+BoundaryValues<dim>::value_list (const vector<Point<dim> > &points,
+                                vector<double>            &values,
+                                const unsigned int         component) const 
+{
+  Assert (values.size() == points.size(),
+         ExcVectorHasWrongSize (values.size(), points.size()));
+  
+  for (unsigned int i=0; i<points.size(); ++i)
+    values[i] = BoundaryValues<dim>::value (points[i], component);
+};
+
+
+
+template <int dim>
+AdvectionProblem<dim>::AdvectionProblem () :
+               dof_handler (triangulation)
+{};
+
+
+
+template <int dim>
+AdvectionProblem<dim>::~AdvectionProblem () 
+{
+  dof_handler.clear ();
+};
+
+
+
+template <int dim>
+void AdvectionProblem<dim>::setup_system ()
+{
+  dof_handler.distribute_dofs (fe);
+
+  hanging_node_constraints.clear ();
+  DoFTools::make_hanging_node_constraints (dof_handler,
+                                          hanging_node_constraints);
+  hanging_node_constraints.close ();
+
+  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);
+
+  hanging_node_constraints.condense (sparsity_pattern);
+
+  sparsity_pattern.compress();
+
+  system_matrix.reinit (sparsity_pattern);
+
+  solution.reinit (dof_handler.n_dofs());
+  system_rhs.reinit (dof_handler.n_dofs());
+};
+
+
+
+template <int dim>
+void AdvectionProblem<dim>::assemble_system () 
+{
+  AdvectionField<dim> advection_field;
+  RightHandSide<dim>  right_hand_side;
+  BoundaryValues<dim> boundary_values;
+  
+  QGauss3<dim>   quadrature_formula;
+  QGauss3<dim-1> face_quadrature_formula;
+  
+  FEValues<dim> fe_values (fe, quadrature_formula, 
+                          UpdateFlags(update_values    |
+                                      update_gradients |
+                                      update_q_points  |
+                                      update_JxW_values));
+  FEFaceValues<dim> fe_face_values (fe, face_quadrature_formula,
+                                   UpdateFlags (update_values     |
+                                                update_q_points   |
+                                                update_JxW_values |
+                                                update_normal_vectors));
+
+  const unsigned int   dofs_per_cell   = fe.dofs_per_cell;
+  const unsigned int   n_q_points      = quadrature_formula.n_quadrature_points;
+  const unsigned int   n_face_q_points = face_quadrature_formula.n_quadrature_points;
+
+  FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
+  Vector<double>       cell_rhs (dofs_per_cell);
+
+  vector<unsigned int> local_dof_indices (dofs_per_cell);
+
+  vector<double>       rhs_values (n_q_points);
+  vector<Point<dim> >  advection_directions (n_q_points);
+  vector<double>       face_boundary_values (n_face_q_points);
+  vector<Point<dim> >  face_advection_directions (n_face_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 ();
+
+      fe_values.reinit (cell);
+      const FullMatrix<double> 
+       & shape_values = fe_values.get_shape_values();
+      const vector<vector<Tensor<1,dim> > >
+       & shape_grads  = fe_values.get_shape_grads();
+      const vector<double>
+       & JxW_values   = fe_values.get_JxW_values();
+      const vector<Point<dim> >
+       & q_points     = fe_values.get_quadrature_points();
+
+      advection_field.value_list (q_points, advection_directions);
+      right_hand_side.value_list (q_points, rhs_values);
+
+      const double delta = 0.1 * cell->diameter ();
+
+      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) += ((advection_directions[q_point] *
+                                   shape_grads[j][q_point]    *
+                                   (shape_values(i,q_point) +
+                                    delta *
+                                    advection_directions[q_point] *
+                                    shape_grads[i][q_point])) *
+                                  JxW_values[q_point]);
+
+           cell_rhs(i) += ((shape_values (i,q_point) +
+                            delta *
+                            advection_directions[q_point] *
+                            shape_grads[i][q_point]        ) *
+                           rhs_values[i] *
+                           fe_values.JxW (q_point));
+         };
+
+      for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+       {
+         fe_face_values.reinit (cell, face);
+         
+         const FullMatrix<double> 
+           & face_shape_values = fe_face_values.get_shape_values();
+         const vector<double>
+           & face_JxW_values   = fe_face_values.get_JxW_values();
+         const vector<Point<dim> >
+           & face_q_points     = fe_face_values.get_quadrature_points();
+         const vector<Point<dim> >
+           & normal_vectors    = fe_face_values.get_normal_vectors();
+
+         boundary_values.value_list (face_q_points, face_boundary_values);
+         advection_field.value_list (face_q_points, face_advection_directions);
+
+         for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
+           if (cell->face(face)->at_boundary () &&
+               (normal_vectors[q_point] * face_advection_directions[q_point] < 0))
+             for (unsigned int i=0; i<dofs_per_cell; ++i)
+               {
+                 for (unsigned int j=0; j<dofs_per_cell; ++j)
+                   cell_matrix(i,j) -= (face_advection_directions[q_point] *
+                                        normal_vectors[q_point] *
+                                        face_shape_values(i,q_point) *
+                                        face_shape_values(j,q_point) *
+                                        face_JxW_values[q_point]);
+                 
+                 cell_rhs(i) -= (face_advection_directions[q_point] *
+                                 normal_vectors[q_point] *
+                                 face_boundary_values[q_point] *
+                                 face_shape_values(i,q_point) *
+                                 face_JxW_values[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);
+       };
+    };
+
+  hanging_node_constraints.condense (system_matrix);
+  hanging_node_constraints.condense (system_rhs);
+                                  // no bdr val
+};
+
+
+
+template <int dim>
+void AdvectionProblem<dim>::solve () 
+{
+  SolverControl           solver_control (1000, 1e-12);
+  PrimitiveVectorMemory<> vector_memory;
+  SolverBicgstab<>        bicgstab (solver_control, vector_memory);
+
+  PreconditionRelaxation<>
+    preconditioner(system_matrix,
+                  &SparseMatrix<double>::template precondition_Jacobi<double>,
+                  1.0);
+
+  bicgstab.solve (system_matrix, solution, system_rhs,
+                 preconditioner);
+
+  hanging_node_constraints.distribute (solution);
+};
+
+
+template <int dim>
+void AdvectionProblem<dim>::refine_grid ()
+{
+  Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+
+  KellyErrorEstimator<dim>::FunctionMap neumann_boundary;
+
+  KellyErrorEstimator<dim>::estimate (dof_handler,
+                                     QGauss3<dim-1>(),
+                                     neumann_boundary,
+                                     solution,
+                                     estimated_error_per_cell);
+
+  triangulation.refine_and_coarsen_fixed_number (estimated_error_per_cell,
+                                                0.3, 0.03);
+
+  triangulation.execute_coarsening_and_refinement ();
+};
+
+
+
+template <int dim>
+void AdvectionProblem<dim>::output_results (const unsigned int cycle) const
+{
+  string filename = "grid-";
+  filename += ('0' + cycle);
+  Assert (cycle < 10, ExcInternalError());
+  
+  filename += ".eps";
+  ofstream output (filename.c_str());
+
+  GridOut grid_out;
+  grid_out.write_eps (triangulation, output);
+};
+
+
+
+template <int dim>
+void AdvectionProblem<dim>::run () 
+{
+  for (unsigned int cycle=0; cycle<6; ++cycle)
+    {
+      cout << "Cycle " << cycle << ':' << endl;
+
+      if (cycle == 0)
+       {
+         GridGenerator::hyper_cube (triangulation, -1, 1);
+         triangulation.refine_global (4);
+       }
+      else
+       {
+         refine_grid ();
+       };
+      
+
+      cout << "   Number of active cells:       "
+          << triangulation.n_active_cells()
+          << endl;
+
+      setup_system ();
+
+      cout << "   Number of degrees of freedom: "
+          << dof_handler.n_dofs()
+          << endl;
+      
+      assemble_system ();
+      solve ();
+      output_results (cycle);
+    };
+
+  DataOut<dim>::EpsFlags eps_flags;
+  eps_flags.z_scaling = 4;
+  
+  DataOut<dim> data_out;
+  data_out.set_flags (eps_flags);
+
+  data_out.attach_dof_handler (dof_handler);
+  data_out.add_data_vector (solution, "solution");
+  data_out.build_patches ();
+  
+  ofstream output ("final-solution.gmv");
+  data_out.write_gmv (output);
+};
+
+
+
+                                // The ``main'' function is exactly
+                                // like in previous examples, with
+                                // the only difference in the name of
+                                // the main class that actually does
+                                // the computation.
+int main () 
+{
+  try
+    {
+      deallog.depth_console (0);
+
+      AdvectionProblem<2> advection_problem_2d;
+      advection_problem_2d.run ();
+    }
+  catch (exception &exc)
+    {
+      cerr << endl << endl
+          << "----------------------------------------------------"
+          << endl;
+      cerr << "Exception on processing: " << endl
+          << exc.what() << endl
+          << "Aborting!" << endl
+          << "----------------------------------------------------"
+          << endl;
+      return 1;
+    }
+  catch (...) 
+    {
+      cerr << endl << endl
+          << "----------------------------------------------------"
+          << endl;
+      cerr << "Unknown exception!" << endl
+          << "Aborting!" << endl
+          << "----------------------------------------------------"
+          << 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.