From 24945ea228caace81e31f32625e3f5123ef32ae4 Mon Sep 17 00:00:00 2001 From: wolf 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 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 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 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 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::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 (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 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 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 (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