From 24945ea228caace81e31f32625e3f5123ef32ae4 Mon Sep 17 00:00:00 2001
From: wolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Fri, 14 Sep 2001 08:27:02 +0000
Subject: [PATCH] Minor updates.

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

diff --git a/deal.II/examples/step-10/step-10.cc b/deal.II/examples/step-10/step-10.cc
index 5506b4cd3a..d597be3be9 100644
--- a/deal.II/examples/step-10/step-10.cc
+++ b/deal.II/examples/step-10/step-10.cc
@@ -255,29 +255,41 @@ void gnuplot_output()
 				 // 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 the
-				 // 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
+				 // dx=\int_{\hat K} 1 det J(\hat x)
+				 // d\hat x \approx \sum det J(\hat
+				 // xi)w(\hat xi)$, where the sum
+				 // extends over all quadrature points
+				 // on all active cells in the
+				 // triangulation, with $w(xi)$ being
+				 // the weight of quadrature point
+				 // $xi$. 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. (Note that
+				 // `JxW' is meant to abbreviate
+				 // ``Jacobian determinant times
+				 // weight''; since in numerical
+				 // quadrature the two factors always
+				 // occur at the same places, we only
+				 // offer the combined quantity,
+				 // rather than two separate ones.) 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.
+				 // purpose, i.e. for the computation
+				 // of values of basis functions of a
+				 // specific finite element at certain
+				 // quadrature points. Rather, 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. The actual finite
+				 // element given to the FEValues
+				 // object is not used at all, so we
+				 // could give any.
 template <int dim>
 void compute_pi_by_area ()
 {
@@ -287,24 +299,34 @@ void compute_pi_by_area ()
 				   // 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
+				   // 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.
+				   // mappings employed. Note that the
+				   // integrand, the Jacobian
+				   // determinant, is not a polynomial
+				   // function (rather, it is a
+				   // rational one), so we do not use
+				   // Gauss quadrature in order to get
+				   // the exact value of the integral
+				   // as done often in finite element
+				   // computations, but could as well
+				   // have used any quadrature formula
+				   // of like order instead.
   const QGauss4<dim> quadrature;
 
 				   // Now start by looping over
-				   // degrees=1..4
+				   // polynomial mapping degrees=1..4:
   for (unsigned int degree=1; degree<5; ++degree)
     {
       std::cout << "Degree = " << degree << std::endl;
 
-				       // Then we generate the
+				       // First generate the
 				       // triangulation, the Boundary
 				       // and the Mapping object as
 				       // already seen.
@@ -318,16 +340,29 @@ void compute_pi_by_area ()
 
 				       // 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.
+				       // choose any finite element,
+				       // as we are only interested in
+				       // the `JxW' values provided by
+				       // the FEValues object
+				       // below. Nevertheless, we have
+				       // to provide a finite element
+				       // since in this example we
+				       // abuse the FEValues class a
+				       // little in that we only ask
+				       // it to provide us with the
+				       // weights of certain
+				       // quadrature points, in
+				       // contrast to the usual
+				       // purpose (and name) of the
+				       // FEValues class which is to
+				       // provide the values of finite
+				       // elements at these points.
       const FE_Q<dim>     dummy_fe (1);
 
-				       // Then we create a DofHandler
-				       // object. This object will
-				       // provide us with
+				       // Likewise, we need to create
+				       // a DofHandler object. We do
+				       // not actually use it, but it
+				       // will provide us with
 				       // `active_cell_iterators' that
 				       // are needed to reinit the
 				       // FEValues object on each cell
@@ -339,19 +374,40 @@ void compute_pi_by_area ()
 				       // the dummy finite element and
 				       // the quadrature object to the
 				       // constructor, together with
-				       // the UpdateFlag asking for
+				       // the UpdateFlags asking for
 				       // the `JxW' values at the
-				       // quadrature points only.
+				       // quadrature points only. This
+				       // tells the FEValues object
+				       // that it needs not compute
+				       // other quantities upon
+				       // calling the ``reinit''
+				       // function, thus saving
+				       // computation time.
+				       //
+				       // The most important
+				       // difference in the
+				       // construction of the FEValues
+				       // object compared to previous
+				       // example programs is that we
+				       // pass a mapping object as
+				       // first argument, which is to
+				       // be used in the computation
+				       // of the mapping from unit to
+				       // real cell. In previous
+				       // examples, this argument was
+				       // omitted, resulting in the
+				       // implicit use of an object of
+				       // type MappingQ1.
       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
+				       // like the approximated values
+				       // for pi and the error with
+				       // respect to the true value of
+				       // pi. We will also use
+				       // functions provided by the
 				       // ConvergenceTable class to
 				       // compute convergence rates of
 				       // the approximations to pi.
@@ -371,12 +427,12 @@ void compute_pi_by_area ()
 					   // automatically creates a
 					   // table column with
 					   // superscription `cells',
-					   // for the case this column
-					   // was not created before.
+					   // in case this column was
+					   // not created before.
 	  table.add_value("cells", triangulation.n_active_cells());
 
 					   // Then we distribute the
-					   // degrees of freedoms for
+					   // degrees of freedom for
 					   // the dummy finite
 					   // element. Strictly
 					   // speaking we do not need
@@ -397,10 +453,11 @@ void compute_pi_by_area ()
 	  long double area = 0;
 
 					   // Now we loop over all
-					   // cells, reinit the
+					   // cells, reinitialize the
 					   // FEValues object for each
-					   // cell, add all `JxW'
-					   // values to `area'
+					   // cell, and add up all the
+					   // `JxW' values for this
+					   // cell to `area'...
 	  typename DoFHandler<dim>::active_cell_iterator
 	    cell = dof_handler.begin_active(),
 	    endc = dof_handler.end();
@@ -411,13 +468,13 @@ void compute_pi_by_area ()
 		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
+					   // ...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));
@@ -433,16 +490,17 @@ void compute_pi_by_area ()
 				       // `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);
+      table.evaluate_all_convergence_rates(ConvergenceTable::reduction_rate_log2);
 
 				       // Finally we set the precision
-				       // and the scientific mode
+				       // and scientific mode for
+				       // output of some of the
+				       // quantities...
       table.set_precision("eval.pi", 16);
       table.set_scientific("error", true);
 
-				       // and write the whole table to
-				       // cout.
+				       // ...and write the whole table
+				       // to std::cout.
       table.write_text(std::cout);
 
       std::cout << std::endl;
@@ -450,9 +508,9 @@ void compute_pi_by_area ()
 };
 
 
-				 // The following function also
+				 // The following, second function also
 				 // computes an approximation of pi
-				 // but this time via the diameter
+				 // but this time via the perimeter
 				 // 2*pi*radius of the domain instead
 				 // of the area. This function is only
 				 // a variation of the previous
@@ -491,10 +549,12 @@ void compute_pi_by_perimeter ()
 
       DoFHandler<dim> dof_handler (triangulation);
 
-				       // Then we create a FEFaceValues
-				       // object instead of a FEValues
-				       // object as in the previous
-				       // function.
+				       // Then we create a
+				       // FEFaceValues object instead
+				       // of a FEValues object as in
+				       // the previous
+				       // function. Again, we pass a
+				       // mapping as first argument.
       FEFaceValues<dim> fe_face_values (mapping, fe, quadrature, update_JxW_values);
       ConvergenceTable table;
 
@@ -531,19 +591,17 @@ void compute_pi_by_perimeter ()
 		  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
+					   // Then 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.
+				       // ...and end this function as
+				       // we did in the previous one:
       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);
+      table.evaluate_all_convergence_rates(ConvergenceTable::reduction_rate_log2);
 
       table.set_precision("eval.pi", 16);
       table.set_scientific("error", true);
@@ -556,8 +614,8 @@ void compute_pi_by_perimeter ()
 
 
 				 // The following main function just
-				 // calles the above functions in the
-				 // order of appearance.
+				 // calls the above functions in the
+				 // order of their appearance.
 int main () 
 {
   std::cout.precision (16);
-- 
2.39.5