]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add somethings
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 7 Sep 2006 23:02:51 +0000 (23:02 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 7 Sep 2006 23:02:51 +0000 (23:02 +0000)
git-svn-id: https://svn.dealii.org/trunk@13854 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-23/doc/intro.dox
deal.II/examples/step-23/step-23.cc

index bfb31d8df4090000c19b80b4cbb8fe93d5d03974..78172918f5c50692e9c9aad565c3818b928a8c4e 100644 (file)
@@ -193,6 +193,10 @@ time stepping schemes do not conserve it after time
 discretization). Since $v^n$ also appears in the equation for $u^n$,
 the Crank-Nicolson scheme is also implicit.
 
+In the program, we will leave $\theta$ as a parameter, so that it will
+be easy to play with it. The results section will show some numerical
+evidence comparing the different schemes.
+
 The equations above (called the <i>semidiscretized</i> equations
 because we have only discretized the time, but not space), can be
 simplified a bit by eliminating $v^n$ from the first equation and
@@ -306,6 +310,10 @@ boundary values for the first equation above (i.e. the one for $U^n$),
 whereas we do not have to do that for the second one.
 
 
+<h3>Who are Courant, Friedrichs, and Levy?<h3>
+
+CFL condition
+
 <h3>How the program works</h3>
 
 Given the above formulation, ....
index 4afa1129b86444086f0dbf8b6fc7ccb4fb024bcb..c7a9739b69957e107ae19e1be38ca5c602e8b27b 100644 (file)
@@ -1,10 +1,10 @@
 /* $Id: step-4.cc,v 1.34 2006/02/06 21:33:10 wolf Exp $ */
-/* Author: Wolfgang Bangerth, University of Heidelberg, 1999 */
+/* Author: Wolfgang Bangerth, Texas A&M University, 2006 */
 
 /*    $Id: step-4.cc,v 1.34 2006/02/06 21:33:10 wolf Exp $       */
 /*    Version: $Name:  $                                          */
 /*                                                                */
-/*    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors */
+/*    Copyright (C) 2006 by the deal.II authors */
 /*                                                                */
 /*    This file is subject to QPL and may not be  distributed     */
 /*    without copyright and license information. Please refer     */
 /*    further information on this license.                        */
 
 
-#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_q.h>
-#include <dofs/dof_tools.h>
-#include <fe/fe_values.h>
+                                // @sect3{Include files}
+
+                                // We start with the usual assortment
+                                // of include files that we've seen
+                                // in so many of the previous tests:
 #include <base/quadrature_lib.h>
 #include <base/function.h>
-#include <numerics/vectors.h>
-#include <numerics/matrices.h>
+#include <base/logstream.h>
+
 #include <lac/vector.h>
 #include <lac/full_matrix.h>
 #include <lac/sparse_matrix.h>
 #include <lac/solver_cg.h>
 #include <lac/precondition.h>
+
+#include <grid/tria.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+
+#include <dofs/dof_handler.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
 #include <dofs/dof_constraints.h>
 
-#include <numerics/matrices.h>
-#include <numerics/vectors.h>
+#include <fe/fe_q.h>
+#include <fe/fe_values.h>
 
 #include <numerics/data_out.h>
+
 #include <fstream>
 #include <iostream>
 #include <sstream>
 
-#include <base/logstream.h>
-
+                                // Here are the only three include
+                                // files of some new interest: The
+                                // first one is already used, for
+                                // example, for the
+                                // VectorTools::interpolate_boundary_values
+                                // and
+                                // VectorTools::apply_boundary_values
+                                // functions. However, we here use
+                                // another function in that class,
+                                // VectorTools::project to compute
+                                // our initial values as the $L^2$
+                                // projection of the continuous
+                                // initial values. Furthermore, we
+                                // use
+                                // VectorTools::create_right_hand_side
+                                // to generate the integrals
+                                // $(f^n,\phi^n_i)$. These were
+                                // previously always generated by
+                                // hand in
+                                // <code>assemble_system</code> or
+                                // similar functions in application
+                                // code. However, we're too lazy to
+                                // do that here, so simply use a
+                                // library function:
+#include <numerics/vectors.h>
 
+                                // In a very similar vein, we are
+                                // also too lazy to write the code to
+                                // assemble mass and Laplace
+                                // matrices, although it would have
+                                // only taken copying the relevant
+                                // code from any number of previous
+                                // tutorial programs. Rather, we want
+                                // to focus on the things that are
+                                // truly new to this program and
+                                // therefore use the
+                                // MatrixTools::create_mass_matrix
+                                // and
+                                // MatrixTools::create_laplace_matrix
+                                // functions. They are declared here:
+#include <numerics/matrices.h>
 
