From 2100652f6f923d757b51abd064026c705c441631 Mon Sep 17 00:00:00 2001
From: hartmann <hartmann@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Thu, 13 Sep 2001 12:26:32 +0000
Subject: [PATCH] Write all remaining documentation. Change variable order to
 degree where necessary.

git-svn-id: https://svn.dealii.org/trunk@4995 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/examples/step-10/step-10.cc | 254 ++++++++++++++++++++++++----
 1 file changed, 224 insertions(+), 30 deletions(-)

diff --git a/deal.II/examples/step-10/step-10.cc b/deal.II/examples/step-10/step-10.cc
index 3b067f4e4b..598dbf593a 100644
--- a/deal.II/examples/step-10/step-10.cc
+++ b/deal.II/examples/step-10/step-10.cc
@@ -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);
-- 
2.39.5