]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Write all remaining documentation. Change variable order to degree where necessary.
authorhartmann <hartmann@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 13 Sep 2001 12:26:32 +0000 (12:26 +0000)
committerhartmann <hartmann@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 13 Sep 2001 12:26:32 +0000 (12:26 +0000)
git-svn-id: https://svn.dealii.org/trunk@4995 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-10/step-10.cc

index 3b067f4e4b1a041d0a231eade427838f509d9be0..598dbf593ada2f827cc38f1c5daec40cf514a974 100644 (file)
@@ -79,7 +79,11 @@ void gnuplot_output()
                                   // So first generate a coarse
                                   // triangulation of the circle and
                                   // associate a suitable boundary
-                                  // description to it:
+                                  // description to it. Note that the
+                                  // default values of the
+                                  // HyperBallBoundary constructor
+                                  // are a center at the origin and a
+                                  // radius equals one.
   Triangulation<dim> triangulation;
   GridGenerator::hyper_ball (triangulation);
   static const HyperBallBoundary<dim> boundary;
@@ -128,9 +132,9 @@ void gnuplot_output()
 
                                       // Then output the present grid
                                       // for Q1, Q2, and Q3 mappings:
-      for (unsigned int order=1; order<4; ++order)
+      for (unsigned int degree=1; degree<4; ++degree)
        {
-         std::cout << "Order = " << order << std::endl;
+         std::cout << "Degree = " << degree << std::endl;
 
                                           // For this, first set up
                                           // an object describing the
@@ -139,9 +143,9 @@ void gnuplot_output()
                                           // class, which takes as
                                           // argument to the
                                           // constructor the
-                                          // polynomial order which
+                                          // polynomial degree which
                                           // it shall use.
-         const MappingQ<dim> mapping (order);
+         const MappingQ<dim> mapping (degree);
                                           // We note one interesting
                                           // fact: if you want a
                                           // piecewise linear
@@ -156,7 +160,7 @@ void gnuplot_output()
                                           // called ``MappingQ1''
                                           // which does exactly the
                                           // same is if you gave an
-                                          // order of ``1'' to the
+                                          // degree of ``1'' to the
                                           // ``MappingQ'' class, but
                                           // does so significantly
                                           // faster. ``MappingQ1'' is
@@ -169,7 +173,7 @@ void gnuplot_output()
                                           // explicitly.
 
 
-                                          // In order to actually
+                                          // In degree to actually
                                           // write out the present
                                           // grid with this mapping,
                                           // we set up an object
@@ -185,7 +189,7 @@ void gnuplot_output()
                                           // but since we want to
                                           // explicitely see the
                                           // effect of the mapping,
-                                          // we want to have teh
+                                          // we want to have the
                                           // faces in more
                                           // detail. This can be done
                                           // by passing the output
@@ -213,7 +217,7 @@ void gnuplot_output()
                                           // output using the same
                                           // evil hack as above:
          std::string filename = filename_base+"_mapping_q";
-         filename += ('0'+order);
+         filename += ('0'+degree);
          filename += ".dat";
          std::ofstream gnuplot_file (filename.c_str());
 
@@ -239,61 +243,206 @@ void gnuplot_output()
     }
 }
 