+                                // Finally, here is an include file
+                                // that contains all sorts of tool
+                                // functions that one sometimes
+                                // needs. In particular, we need the
+                                // Utilities::int_to_string class
+                                // that, given an integer argument,
+                                // returns a string representation of
+                                // it. It is particularly useful
+                                // since it allows for a second
+                                // parameter indicating the number of
+                                // digits to which we want the result
+                                // padded with leading zeros. We will
+                                // use this to write output files
+                                // that have the form
+                                // <code>solution-XXX.gnuplot</code>
+                                // where <code>XXX</code> denotes the
+                                // number of the time step and always
+                                // consists of three digits even if
+                                // we are still in the single or
+                                // double digit time steps.
+#include <base/utilities.h>
+
+
+                                // @sect3{The <code>WaveEquation</code> class}
+
+                                // Next comes the declaration of the
+                                // main class. It's public interface
+                                // of functions is like in most of
+                                // the other tutorial programs. Worth
+                                // mentioning is that we now have to
+                                // store three matrices instead of
+                                // one: the mass matrix $M$, the
+                                // Laplace matrix $A$, and the system
+                                // matrix $M+k^2\theta^2A$ used when
+                                // solving for $U^n$. Likewise, we
+                                // need solution vectors for
+                                // $U^n,V^n$ as well as for the
+                                // corresponding vectors at the
+                                // previous time step,
+                                // $U^{n-1},V^{n-1}$. The
+                                // <code>system_rhs</code> will be
+                                // used for whatever right hand side
+                                // vector we have when solving one of
+                                // the two linear systems we have to
+                                // solve in each time step. These
+                                // will be solved in the two
+                                // functions <code>solve_u</code> and
+                                // <code>solve_v</code>.
+                                //
+                                // Finally, the variable
+                                // <code>theta</code> is used to
+                                // indicate the parameter $\theta$
+                                // that is used to define which time
+                                // stepping scheme to use. The rest
+                                // is self-explanatory.
 template <int dim>
-class WaveEquationProblem 
+class WaveEquation 
 {
   public:
-    WaveEquationProblem ();
+    WaveEquation ();
     void run ();
     
   private:
-    void make_grid_and_dofs ();
-    void assemble_system ();
+    void setup_system ();
     void solve_u ();
     void solve_v ();
-    void output_results (const unsigned int timestep_number) const;
+    void output_results () const;
 
     Triangulation<dim>   triangulation;
     FE_Q<dim>            fe;
     DoFHandler<dim>      dof_handler;
 
+    ConstraintMatrix constraints;
+    
     SparsityPattern      sparsity_pattern;
     SparseMatrix<double> system_matrix;
     SparseMatrix<double> mass_matrix;
     SparseMatrix<double> laplace_matrix;
 
-    double time_step;
-    double theta;
-
     Vector<double>       solution_u, solution_v;
     Vector<double>       old_solution_u, old_solution_v;
     Vector<double>       system_rhs;
+
+    double time, time_step;
+    unsigned int timestep_number;
+    const double theta;
 };
 
 
 
+                                // @sect3{Equation data}
+
+                                // Before we go on filling in the
+                                // details of the main class, let us
+                                // define the equation data
+                                // corresponding to the problem,
+                                // i.e. initial and boundary values
+                                // as well as a right hand side
+                                // class. We do so using classes
+                                // derived from the Function class
+                                // template that has been used many
+                                // times before, so the following
+                                // should not be a surprise.
+                                //
+                                // Let's start with initial values
+                                // and choose zero for both the value
+                                // $u$ as well as its time
+                                // derivative, the velocity $v$:
 template <int dim>
