From 2100652f6f923d757b51abd064026c705c441631 Mon Sep 17 00:00:00 2001 From: hartmann 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 triangulation; GridGenerator::hyper_ball (triangulation); static const HyperBallBoundary 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 mapping (order); + const MappingQ 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 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 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 triangulation; GridGenerator::hyper_ball (triangulation); static const HyperBallBoundary boundary; triangulation.set_boundary (0, boundary); - const MappingQ mapping (order); - const FE_Q fe (1); - + const MappingQ 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 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 dof_handler (triangulation); - - FEValues 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 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::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 (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 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 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 triangulation; GridGenerator::hyper_ball (triangulation); static const HyperBallBoundary boundary; triangulation.set_boundary (0, boundary); - const MappingQ mapping (order); + const MappingQ mapping (degree); const FE_Q fe (1); DoFHandler dof_handler (triangulation); - + + // Then we create a FEFaceValues + // object instead of a FEValues + // object as in the previous + // function. FEFaceValues 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::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::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 (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