-
-
+                                // Now we proceed with the main part
+                                // of the code, the approximation of
+                                // pi. The area of a circle is given
+                                // by pi*radius^2, so having a circle
+                                // of radius 1, the area represents
+                                // just the number that is searched
+                                // for. The numerical computation of
+                                // the area is performed by
+                                // integrating the constant function
+                                // of value 1 over the whole
+                                // computational domain, i.e. by
+                                // computing the areas $\int_K 1
+                                // dx=\int_{\hat K} 1 J(\hat x) d\hat
+                                // x\approx\sum J(\hat x)w(\hat x)$
+                                // of all active cells of
+                                // triangulation and summing up these
+                                // contributions to gain the area of
+                                // the overall domain. The integrals
+                                // on each cell are approximated by
+                                // numerical quadrature, hence the
+                                // only additional ingredient we need
+                                // is to set up a FEValues object
+                                // that provides the corresponding
+                                // `JxW' values of each cell. We note
+                                // that here we won't use the
+                                // FEValues object in its original
+                                // purpose that is computing the
+                                // values of basis functions of a
+                                // specific finite element. But here
+                                // we use it only to gain the `JxW'
+                                // at the quadrature points,
+                                // irrespective of the (dummy) finite
+                                // element we will give to the
+                                // constructor of the FEValues
+                                // object.
 template <int dim>
 void compute_pi_by_area ()
 {
   std::cout << "Computation of Pi by the area:" << std::endl
            << "==============================" << std::endl;
-  
+
+                                  // For the numerical quadrature on
+                                  // all cells we employ a quadrature
+                                  // rule of sufficiently high
+                                  // degree. We choose QGauss4 that is
+                                  // of order 8, to be sure that the
+                                  // errors due to numerical
+                                  // quadrature are of higher order
+                                  // than the order (maximal 6) that
+                                  // will occur due to the order of
+                                  // the approximation of the
+                                  // boundary, i.e. the order of the
+                                  // mappings employed.
   const QGauss4<dim> quadrature;
-  for (unsigned int order=1; order<5; ++order)
+
+                                  // Now start by looping over
+                                  // degrees=1..4
+  for (unsigned int degree=1; degree<5; ++degree)
     {
-      std::cout << "Order = " << order << std::endl;
+      std::cout << "Degree = " << degree << std::endl;
+
+                                      // Then we generate the
+                                      // triangulation, the Boundary
+                                      // and the Mapping object as
+                                      // already seen.
       Triangulation<dim> triangulation;
       GridGenerator::hyper_ball (triangulation);
   
       static const HyperBallBoundary<dim> boundary;
       triangulation.set_boundary (0, boundary);
 
-      const MappingQ<dim> mapping (order);
-      const FE_Q<dim>     fe (1);
-
+      const MappingQ<dim> mapping (degree);
+
+                                      // We now create a dummy finite
+                                      // element. Here we could
+                                      // choose a finite element no
+                                      // matter which, as we are only
+                                      // interested in the `JxW'
+                                      // values provided by the
+                                      // FEValues object below.
+      const FE_Q<dim>     dummy_fe (1);
+
+                                      // Then we create a DofHandler
+                                      // object. This object will
+                                      // provide us with
+                                      // `active_cell_iterators' that
+                                      // are needed to reinit the
+                                      // FEValues object on each cell
+                                      // of the triangulation.
       DoFHandler<dim> dof_handler (triangulation);
-  
-      FEValues<dim> fe_values (mapping, fe, quadrature, update_JxW_values);
+
+                                      // Now we set up the FEValues
+                                      // object, giving the Mapping,
+                                      // the dummy finite element and
+                                      // the quadrature object to the
+                                      // constructor, together with
+                                      // the UpdateFlag asking for
+                                      // the `JxW' values at the
+                                      // quadrature points only.
+      FEValues<dim> fe_values (mapping, dummy_fe, quadrature, update_JxW_values);
+
+                                      // We employ an object of the
+                                      // ConvergenceTable class to
+                                      // store all important data
+                                      // like the approximative
+                                      // values for pi and the error
+                                      // wrt. the true value of
+                                      // pi. We will use functions
+                                      // provided by the
+                                      // ConvergenceTable class to
+                                      // compute convergence rates of
+                                      // the approximations to pi.
       ConvergenceTable table;
-      
+
+                                      // Now we loop over several
+                                      // refinement steps of the
+                                      // triangulation.
       for (unsigned int refinement=0; refinement<6;
           ++refinement, triangulation.refine_global (1))
        {
+                                          // In this loop we first
+                                          // add the number of active
+                                          // cells of the current
+                                          // triangulation to the
+                                          // table. This function
+                                          // automatically creates a
+                                          // table column with
+                                          // superscription `cells',
+                                          // for the case this column
+                                          // was not created before.
          table.add_value("cells", triangulation.n_active_cells());
-           
-         dof_handler.distribute_dofs (fe);  
-         
+
+                                          // Then we distribute the
+                                          // degrees of freedoms for
+                                          // the dummy finite
+                                          // element. Strictly
+                                          // speaking we do not need
+                                          // this function call in
+                                          // our special case but we
+                                          // call it to make the
+                                          // DoFHandler happy --
+                                          // otherwise it would throw
+                                          // an assertion in the
+                                          // FEValues::reinit
+                                          // function below.
+         dof_handler.distribute_dofs (dummy_fe);
+
+                                          // We define the variable
+                                          // area as `long double'
+                                          // like we did for the pi
+                                          // variable before.
+         long double area = 0;
+
+                                          // Now we loop over all
+                                          // cells, reinit the
+                                          // FEValues object for each
+                                          // cell, add all `JxW'
+                                          // values to `area'
          typename DoFHandler<dim>::active_cell_iterator
            cell = dof_handler.begin_active(),
            endc = dof_handler.end();
-         long double area = 0;
          for (; cell!=endc; ++cell)
            {
              fe_values.reinit (cell);
              for (unsigned int i=0; i<fe_values.n_quadrature_points; ++i)
                area += fe_values.JxW (i);
            };
+
+                                          // and store the resulting
+                                          // area values and the
+                                          // errors in the table. We
+                                          // need a static cast to
+                                          // double as there is no
+                                          // add_value(string, long
+                                          // double) function
+                                          // implemented.
          table.add_value("eval.pi", static_cast<double> (area));
          table.add_value("error", fabs(area-pi));
        };
 
+                                      // We want to compute
+                                      // the convergence rates of the
+                                      // `error' column. Therefore we
+                                      // need to omit the other
+                                      // columns from the convergence
+                                      // rate evaluation before
+                                      // calling
+                                      // `evaluate_all_convergence_rates'
       table.omit_column_from_convergence_rate_evaluation("cells");
       table.omit_column_from_convergence_rate_evaluation("eval.pi");
       table.evaluate_all_convergence_rates(
        ConvergenceTable::reduction_rate_log2);
 
+                                      // Finally we set the precision
+                                      // and the scientific mode
       table.set_precision("eval.pi", 16);
       table.set_scientific("error", true);
 
+                                      // and write the whole table to
+                                      // cout.
       table.write_text(std::cout);
 
       std::cout << std::endl;
@@ -301,28 +450,51 @@ void compute_pi_by_area ()
 };
 
 
-
+                                // The following function also
+                                // computes an approximation of pi
+                                // but this time via the diameter
+                                // 2*pi*radius of the domain instead
+                                // of the area. This function is only
+                                // a variation of the previous
+                                // function. So we will mainly give
+                                // documentation for the differences.
 template <int dim>
 void compute_pi_by_perimeter ()
 {
   std::cout << "Computation of Pi by the perimeter:" << std::endl
            << "===================================" << std::endl;
 
+                                  // We take the same order of
+                                  // quadrature but this time a
+                                  // `dim-1' dimensional quadrature
+                                  // as we will integrate over
+                                  // (boundary) lines rather than
+                                  // over cells.
   const QGauss4<dim-1> quadrature;
-  for (unsigned int order=1; order<5; ++order)
+
+                                  // We loop over all degrees, create
+                                  // the Triangulation, the Boundary,
+                                  // the Mapping, the dummy
+                                  // FiniteElement and the DoFHandler
+                                  // object as seen before.
+  for (unsigned int degree=1; degree<5; ++degree)
     {
-      std::cout << "Order = " << order << std::endl;
+      std::cout << "Degree = " << degree << std::endl;
       Triangulation<dim> triangulation;
       GridGenerator::hyper_ball (triangulation);
   
       static const HyperBallBoundary<dim> boundary;
       triangulation.set_boundary (0, boundary);
 
-      const MappingQ<dim> mapping (order);
+      const MappingQ<dim> mapping (degree);
       const FE_Q<dim>     fe (1);
 
       DoFHandler<dim> dof_handler (triangulation);
-      
+
+                                      // Then we create a FEFaceValues
+                                      // object instead of a FEValues
+                                      // object as in the previous
+                                      // function.
       FEFaceValues<dim> fe_face_values (mapping, fe, quadrature, update_JxW_values);
       ConvergenceTable table;
 
@@ -332,7 +504,15 @@ void compute_pi_by_perimeter ()
          table.add_value("cells", triangulation.n_active_cells());
 
          dof_handler.distribute_dofs (fe);
-         
+
+                                          // Now we run over all
+                                          // cells and over all faces
+                                          // of each cell. Only the
+                                          // contributions of the
+                                          // `JxW' values on boundary
+                                          // faces are added to the
+                                          // long double variable
+                                          // `perimeter'.
          typename DoFHandler<dim>::active_cell_iterator
            cell = dof_handler.begin_active(),
            endc = dof_handler.end();
@@ -341,14 +521,25 @@ void compute_pi_by_perimeter ()
            for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
              if (cell->face(face_no)->at_boundary())
                {
+                                                  // We reinit the
+                                                  // FEFaceValues
+                                                  // object with the
+                                                  // cell iterator
+                                                  // and the number
+                                                  // of the face.
                  fe_face_values.reinit (cell, face_no);
                  for (unsigned int i=0; i<fe_face_values.n_quadrature_points; ++i)
                    perimeter += fe_face_values.JxW (i);
                };
+                                          // We store the evaluated
+                                          // values in the table
          table.add_value("eval.pi", static_cast<double> (perimeter/2.));
          table.add_value("error", fabs(perimeter/2.-pi));
        };
 
+                                      // and we end this function as
+                                      // we did in the previous
+                                      // function.
       table.omit_column_from_convergence_rate_evaluation("cells");
       table.omit_column_from_convergence_rate_evaluation("eval.pi");
       table.evaluate_all_convergence_rates(
@@ -364,6 +555,9 @@ void compute_pi_by_perimeter ()
 };
 
 
+                                // The following main function just
+                                // calles the above functions in the
+                                // order of appearance.
 int main () 
 {
   std::cout.precision (16);

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.