-class RightHandSide : public Function<dim> 
+class InitialValuesU : public Function<dim> 
 {
   public:
-    RightHandSide () : Function<dim>() {};
+    InitialValuesU () : Function<dim>() {};
     
     virtual double value (const Point<dim>   &p,
                          const unsigned int  component = 0) const;
@@ -89,10 +208,10 @@ class RightHandSide : public Function<dim>
 
 
 template <int dim>
-class InitialValuesU : public Function<dim> 
+class InitialValuesV : public Function<dim> 
 {
   public:
-    InitialValuesU () : Function<dim>() {};
+    InitialValuesV () : Function<dim>() {};
     
     virtual double value (const Point<dim>   &p,
                          const unsigned int  component = 0) const;
@@ -101,10 +220,33 @@ class InitialValuesU : public Function<dim>
 
 
 template <int dim>
-class BoundaryValues : public Function<dim> 
+double InitialValuesU<dim>::value (const Point<dim>  &/*p*/,
+                                  const unsigned int component) const 
+{
+  Assert (component == 0, ExcInternalError());
+  return 0;
+}
+
+
+
+template <int dim>
+double InitialValuesV<dim>::value (const Point<dim>  &/*p*/,
+                                  const unsigned int component) const 
+{
+  Assert (component == 0, ExcInternalError());
+  return 0;
+}
+
+
+
+                                // Secondly, we have the right hand
+                                // side forcing term. Boring as we
+                                // are, we choose zero here as well:
+template <int dim>
+class RightHandSide : public Function<dim> 
 {
   public:
-    BoundaryValues () : Function<dim>() {};
+    RightHandSide () : Function<dim>() {};
     
     virtual double value (const Point<dim>   &p,
                          const unsigned int  component = 0) const;
@@ -112,45 +254,55 @@ class BoundaryValues : public Function<dim>
 
 
 
-
 template <int dim>
-double RightHandSide<dim>::value (const Point<dim> &/*p*/,
-                                 const unsigned int /*component*/) const 
+double RightHandSide<dim>::value (const Point<dim>  &/*p*/,
+                                 const unsigned int component) const 
 {
-//   if (get_time() <= 0.25)
-//     if ((p[0] <=0) && (p[1] <= 0))
-//       return 1;
-
+  Assert (component == 0, ExcInternalError());
   return 0;
 }
 
 
+
+                                // Finally, we have boundary
+                                // values. They are as described in
+                                // the introduction:
 template <int dim>
-double InitialValuesU<dim>::value (const Point<dim> &p,
-                                 const unsigned int /*component*/) const 
+class BoundaryValues : public Function<dim> 
 {
-  //  return std::sqrt(p.square()) * std::exp (-p.square()) / 3;
-  if ((p[0] <=0) && (p[1] <= 0))
-    return 1;
+  public:
+    BoundaryValues () : Function<dim>() {};
+    
+    virtual double value (const Point<dim>   &p,
+                         const unsigned int  component = 0) const;
+};
+
 
-  return 0;}
 
 
 template <int dim>
-double BoundaryValues<dim>::value (const Point<dim> &/*p*/,
-                                  const unsigned int /*component*/) const 
+double BoundaryValues<dim>::value (const Point<dim> &p,
+                                  const unsigned int component) const 
 {
-  return 0;
+  Assert (component == 0, ExcInternalError());
+
+  if ((this->get_time() <= 1) &&
+      (p[0] < 1) &&
+      (p[1] < 1./3) &&
+      (p[1] > -1./3))
+    return std::sin (this->get_time() * 2 * deal_II_numbers::PI);
+  else
+    return 0;
 }
 
 
 
 
-
+                                // @sect3{Implementation of the <code>WaveEquation</code> class}
 
 
 template <int dim>
-WaveEquationProblem<dim>::WaveEquationProblem () :
+WaveEquation<dim>::WaveEquation () :
                 fe (1),
                dof_handler (triangulation),
                time_step (1./64),
@@ -160,17 +312,18 @@ WaveEquationProblem<dim>::WaveEquationProblem () :
 
 
 template <int dim>
-void WaveEquationProblem<dim>::make_grid_and_dofs ()
+void WaveEquation<dim>::setup_system ()
 {
   GridGenerator::hyper_cube (triangulation, -1, 1);
   triangulation.refine_global (7);
   
-  std::cout << "   Number of active cells: "
+  std::cout << "Number of active cells: "
            << triangulation.n_active_cells()
            << std::endl
-           << "   Total number of cells: "
+           << "Total number of cells: "
            << triangulation.n_cells()
-           << std::endl;
+           << std::endl
+           << std::endl;
 
   dof_handler.distribute_dofs (fe);
 
@@ -188,18 +341,6 @@ void WaveEquationProblem<dim>::make_grid_and_dofs ()
   mass_matrix.reinit (sparsity_pattern);
   laplace_matrix.reinit (sparsity_pattern);
 
-  solution_u.reinit (dof_handler.n_dofs());
-  solution_v.reinit (dof_handler.n_dofs());
-  old_solution_u.reinit (dof_handler.n_dofs());
-  old_solution_v.reinit (dof_handler.n_dofs());
-  system_rhs.reinit (dof_handler.n_dofs());
-}
-
-
-
-template <int dim>
-void WaveEquationProblem<dim>::assemble_system () 
-{  
   MatrixCreator::create_mass_matrix (dof_handler, QGauss<dim>(3),
                                     mass_matrix);
   MatrixCreator::create_laplace_matrix (dof_handler, QGauss<dim>(3),
@@ -207,12 +348,20 @@ void WaveEquationProblem<dim>::assemble_system ()
   
   system_matrix.copy_from (mass_matrix);
   system_matrix.add (theta * theta * time_step * time_step, laplace_matrix);
+  
+  solution_u.reinit (dof_handler.n_dofs());
+  solution_v.reinit (dof_handler.n_dofs());
+  old_solution_u.reinit (dof_handler.n_dofs());
+  old_solution_v.reinit (dof_handler.n_dofs());
+  system_rhs.reinit (dof_handler.n_dofs());
+
+  constraints.close ();
 }
 
 
 
 template <int dim>
-void WaveEquationProblem<dim>::solve_u () 
+void WaveEquation<dim>::solve_u () 
 {
   SolverControl           solver_control (1000, 1e-8*system_rhs.l2_norm());
   SolverCG<>              cg (solver_control);
@@ -226,7 +375,7 @@ void WaveEquationProblem<dim>::solve_u ()
 
 
 template <int dim>
-void WaveEquationProblem<dim>::solve_v () 
+void WaveEquation<dim>::solve_v () 
 {
   SolverControl           solver_control (1000, 1e-8*system_rhs.l2_norm());
   SolverCG<>              cg (solver_control);
@@ -241,7 +390,7 @@ void WaveEquationProblem<dim>::solve_v ()
 
 
 template <int dim>
-void WaveEquationProblem<dim>::output_results (const unsigned int timestep_number) const
+void WaveEquation<dim>::output_results () const
 {
   DataOut<dim> data_out;
 
@@ -253,7 +402,7 @@ void WaveEquationProblem<dim>::output_results (const unsigned int timestep_numbe
 
   std::ostringstream filename;
   filename << "solution-"
-          << timestep_number
+          << Utilities::int_to_string (timestep_number, 3)
           << ".gnuplot";
   std::ofstream output (filename.str().c_str());
   data_out.write_gnuplot (output);
@@ -263,24 +412,18 @@ void WaveEquationProblem<dim>::output_results (const unsigned int timestep_numbe
 
 
 template <int dim>
-void WaveEquationProblem<dim>::run () 
+void WaveEquation<dim>::run () 
 {
-  std::cout << "Solving problem in " << dim << " space dimensions." << std::endl;
-  
-  make_grid_and_dofs();
-  assemble_system ();
+  setup_system();
 
-  ConstraintMatrix constraints;
-  constraints.close();
   VectorTools::project (dof_handler, constraints, QGauss<dim>(3),
                        InitialValuesU<dim>(),
                        old_solution_u);
   VectorTools::project (dof_handler, constraints, QGauss<dim>(3),
-                       ZeroFunction<dim>(),
+                       InitialValuesV<dim>(),
                        old_solution_v);
-
-  unsigned int timestep_number = 1;
-  for (double time = time_step; time<=5; time+=time_step, ++timestep_number)
+  
+  for (timestep_number=1, time=time_step; time<=5; time+=time_step, ++timestep_number)
     {
       std::cout << "Time step " << timestep_number
                << " at t=" << time
@@ -308,10 +451,13 @@ void WaveEquationProblem<dim>::run ()
       system_rhs.add (theta * (1-theta) * time_step * time_step, tmp);
 
 
+      BoundaryValues<dim> boundary_values_function;
+      boundary_values_function.set_time (time);
+      
       std::map<unsigned int,double> boundary_values;
       VectorTools::interpolate_boundary_values (dof_handler,
                                                0,
-                                               BoundaryValues<dim>(),
+                                               boundary_values_function,
                                                boundary_values);
       MatrixTools::apply_boundary_values (boundary_values,
                                          system_matrix,
@@ -341,7 +487,7 @@ void WaveEquationProblem<dim>::run ()
 
       solve_v ();
 
-      output_results (timestep_number);
+      output_results ();
 
       old_solution_u = solution_u;
       old_solution_v = solution_v;
@@ -354,8 +500,8 @@ int main ()
 {
   deallog.depth_console (0);
   {
-    WaveEquationProblem<2> wave_equation_problem_2d;
-    wave_equation_problem_2d.run ();
+    WaveEquation<2> wave_equation_solver;
+    wave_equation_solver.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.