From: Wolfgang Bangerth Date: Thu, 9 Feb 2006 07:06:33 +0000 (+0000) Subject: Look through most of the program. X-Git-Tag: v8.0.0~12402 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=4bc416c69d7b38309dffb22f74ee1669c2ab2908;p=dealii.git Look through most of the program. git-svn-id: https://svn.dealii.org/trunk@12277 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-7/step-7.cc b/deal.II/examples/step-7/step-7.cc index 1ca3de6d01..fc9f0b79e5 100644 --- a/deal.II/examples/step-7/step-7.cc +++ b/deal.II/examples/step-7/step-7.cc @@ -40,46 +40,43 @@ #include #include - // In this example, we will not use - // the numeration scheme which is - // used per default by the - // ``DoFHandler'' class, but will - // renumber them using the - // Cuthill-McKee algorithm. The - // necessary functions are declared - // in the following file: + // In this example, we will not use the + // numeration scheme which is used per + // default by the ``DoFHandler'' class, but + // will renumber them using the Cuthill-McKee + // algorithm. As has already been explained + // in step-2, the necessary functions are + // declared in the following file: #include // Then we will show a little trick // how we can make sure that objects // are not deleted while they are // still in use. For this purpose, - // there is the ``SmartPointer'' + // deal.II has the ``SmartPointer'' // helper class, which is declared in // this file: #include - // Then we will want to use the + // Next, we will want to use the // ``integrate_difference'' function - // mentioned in the introduction. It - // comes from this file: + // mentioned in the introduction, and we are + // going to use a ``ConvergenceTable'' that + // collects all important data during a run + // and prints it at the end as a table. These + // comes from the following two files: #include - // We are going to use a - // ``ConvergenceTable'' that collects - // all important data during a run - // and prints it at the end as a - // table. #include // And finally, we need to use the // ``FEFaceValues'' class, which is - // declare in the same file as the + // declared in the same file as the // ``FEValues'' class: #include - // We need one more include from - // standard C++, which is necessary - // when we try to find out the actual - // type behind a pointer to a base - // class. We will explain this in - // slightly more detail below. + // We need one more include from standard + // C++, which is necessary when we try to + // find out the actual type behind a pointer + // to a base class. We will explain this in + // slightly more detail below. The other two + // include files are obvious then: #include #include #include @@ -87,42 +84,36 @@ // @sect3{Equation data} - // Before implementing the classes - // that actually solve something, we - // first declare and define some - // function classes that represent - // right hand side and solution - // classes. Since we want to compare - // the exactly known continuous - // solution to the computed one, we - // need a function object which - // represents the continuous - // solution. On the other hand, we - // need the right hand side function, - // and that one of course shares some - // characteristics with the - // solution. In order to reduce - // dependencies which arise if we - // have to change something in both - // classes at the same time, we - // exclude the common characteristics - // of both functions into a base - // class. + // Before implementing the classes that + // actually solve something, we first declare + // and define some function classes that + // represent right hand side and solution + // classes. Since we want to compare the + // numerically obtained solution to the exact + // continuous one, we need a function object + // that represents the continuous + // solution. On the other hand, we need the + // right hand side function, and that one of + // course shares some characteristics with + // the solution. In order to reduce + // dependencies which arise if we have to + // change something in both classes at the + // same time, we move the common + // characteristics of both functions into a + // base class. // - // The common characteristics for the - // given solution, which as explained - // in the introduction is a sum of - // three exponentials, are here: the - // number of exponentials, their - // centers, and their half width. We - // declare them in the following - // class. Since the number of - // exponentials is a constant scalar - // integral quantity, C++ allows its - // definition (i.e. assigning a - // value) right at the place of - // declaration (i.e. where we declare - // that such a variable exists). + // The common characteristics for solution + // (as explained in the introduction, we + // choose a sum of three exponentials) and + // right hand side, are these: the number of + // exponentials, their centers, and their + // half width. We declare them in the + // following class. Since the number of + // exponentials is a constant scalar integral + // quantity, C++ allows its definition + // (i.e. assigning a value) right at the + // place of declaration (i.e. where we + // declare that such a variable exists). template class SolutionBase { @@ -147,10 +138,25 @@ class SolutionBase // the 1d case for exposition of a // useful technique. // - // First we assign values to the - // centers for the 1d case, where we - // place the centers equidistantly at - // -1/3, 0, and 1/3: + // First we assign values to the centers for + // the 1d case, where we place the centers + // equidistantly at -1/3, 0, and 1/3. The + // ``template <>'' header for this definition + // indicates an explicit specialization. This + // means, that the variable belongs to a + // template, but that instead of providing + // the compiler with a template from which it + // can specialize a concrete variable by + // substituting ``dim'' with some concrete + // value, we provide a specialization + // ourselves, in this case for ``dim=1''. If + // the compiler then sees a reference to this + // variable in a place where the template + // argument equals one, it knows that it + // doesn't have to generate the variable from + // a template by substituting ``dim'', but + // can immediately use the following + // definition: template <> const Point<1> SolutionBase<1>::source_centers[SolutionBase<1>::n_source_centers] @@ -158,8 +164,9 @@ SolutionBase<1>::source_centers[SolutionBase<1>::n_source_centers] Point<1>(0.0), Point<1>(+1.0 / 3.0) }; - // Then we place the centers for the - // 2d case as follows: + // Likewise, we can provide an explicit + // specialization for ``dim=2''. We place the + // centers for the 2d case as follows: template <> const Point<2> SolutionBase<2>::source_centers[SolutionBase<2>::n_source_centers] @@ -167,11 +174,14 @@ SolutionBase<2>::source_centers[SolutionBase<2>::n_source_centers] Point<2>(-0.5, -0.5), Point<2>(+0.5, -0.5) }; - // There remains to assign a value to - // the half-width of the - // exponentials. We would like to use - // the same value for all dimensions, - // so here is how that works: + // There remains to assign a value to the + // half-width of the exponentials. We would + // like to use the same value for all + // dimensions. In this case, we simply + // provide the compiler with a template from + // which it can generate a concrete + // instantiation by substituting ``dim'' with + // a concrete value: template const double SolutionBase::width = 1./3.; @@ -188,17 +198,26 @@ const double SolutionBase::width = 1./3.; // the characteristics defined in the // ``SolutionBase'' class. // - // The actual classes are declared in - // the following. Note that in order - // to compute the error of the - // numerical solution against the - // continuous one in the L2 and H1 - // norms, we have to export value and - // gradient of the exact solution, - // which is done by overloading the - // respective virtual member - // functions in the ``Function'' base - // class. + // The actual classes are declared in the + // following. Note that in order to compute + // the error of the numerical solution + // against the continuous one in the L2 and + // H1 norms, we have to provide value and + // gradient of the exact solution. This is + // more than we have done in previous + // examples, where all we provided was the + // value at one or a list of + // points. Fortunately, the ``Function'' + // class also has virtual functions for the + // gradient, so we can simply overload the + // respective virtual member functions in the + // ``Function'' base class. Note that the + // gradient of a function in ``dim'' space + // dimensions is a vector of size ``dim'', + // i.e. a tensor of rank 1 and dimension + // ``dim''. As for so many other things, the + // library provides a suitable class for + // this. // // Just as in previous examples, we // are forced by the C++ language @@ -214,17 +233,16 @@ class Solution : public Function, virtual double value (const Point &p, const unsigned int component = 0) const; + virtual Tensor<1,dim> gradient (const Point &p, const unsigned int component = 0) const; }; - // The actual definition of the - // values and gradients of the exact - // solution class is according to - // their mathematical definition and - // probably needs not much - // explanation. + // The actual definition of the values and + // gradients of the exact solution class is + // according to their mathematical definition + // and does not need much explanation. // // The only thing that is worth // mentioning is that if we access @@ -238,7 +256,7 @@ class Solution : public Function, // qualification is not necessary if // the base class is not template // dependent, and also that the gcc - // compilers, among others, don't + // compilers prior to version 3.4 don't // enforce this requirement of the // C++ standard. The reason why this // is necessary is complicated; some @@ -253,70 +271,74 @@ double Solution::value (const Point &p, double return_value = 0; for (unsigned int i=0; in_source_centers; ++i) { - // One of the few things worth - // mentioning is the following - // variables, which represents - // the vector (x-x_i). It is - // computed in the way that one - // would intuitively expect: - const Point shifted_point = p-this->source_centers[i]; - - // The ``Point'' class - // offers a member function - // ``square'' that does what - // it's name suggests. - return_value += std::exp(-shifted_point.square() / + const Point x_minus_xi = p - this->source_centers[i]; + return_value += std::exp(-x_minus_xi.square() / (this->width * this->width)); - }; + } return return_value; } - + // Likewise, this is the computation of the + // gradient of the solution. In order to + // accumulate the gradient from the + // contributions of the exponentials, we + // allocate an object ``return_value'' that + // denotes the mathematical quantity of a + // tensor of rank ``1'' and dimension + // ``dim''. Its default constructor sets it + // to the vector containing only zeroes, so + // we need not explicitly care for its + // initialization. + // + // Note that we could as well have taken the + // type of the object to be ``Point'' + // instead of ``Tensor<1,dim>''. Tensors of + // rank 1 and points are almost exchangeable, + // and have only very slightly different + // mathematical meanings. In fact, the + // ``Point'' class is derived from the + // ``Tensor<1,dim>'' class, which makes up + // for their mutual exchange ability. Their + // main difference is in what they logically + // mean: points are points in space, such as + // the location at which we want to evaluate + // a function (see the type of the first + // argument of this function for example). On + // the other hand, tensors of rank 1 share + // the same transformation properties, for + // example that they need to be rotated in a + // certain way when we change the coordinate + // system; however, they do not share the + // same connotation that points have and are + // only objects in a more abstract space than + // the one spanned by the coordinate + // directions. (In fact, gradients live in + // `reciprocal' space, since the dimension of + // their components is not that of a length, + // but one over length). template Tensor<1,dim> Solution::gradient (const Point &p, const unsigned int) const { - // In order to accumulate the - // gradient from the contributions - // of the exponentials, we allocate - // an object which denotes the - // mathematical quantity of a - // tensor of rank ``1'' and - // dimension ``dim''. Its default - // constructor sets it to the - // vector containing only zeroes, - // so we need not explicitly care - // for its initialization. Tensor<1,dim> return_value; - // Note that we could as well have - // taken the type of the object to - // be ``Point''. Tensors of - // rank 1 and points are almost - // exchangeable, and have only very - // slightly different mathematical - // meanings. In fact, the - // ``Point'' class is derived - // from the ``Tensor<1,dim>'' - // class, which makes up for their - // mutual exchange ability. for (unsigned int i=0; in_source_centers; ++i) { - const Point shifted_point = p-this->source_centers[i]; + const Point x_minus_xi = p - this->source_centers[i]; // For the gradient, note that - // it's direction is along + // its direction is along // (x-x_i), so we add up // multiples of this distance // vector, where the factor is // given by the exponentials. return_value += (-2 / (this->width * this->width) * - std::exp(-shifted_point.square() / + std::exp(-x_minus_xi.square() / (this->width * this->width)) * - shifted_point); - }; + x_minus_xi); + } return return_value; } @@ -357,20 +379,20 @@ double RightHandSide::value (const Point &p, double return_value = 0; for (unsigned int i=0; in_source_centers; ++i) { - const Point shifted_point = p-this->source_centers[i]; + const Point x_minus_xi = p - this->source_centers[i]; // The first contribution is // the Laplacian: - return_value += ((2*dim - 4*shifted_point.square()/ + return_value += ((2*dim - 4*x_minus_xi.square()/ (this->width * this->width)) / (this->width * this->width) * - std::exp(-shifted_point.square() / + std::exp(-x_minus_xi.square() / (this->width * this->width))); // And the second is the // solution itself: - return_value += std::exp(-shifted_point.square() / + return_value += std::exp(-x_minus_xi.square() / (this->width * this->width)); - }; + } return return_value; } @@ -378,70 +400,58 @@ double RightHandSide::value (const Point &p, // @sect3{The Helmholtz solver class} - // Then we need the class that does - // all the work. It is mostly the - // same as in previous examples, and - // we will discuss the differences - // only when we declare the - // respective functions or variables - // below. + // Then we need the class that does all the + // work. Except for its name, its interface + // is mostly the same as in previous + // examples. + // + // One of the differences is that we will use + // this class in several modes: for different + // finite elements, as well as for adaptive + // and global refinement. The decision + // whether global or adaptive refinement + // shall be used is communicated to the + // constructor of this class through an + // enumeration type declared at the top of + // the class. The constructor then takes a + // finite element object and the refinement + // mode as arguments. + // + // The rest of the member functions are as + // before except for the ``process_solution'' + // function: After the solution has been + // computed, we perform some analysis on it, + // such as computing the error in various + // norms. To enable some output, it requires + // the number of the refinement cycle, and + // consequently gets it as an argument. template class HelmholtzProblem { public: - // We will use this class in - // several modes: for different - // finite elements, as well as - // for adaptive and global - // refinement. The decision - // whether global or adaptive - // refinement shall be used is - // communicated to the - // constructor of this class - // through an enumeration type, - // which we declare here: enum RefinementMode { global_refinement, adaptive_refinement }; - // This is the constructor of the - // class, it takes the finite - // element and the refinement - // mode as parameter and stores - // them in local variables. HelmholtzProblem (const FiniteElement &fe, const RefinementMode refinement_mode); - // The following two functions - // are the same as in previous - // examples. ~HelmholtzProblem (); void run (); private: - // As are these: void setup_system (); void assemble_system (); void solve (); void refine_grid (); - - // After the solution has been - // computed, we perform some - // analysis on it, such as - // computing the error in various - // norms. This is done in the - // following function. To enable - // some output, we pass it the - // number of the refinement - // cycle. void process_solution (const unsigned int cycle); // Now for the data elements of - // this class: - Triangulation triangulation; - DoFHandler dof_handler; - + // this class. Among the variables + // that we have already used in + // previous examples, only the + // finite element object differs: // The finite elements which the // objects of this class operate // on are passed to the @@ -478,7 +488,7 @@ class HelmholtzProblem // which the ``DoFHandler'' uses, // live at least as long as they // are in use? This means that - // the ``DoFHandler'' must have a + // the ``DoFHandler'' must have some // kind of lock on the // destruction of the other // objects, and it can only @@ -518,20 +528,22 @@ class HelmholtzProblem // subscribe to, thus the name of // the class. Whenever we // initialize a pointer to that - // object, we can increase it use + // object, we can increase its use // counter, and when we move away // our pointer or do not need it // any more, we decrease the // counter again. This way, we // can always check how many // objects still use that - // object. If an object of a + // object. + // + // On the other hand, if an object of a // class that is derived from the - // ``Subscriptor'' class is - // destroyed, it also has to call - // the destructor of the - // ``Subscriptor'' class; this - // will then check whether the + // ``Subscriptor'' class is destroyed, it + // also has to call the destructor of the + // ``Subscriptor'' class. In this + // destructor, there + // will then be a check whether the // counter is really zero. If // yes, then there are no active // references to this object any @@ -542,9 +554,9 @@ class HelmholtzProblem // stale and thus potentially // dangerous pointers, and we // rather throw an exception to - // alert the programmer that she + // alert the programmer that this // is doing something dangerous - // and better had her program + // and the program better be // fixed. // // While this certainly all @@ -560,12 +572,12 @@ class HelmholtzProblem // difficult to find bugs, since // the place where we have // forgotten something may be - // very far away from the place + // far away from the place // where the check for zeroness // of the counter upon // destruction actually // fails. This kind of bug is - // very annoying and usually very + // rather annoying and usually very // hard to fix. // // The solution to this problem @@ -602,29 +614,30 @@ class HelmholtzProblem // as long as it is derived from // the ``Subscriptor'' class. // - // In the present example - // program, we protect object - // using the pointer to the - // finite element, i.e. the - // following member variable, - // from the situation that for - // some reason the finite element - // pointed to is destroyed while - // still in use. Note that the - // pointer is assigned at - // construction time of this - // object, and destroyed upon - // destruction of this object, so - // the lock on the destruction of - // the finite element object is - // basically all through the - // lifetime of this object. + // In the present example program, we + // want to protect the finite element + // object from the situation that for + // some reason the finite element pointed + // to is destroyed while still in use. We + // therefore use a ``SmartPointer'' to + // the finite element object; since the + // finite element object is actually + // never changed in our computations, we + // pass a ``const FiniteElement'' as + // template argument to the + // ``SmartPointer'' class. Note that the + // pointer so declared is assigned at + // construction time of the solve object, + // and destroyed upon destruction, so the + // lock on the destruction of the finite + // element object extends throughout the + // lifetime of this ``HelmholtzProblem'' + // object. + Triangulation triangulation; + DoFHandler dof_handler; + SmartPointer > fe; - // The next few member variables - // are unspectacular, since they - // have already been discussed in - // detail: ConstraintMatrix hanging_node_constraints; SparsityPattern sparsity_pattern; @@ -633,7 +646,7 @@ class HelmholtzProblem Vector solution; Vector system_rhs; - // The second last variable + // The second to last variable // stores the refinement mode // passed to the // constructor. Since it is only @@ -646,30 +659,30 @@ class HelmholtzProblem // chance). const RefinementMode refinement_mode; - // For each refinement level some - // important data (like the - // number of cells, or the L2 - // error of the numerical - // solution) is printed out. The - // ``TableHandler'' can be used - // to collect all this data and - // to output it at the end of the - // run as a table in a simple - // text format or in Tex - // format. Here we don't only use - // the ``TableHandler'' but we - // use the derived class - // ``ConvergenceTable'' that - // additionally evaluates rates - // of convergence. + // For each refinement level some data + // (like the number of cells, or the L2 + // error of the numerical solution) will + // be generated and later printed. The + // ``TableHandler'' can be used to + // collect all this data and to output it + // at the end of the run as a table in a + // simple text or in LaTeX + // format. Here we don't only use the + // ``TableHandler'' but we use the + // derived class ``ConvergenceTable'' + // that additionally evaluates rates of + // convergence: ConvergenceTable convergence_table; }; + // @sect3{The ``HelmholtzProblem'' class implementation} + + // @sect4{HelmholtzProblem::HelmholtzProblem} // In the constructor of this class, // we only set the variables passed - // to this object, and associate the + // as arguments, and associate the // DoF handler object with the // triangulation (which is empty at // present, however). @@ -682,7 +695,9 @@ HelmholtzProblem::HelmholtzProblem (const FiniteElement &fe, {} + // @sect4{HelmholtzProblem::~HelmholtzProblem} + // This is no different than before: template HelmholtzProblem::~HelmholtzProblem () { @@ -690,46 +705,61 @@ HelmholtzProblem::~HelmholtzProblem () } - // The following function sets up the + // @sect4{HelmholtzProblem::setup_system} + + // The following function sets up the // degrees of freedom, sizes of // matrices and vectors, etc. Most of // its functionality has been showed // in previous examples, the only // difference being the renumbering - // step. + // step immediately after first + // distributing degrees of freedom. + // + // Renumbering the degrees of + // freedom is not overly difficult, + // as long as you use one of the + // algorithms included in the + // library. It requires only a single + // line of code. Some more information + // on this can be found in step-2. + // + // Note, however, that when you + // renumber the degrees of freedom, + // you must do so immediately after + // distributing them, since such + // things as hanging nodes, the + // sparsity pattern etc. depend on + // the absolute numbers which are + // altered by renumbering. + // + // The reason why we introduce renumbering + // here is that it is a relatively cheap + // operation but often has a beneficial + // effect: While the CG iteration itself is + // independent of the actual ordering of + // degrees of freedom, we will use SSOR as a + // preconditioner. SSOR goes through all + // degrees of freedom and does some + // operations that depend on what happened + // before; the SSOR operation is therefore + // not independent of the numbering of + // degrees of freedom, and it is known that + // its performance improves by using + // renumbering techniques. A little + // experiment shows that indeed, for example, + // the number of CG iterations for the fifth + // refinement cycle of adaptive refinement + // with the Q1 program used here is 40 + // without, but 36 with renumbering. Similar + // savings can generally be observed for all + // the computations in this program. template void HelmholtzProblem::setup_system () { dof_handler.distribute_dofs (*fe); - // Renumbering the degrees of - // freedom is not overly difficult, - // as long as you use one of the - // algorithms included in the - // library. It requires just one - // line of code, namely the - // following: DoFRenumbering::Cuthill_McKee (dof_handler); - // Note, however, that when you - // renumber the degrees of freedom, - // you must do so immediately after - // distributing them, since such - // things as hanging nodes, the - // sparsity pattern etc. depend on - // the absolute numbers which are - // altered by renumbering. - // - // Renumbering does not serve any - // specific purpose in this - // example, it is done only for - // exposition of the technique. To - // see the effect of renumbering on - // the sparsity pattern of the - // matrix, refer to the second - // example program. - - // The rest of the function is - // almost identically taken over - // from previous examples: + hanging_node_constraints.clear (); DoFTools::make_hanging_node_constraints (dof_handler, hanging_node_constraints); @@ -749,6 +779,7 @@ void HelmholtzProblem::setup_system () } + // @sect4{HelmholtzProblem::assemble_system} // Assembling the system of equations // for the problem at hand is mostly @@ -756,25 +787,34 @@ void HelmholtzProblem::setup_system () // before. However, some things have // changed anyway, so we comment on // this function fairly extensively. + // + // At the top of the function you will find + // the usual assortment of variable + // declarations. Compared to previous + // programs, of importance is only that we + // expect to solve problems also with + // bi-quadratic elements and therefore have + // to use sufficiently accurate quadrature + // formula. In addition, we need to compute + // integrals over faces, i.e. ``dim-1'' + // dimensional objects. The declaration of a + // face quadrature formula is then + // straightforward: template void HelmholtzProblem::assemble_system () { - // First we need to define objects - // which will be used as quadrature - // formula for domain and face - // integrals. - // - // Note the way in which we define - // a quadrature rule for the faces: - // it is simply a quadrature rule - // for one dimension less! QGauss quadrature_formula(3); QGauss face_quadrature_formula(3); - // For simpler use later on, we - // alias the number of quadrature - // points to local variables: + const unsigned int n_q_points = quadrature_formula.n_quadrature_points; const unsigned int n_face_q_points = face_quadrature_formula.n_quadrature_points; + + const unsigned int dofs_per_cell = fe->dofs_per_cell; + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); // Then we need objects which can // evaluate the values, gradients, @@ -791,11 +831,11 @@ void HelmholtzProblem::assemble_system () // requires the measure of the face // in a lower-dimensional // manifold. Internally these two - // classes are rooted on a common + // classes are rooted in a common // base class which does most of - // the work; that, however, is - // something that you need not - // worry about. + // the work and offers the same + // interface to both domain and + // interface integrals. // // For the domain integrals in the // bilinear form for Helmholtz's @@ -806,11 +846,11 @@ void HelmholtzProblem::assemble_system () // quadrature points on the real // cell (rather than on the unit // cell) to evaluate the right hand - // side function. - FEValues fe_values (*fe, quadrature_formula, - update_values | update_gradients | - update_q_points | update_JxW_values); - + // side function. The object we use + // to get at this information is + // the ``FEValues'' class discussed + // previously. + // // For the face integrals, we only // need the values of the shape // functions, as well as the @@ -819,19 +859,17 @@ void HelmholtzProblem::assemble_system () // the real cell since we want to // determine the Neumann values // from the exact solution object - // (see below). + // (see below). The class that gives + // us this information is called + // ``FEFaceValues'': + FEValues fe_values (*fe, quadrature_formula, + update_values | update_gradients | + update_q_points | update_JxW_values); + FEFaceValues fe_face_values (*fe, face_quadrature_formula, update_values | update_q_points | update_normal_vectors | update_JxW_values); - // In order to make programming - // more readable below, we alias - // the number of degrees of freedom - // per cell to a local variable, as - // already done for the number of - // quadrature points above: - const unsigned int dofs_per_cell = fe->dofs_per_cell; - // Then we need some objects // already known from previous // examples: An object denoting the @@ -841,38 +879,40 @@ void HelmholtzProblem::assemble_system () // right hand side, and the indices // of the degrees of freedom on a // cell. - RightHandSide right_hand_side; - std::vector rhs_values (n_q_points); - - FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); - Vector cell_rhs (dofs_per_cell); - - std::vector local_dof_indices (dofs_per_cell); - - // Then we define an object + // + // Note that the operations we will do with + // the right hand side object are only + // querying data, never changing the + // object. We can therefore declare it + // ``const'': + const RightHandSide right_hand_side; + std::vector rhs_values (n_q_points); + + // Finally we define an object // denoting the exact solution // function. We will use it to // compute the Neumann values at // the boundary from it. Usually, // one would of course do so using // a separate object, in particular - // since the exact solution is not - // known while the Neumann values + // since the exact solution is generally + // unknown while the Neumann values // are prescribed. We will, // however, be a little bit lazy // and use what we already have in // information. Real-life programs // would to go other ways here, of // course. - Solution exact_solution; + const Solution exact_solution; // Now for the main loop over all // cells. This is mostly unchanged // from previous examples, so we // only comment on the things that // have changed. - typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), - endc = dof_handler.end(); + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); for (; cell!=endc; ++cell) { cell_matrix = 0; @@ -880,7 +920,8 @@ void HelmholtzProblem::assemble_system () fe_values.reinit (cell); - right_hand_side.value_list (fe_values.get_quadrature_points(), rhs_values); + right_hand_side.value_list (fe_values.get_quadrature_points(), + rhs_values); for (unsigned int q_point=0; q_point::assemble_system () // now contains the // additional term from // the Helmholtz - // equation, namely the - // scalar products of - // the two function - // values, rather than - // their gradients, - // which is the second - // term below: + // equation: cell_matrix(i,j) += ((fe_values.shape_grad(i,q_point) * - fe_values.shape_grad(j,q_point) * - fe_values.JxW(q_point)) - + - (fe_values.shape_value(i,q_point) * - fe_values.shape_value(j,q_point) * - fe_values.JxW(q_point))); + fe_values.shape_grad(j,q_point) + + + fe_values.shape_value(i,q_point) * + fe_values.shape_value(j,q_point)) * + fe_values.JxW(q_point)); cell_rhs(i) += (fe_values.shape_value(i,q_point) * rhs_values [q_point] * fe_values.JxW(q_point)); - }; + } // Then there is that second // term on the right hand side, // the contour integral. First // we have to find out whether - // the intersection of the face + // the intersection of the faces // of this cell with the // boundary part Gamma2 is // nonzero. To this end, we @@ -926,20 +960,18 @@ void HelmholtzProblem::assemble_system () // which is the value that we // have assigned to that // portions of the boundary - // composing Gamma2 in a - // function further below. The + // composing Gamma2 in the + // ``run()'' function further + // below. (The // default value of boundary - // indicators is ``0'' for - // external faces, and ``255'' - // for internal faces (the - // latter value should never be - // changed, and there is also - // no need to do so), so faces + // indicators is ``0'', so faces // can only have an indicator // equal to ``1'' if we have - // explicitly set it. + // explicitly set it.) for (unsigned int face=0; face::faces_per_cell; ++face) - if (cell->face(face)->boundary_indicator() == 1) + if (cell->face(face)->at_boundary() + && + (cell->face(face)->boundary_indicator() == 1)) { // If we came into here, // then we have found an @@ -965,42 +997,39 @@ void HelmholtzProblem::assemble_system () // integration by using a // loop over all // quadrature points. + // + // On each quadrature point, we + // first compute the value of the + // normal derivative. We do so + // using the gradient of the + // exact solution and the normal + // vector to the face at the + // present quadrature point + // obtained from the + // ``fe_face_values'' + // object. This is then used to + // compute the additional + // contribution of this face to + // the right hand side: for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); for (unsigned int i=0; i::assemble_system () cell_matrix(i,j)); system_rhs(local_dof_indices[i]) += cell_rhs(i); - }; - }; - - // The rest of the function has - // also been shown previously: - hanging_node_constraints.condense (system_matrix); - hanging_node_constraints.condense (system_rhs); + } + } - // Only with the interpolation of - // boundary values, there is one - // notable thing, namely that now + // Likewise, elimination and treatment of + // boundary values has been shown + // previously. + // + // We note, however that now // the boundary indicator for which // we interpolate boundary values // (denoted by the second parameter // to // ``interpolate_boundary_values'') // does not represent the whole - // boundary an more. Rather, it is + // boundary any more. Rather, it is // that portion of the boundary // which we have not assigned // another indicator (see @@ -1035,7 +1061,11 @@ void HelmholtzProblem::assemble_system () // at the boundary that do not // belong to Gamma1 are therefore // excluded from the interpolation - // of boundary values. + // of boundary values, just as + // we want. + hanging_node_constraints.condense (system_matrix); + hanging_node_constraints.condense (system_rhs); + std::map boundary_values; VectorTools::interpolate_boundary_values (dof_handler, 0, @@ -1048,8 +1078,10 @@ void HelmholtzProblem::assemble_system () } - // Solving the system of equations is - // done in the same way as before. + // @sect4{HelmholtzProblem::solve} + + // Solving the system of equations is + // done in the same way as before: template void HelmholtzProblem::solve () { @@ -1066,44 +1098,79 @@ void HelmholtzProblem::solve () } - // Now for the function doing grid + // @sect4{HelmholtzProblem::refine_grid} + + // Now for the function doing grid // refinement. Depending on the // refinement mode passed to the // constructor, we do global or // adaptive refinement. + // + // Global refinement is simple, + // so there is + // not much to comment on. + // In case of adaptive + // refinement, we use the same + // functions and classes as in + // the previous example + // program. Note that one + // could treat Neumann + // boundaries differently than + // Dirichlet boundaries, and + // one should in fact do so + // here since we have Neumann + // boundary conditions on part + // of the boundaries, but + // since we don't have a + // function here that + // describes the Neumann + // values (we only construct + // these values from the exact + // solution when assembling + // the matrix), we omit this + // detail even though they would + // not be hard to add. + // + // At the end of the switch, we have a + // default case that looks slightly strange: + // an ``Assert'' statement with a ``false'' + // condition. Since the ``Assert'' macro + // raises an error whenever the condition is + // false, this means that whenever we hit + // this statement the program will be + // aborted. This in intentional: Right now we + // have only implemented two refinement + // strategies (global and adaptive), but + // someone might want to add a third strategy + // (for example adaptivity with a different + // refinement criterion) and add a third + // member to the enumeration that determines + // the refinement mode. If it weren't for the + // default case of the switch statement, this + // function would simply run to its end + // without doing anything. This is most + // likely not what was intended. One of the + // defensive programming techniques that you + // will find all over the deal.II library is + // therefore to always have default cases + // that abort, to make sure that values not + // considered when listing the cases in the + // switch statement are eventually caught, + // and forcing programmers to add code to + // handle them. We will use this same + // technique in other places further down as + // well. template void HelmholtzProblem::refine_grid () { switch (refinement_mode) { - // If global refinement is - // required, this is simple: case global_refinement: { triangulation.refine_global (1); break; - }; - - // In case of adaptive - // refinement, we use the same - // functions and classes as in - // the previous example - // program. Note that one - // could treat Neumann - // boundaries differently than - // Dirichlet boundaries, and - // one should in fact do so - // here since we have Neumann - // boundary conditions on part - // of the boundaries, but - // since we don't have a - // function here that - // describes the Neumann - // values (we only construct - // these values from the exact - // solution when assembling - // the matrix), we omit this - // detail here. + } + case adaptive_refinement: { Vector estimated_error_per_cell (triangulation.n_active_cells()); @@ -1122,24 +1189,30 @@ void HelmholtzProblem::refine_grid () triangulation.execute_coarsening_and_refinement (); break; - }; - }; + } + + default: + { + Assert (false, ExcNotImplemented()); + } + } } + // @sect4{HelmholtzProblem::process_solution} - // Finally process the solution after - // it has been computed. For this, we - // integrate the error in various - // norms, and we generate tables that - // will be later used to display the - // convergence against the continuous - // solution in a nice format. + // Finally we want to process the solution + // after it has been computed. For this, we + // integrate the error in various norms, and + // we generate tables that will later be used + // to display the convergence against the + // continuous solution in a nice format. template void HelmholtzProblem::process_solution (const unsigned int cycle) { - // In order to integrate the - // difference between computed + // Our first task is to compute + // error norms. In order to integrate + // the difference between computed // numerical solution and the // continuous solution (described // by the ``Solution'' class @@ -1152,12 +1225,11 @@ void HelmholtzProblem::process_solution (const unsigned int cycle) // quantities, we save some memory // by using ``float'' instead of // ``double'' values. - Vector difference_per_cell (triangulation.n_active_cells()); - - // Next we use a function from the - // library which computes the error - // in the L2 norm on each cell. We - // have to pass it the DoF handler + // + // The next step is to use a function + // from the library which computes the + // error in the L2 norm on each cell. + // We have to pass it the DoF handler // object, the vector holding the // nodal values of the numerical // solution, the continuous @@ -1171,12 +1243,7 @@ void HelmholtzProblem::process_solution (const unsigned int cycle) // formula with three points in // each space direction, and // compute the L2 norm. - VectorTools::integrate_difference (dof_handler, - solution, - Solution(), - difference_per_cell, - QGauss(3), - VectorTools::L2_norm); + // // Finally, we want to get the // global L2 norm. This can of // course be obtained by summing @@ -1185,11 +1252,21 @@ void HelmholtzProblem::process_solution (const unsigned int cycle) // of that value. This is // equivalent to taking the l2 // (lower case ``l'') norm of the - // vector of norms on each cell: + // vector of norms on each cell: + Vector difference_per_cell (triangulation.n_active_cells()); + VectorTools::integrate_difference (dof_handler, + solution, + Solution(), + difference_per_cell, + QGauss(3), + VectorTools::L2_norm); const double L2_error = difference_per_cell.l2_norm(); - // The same procedure is done to - // get the H1 semi-norm: + // By same procedure we get the H1 + // semi-norm. We re-use the + // ``difference_per_cell'' vector since it + // is no longer used after computing the + // ``L2_error'' variable above. VectorTools::integrate_difference (dof_handler, solution, Solution(), @@ -1200,10 +1277,10 @@ void HelmholtzProblem::process_solution (const unsigned int cycle) // Finally, we compute the maximum // norm. Of course, we can't - // actually use the true maximum, + // actually compute the true maximum, // but only the maximum at the // quadrature points. Since this - // quite sensitively depends on the + // depends quite sensitively on the // quadrature rule being used, and // since we would like to avoid // false results due to @@ -1219,27 +1296,40 @@ void HelmholtzProblem::process_solution (const unsigned int cycle) // that tells it how often it shall // use this rule in each space // direction. - QTrapez<1> q_trapez; - QIterated q_iterated (q_trapez, 5); - - // Using this special quadrature - // rule, we can now try to find the - // maximal error on each cell: + // + // Using this special quadrature rule, we + // can then try to find the maximal error + // on each cell. Finally, we compute the + // global L infinity error from the L + // infinite errors on each cell. Instead of + // summing squares, we now have to take the + // maximum value over all cell-wise + // entries, an operation that is + // conveniently done using the + // ``Vector::linfty'' function: + const QTrapez<1> q_trapez; + const QIterated q_iterated (q_trapez, 5); VectorTools::integrate_difference (dof_handler, solution, Solution(), difference_per_cell, q_iterated, VectorTools::Linfty_norm); - // Obviously, the maximal error - // globally is the maximum over the - // maximal errors on each cell: const double Linfty_error = difference_per_cell.linfty_norm(); // After all these errors have been // computed, we finally write some - // output and put all the data into - // a table. + // output. In addition, we add the + // important data to the + // ``TableHandler'' by specifying + // the key of the column and the value. + // Note that it is not necessary to + // define column keys beforehand -- it is + // sufficient to just add values, + // and columns will be + // introduced into the table in the + // order values are added the + // first time. const unsigned int n_active_cells=triangulation.n_active_cells(); const unsigned int n_dofs=dof_handler.n_dofs(); @@ -1252,180 +1342,109 @@ void HelmholtzProblem::process_solution (const unsigned int cycle) << n_dofs << std::endl; - // Add the important data to the - // ``TableHandler'' by giving the key - // of the column and the value. - // You don't need to define the keys - // beforehand, just add the values, - // and the column will be - // introduced into the table in the - // order the values are added the - // first time. convergence_table.add_value("cycle", cycle); convergence_table.add_value("cells", n_active_cells); convergence_table.add_value("dofs", n_dofs); convergence_table.add_value("L2", L2_error); convergence_table.add_value("H1", H1_error); convergence_table.add_value("Linfty", Linfty_error); - // You may set the precision with - // which the values will be written - // upon output. - convergence_table.set_precision("L2", 3); - convergence_table.set_precision("H1", 3); - convergence_table.set_precision("Linfty", 3); - // The default notation is fixed - // point. For the columns you'd - // like to see in scientific notation - // set the `scientific_flag' `true' - // by the following lines: - convergence_table.set_scientific("L2", true); - convergence_table.set_scientific("H1", true); - convergence_table.set_scientific("Linfty", true); - // For the output of a table into a - // LaTeX file, the default captions of - // the columns are the keys given - // as argument to the ``add_value'' - // functions. If you'd like to - // have TeX captions that differ - // from the default ones you can - // specify them by the following. - convergence_table.set_tex_caption("cells", "\\# cells"); - convergence_table.set_tex_caption("dofs", "\\# dofs"); - convergence_table.set_tex_caption("L2", "$L^2$-error"); - convergence_table.set_tex_caption("H1", "$H^1$-error"); - convergence_table.set_tex_caption("Linfty", "$L^\\infty$-error"); - // Note, that `\\' is reduced to - // `\' by the compiler such that the - // real TeX caption is e.g. - // `$L^\infty$-error'. - // - // The default TeX format of each - // column of the table is `c' - // (centered). To specify a - // different (e.g. `right') one, - // the following function may be - // used: - convergence_table.set_tex_format("cells", "r"); - convergence_table.set_tex_format("dofs", "r"); } + // @sect4{HelmholtzProblem::run} - // The following function is the main - // one which controls the flow of - // execution. The basic layout is as - // in previous examples: an outer - // loop over successively refined - // grids, and in this loop first - // problem setup, assemblage of the + // As in previous example programs, the + // ``run'' function controls controls the + // flow of execution. The basic layout is as + // in previous examples: an outer loop over + // successively refined grids, and in this + // loop first problem setup, assembling the // linear system, solution, and // post-processing. + // + // The first task in the main loop is + // creation and refinement of grids. This is + // as in previous examples, with the only + // difference that we want to have part of + // the boundary marked as Neumann type, + // rather than Dirichlet. + // + // For this, we will use the following + // convention: Faces belonging to Gamma1 will + // have the boundary indicator ``0'' (which + // is the default, so we don't have to set it + // explicitely), and faces belonging to + // Gamma2 will use ``1'' as boundary + // indicator. To set these values, we loop + // over all cells, then over all faces of a + // given cell, check whether it is part of + // the boundary that we want to denote by + // Gamma2, and if so set its boundary + // indicator to ``1''. For the present + // program, we consider the left and bottom + // boundaries as Gamma2. We determine whether + // a face is part of that boundary by asking + // whether the x or y coordinates + // (i.e. vector components 0 and 1) of the + // midpoint of a face equals -1. + // + // It is worth noting that + // we have to loop over all + // cells here, not only the + // active ones. The reason + // is that upon refinement, + // newly created faces + // inherit the boundary + // indicator of their + // parent face. If we now + // only set the boundary + // indicator for active + // faces, coarsen some + // cells and refine them + // later on, they will + // again have the boundary + // indicator of the parent + // cell which we have not + // modified, instead of the + // one we + // intended. Consequently, we + // have to change the + // boundary indicators of + // faces of all cells on Gamma2, + // whether they are active or not. + // Alternatively, we could of + // course have done this job on + // the coarsest mesh (i.e. before + // the first refinement step) and + // refined the mesh only after that. template void HelmholtzProblem::run () { for (unsigned int cycle=0; cycle<7; ++cycle) { - // The first action in each - // iteration of the outer loop - // is setting up the grid on - // which we will solve in this - // iteration. In the first - // iteration, the coarsest grid - // is generated, in later - // iterations it is refined, - // for which we call the - // ``refine_grid'' function. if (cycle == 0) { - // Setting up the coarse - // grid is done as in - // previous examples: we - // first create an initial - // grid, which is the unit - // square [-1,1]x[-1,1] in - // the present case. Then - // we refine it globally a - // specific number of - // times. GridGenerator::hyper_cube (triangulation, -1, 1); triangulation.refine_global (1); - // However, here we have to - // do something else in - // addition: mark those - // faces that belong to the - // different components of - // the boundary, Gamma1 and - // Gamma2. We will use the - // following convention: - // Faces belonging to - // Gamma1 will have the - // boundary indicator ``0'' - // (which is the default, - // so we don't have to set - // it explicitely), and - // faces belonging to - // Gamma2 will use ``1'' as - // boundary indicator. - // - // To set these values, we - // loop over all cells, - // then over all faces of a - // given cell, check - // whether it belongs to - // the boundary Gamma2, and - // if so set its boundary - // indicator to ``1''. - // - // It is worth noting that - // we have to loop over all - // cells here, not only the - // active ones. The reason - // is that upon refinement, - // newly created faces - // inherit the boundary - // indicator of their - // parent face. If we now - // only set the boundary - // indicator for active - // faces, coarsen some - // cells and refine them - // later on, they will - // again have the boundary - // indicator of the parent - // cell which we have not - // modified, instead of the - // one we - // intended. Therefore, we - // have to change the - // boundary indicators of - // all faces on Gamma2, - // irrespective whether - // they are active or not. - typename Triangulation::cell_iterator cell = triangulation.begin (), - endc = triangulation.end(); + typename Triangulation::cell_iterator + cell = triangulation.begin (), + endc = triangulation.end(); for (; cell!=endc; ++cell) - for (unsigned int face=0; face::faces_per_cell; ++face) + for (unsigned int face=0; + face::faces_per_cell; + ++face) if ((cell->face(face)->center()(0) == -1) || (cell->face(face)->center()(1) == -1)) cell->face(face)->set_boundary_indicator (1); } else - { - // If this is not the first - // step, the we call - // ``refine_grid'' to - // actually refine the grid - // according to the - // refinement mode passed to - // the constructor. - refine_grid (); - }; + refine_grid (); - // The next steps you already - // know from previous + // The next steps are already + // known from previous // examples. This is mostly the // basic set-up of every finite // element program: @@ -1436,100 +1455,87 @@ void HelmholtzProblem::run () // The last step in this chain // of function calls is usually - // evaluation of the computed + // the evaluation of the computed // solution for the quantities // one is interested in. This // is done in the following - // function. We pass the number - // of the loop iteration since - // that might be of interest to - // see in the logs which this - // function produces. + // function. Since the function + // generates output that indicates + // the number of the present + // refinement step, we pass this + // number as an argument. process_solution (cycle); - }; + } - // After the last iteration we - // output the solution on the - // finest grid. This is done using - // the following sequence of - // statements which you have - // already seen in previous - // examples: - std::string filename; + // After the last iteration we output the + // solution on the finest grid. This is + // done using the following sequence of + // statements which we have already + // discussed in previous examples. The + // first step is to generate a suitable + // filename (called ``gmv_filename'' here, + // since we want to output data in GMV + // format; we add the prefix to distinguish + // the filename from that used for other + // output files further down below). Here, + // we augment the name by the mesh + // refinement algorithm, and as above we + // make sure that we abort the program if + // another refinement method is added and + // not handled by the following switch + // statement: + std::string gmv_filename; switch (refinement_mode) { case global_refinement: - filename = "solution-global"; + gmv_filename = "solution-global"; break; case adaptive_refinement: - filename = "solution-adaptive"; + gmv_filename = "solution-adaptive"; break; default: - Assert (false, ExcInternalError()); - }; + Assert (false, ExcNotImplemented()); + } - // We augment the filename by a - // postfix denoting the finite - // element which we have used in - // the computation. Finding out - // which finite element we are - // actually using is not that - // simple here, since we only have - // a pointer to the common base - // class of all finite elements, - // which does not know anything - // about polynomial - // degrees. However, we actually - // know that we have generated a - // finite element of class - // ``FE_Q'', so we can use some C++ - // feature to actually get a - // reference to the ``FE_Q'' - // element pointed to by the - // reference and ask it for the - // polynomial degree. Note that if - // for whatever reason the object - // referenced behind the pointer to - // the base class should not be of - // type ``FE_Q'', then the C++ - // language lets the - // ``dynamic_cast'' operator - // applied to a reference type - // throw an exception (if it were a - // pointer type, then a null - // pointer would be returned, which - // would then yield a segmentation - // fault when dereferenced in the - // subsequent call to - // ``get_order''). - switch (dynamic_cast&>(*fe).get_degree()) + // We augment the filename by a postfix + // denoting the finite element which we + // have used in the computation. To this + // end, the finite element base class + // stores the maximal polynomial degree of + // shape functions in each coordinate + // variable as a variable ``degree'', and + // we use for the switch statement (note + // that the polynomial degree of bilinear + // shape functions is really 2, since they + // contain the term ``x*y''; however, the + // polynomial degree in each coordinate + // variable is still only 1). We again use + // the same defensive programming technique + // to safeguard against the case that the + // polynomial degree has an unexpected + // value, using the ``Assert (false, + // ExcNotImplemented())'' idiom in the + // default branch of the switch statement: + switch (fe->degree) { case 1: - filename += "-q1"; + gmv_filename += "-q1"; break; case 2: - filename += "-q2"; + gmv_filename += "-q2"; break; default: - // The finite element is - // neither Q1 nor Q2. This - // should not have happened, - // but maybe someone has tried - // to change this in ``main'', - // so it might happen. We catch - // this case and throw an - // exception, since we don't - // know how to name the - // respective output file - Assert (false, ExcInternalError()); - }; - - - filename += ".gmv"; - - std::ofstream output (filename.c_str()); + Assert (false, ExcNotImplemented()); + } + // Once we have the base name for the + // output file, we add an extension + // appropriate for GMV output, open a file, + // and add the solution vector to the + // object that will do the actual output: + gmv_filename += ".gmv"; + std::ofstream output (gmv_filename.c_str()); DataOut data_out; data_out.attach_dof_handler (dof_handler); @@ -1547,8 +1553,8 @@ void HelmholtzProblem::run () // all output formats only support // bilinear data, the data is // written only bilinear, and - // information is lost - // therefore. Of course, we can't + // information is consequently lost. + // Of course, we can't // change the format in which // graphic programs accept their // inputs, but we can write the @@ -1590,67 +1596,120 @@ void HelmholtzProblem::run () // ``q'' subdivisions, and the // order of the elements is // determined in the same way as - // above in the ``run'' function: - const unsigned int - n_subcells = dynamic_cast&>(*fe).get_degree(); - data_out.build_patches (n_subcells); - - // Finally write out the data in - // GMV format. + // above. + // + // With the intermediate format + // so generated, we can then actually + // write the graphical output in GMV + // format: + data_out.build_patches (fe->degree); data_out.write_gmv (output); - // In each cycle values were added - // to the TableHandler. Now write - // the table to the standard output - // stream `std::cout'. Note, that the - // output in text format is a quite - // simple one and the captions may - // not be printed directly above - // the specific columns. + // After graphical output, we would also + // like to generate tables from the error + // computations we have done in + // ``process_solution''. There, we have + // filled a table object with the number of + // cells for each refinement step as well + // as the errors in different norms. + + // For a nicer textual output of this data, + // one may want to set the precision with + // which the values will be written upon + // output. We use 3 digits for this, which + // is usually sufficient for error + // norms. By default, data is written in + // fixed point notation. However, for + // columns one would like to see in + // scientific notation another function + // call sets the `scientific_flag' to + // `true', leading to floating point + // representation of numbers. + convergence_table.set_precision("L2", 3); + convergence_table.set_precision("H1", 3); + convergence_table.set_precision("Linfty", 3); + + convergence_table.set_scientific("L2", true); + convergence_table.set_scientific("H1", true); + convergence_table.set_scientific("Linfty", true); + + // For the output of a table into a LaTeX + // file, the default captions of the + // columns are the keys given as argument + // to the ``add_value'' functions. To have + // TeX captions that differ from the + // default ones you can specify them by the + // following function calls. + // Note, that `\\' is reduced to + // `\' by the compiler such that the + // real TeX caption is, e.g., + // `$L^\infty$-error'. + convergence_table.set_tex_caption("cells", "\\# cells"); + convergence_table.set_tex_caption("dofs", "\\# dofs"); + convergence_table.set_tex_caption("L2", "$L^2$-error"); + convergence_table.set_tex_caption("H1", "$H^1$-error"); + convergence_table.set_tex_caption("Linfty", "$L^\\infty$-error"); + + // Finally, the default LaTeX format for + // each column of the table is `c' + // (centered). To specify a different + // (e.g. `right') one, the following + // function may be used: + convergence_table.set_tex_format("cells", "r"); + convergence_table.set_tex_format("dofs", "r"); + + // After this, we can finally write the + // table to the standard output stream + // ``std::cout'' (after one extra empty + // line, to make things look + // prettier). Note, that the output in text + // format is quite simple and that + // captions may not be printed directly + // above the specific columns. + std::cout << std::endl; convergence_table.write_text(std::cout); - // The table can also be written - // into a Tex file. The (nicely) - // formatted table can be viewed at - // after calling `latex filename' - // and e.g. `xdvi filename', where - // filename is the name of the file - // which we construct from the name - // of the finite element and the - // refinement mode, as above - if (true) + + // The table can also be written into a + // LaTeX file. The (nicely) formatted + // table can be viewed at after calling + // `latex filename' and e.g. `xdvi + // filename', where filename is the name of + // the file to which we will write output + // now. We construct its name in the same + // way as before, but with a different + // prefix "error": + std::string error_filename = "error"; + switch (refinement_mode) { - std::string filename = "error"; - switch (refinement_mode) - { - case global_refinement: - filename += "-global"; - break; - case adaptive_refinement: - filename += "-adaptive"; - break; - default: - Assert (false, ExcInternalError()); - }; + case global_refinement: + error_filename += "-global"; + break; + case adaptive_refinement: + error_filename += "-adaptive"; + break; + default: + Assert (false, ExcNotImplemented()); + } - switch (dynamic_cast&>(*fe).get_degree()) - { - case 1: - filename += "-q1"; - break; - case 2: - filename += "-q2"; - break; - default: - Assert (false, ExcInternalError()); - }; - - filename += ".tex"; - - std::ofstream table_file(filename.c_str()); - convergence_table.write_tex(table_file); - table_file.close(); + switch (fe->degree) + { + case 1: + error_filename += "-q1"; + break; + case 2: + error_filename += "-q2"; + break; + default: + Assert (false, ExcNotImplemented()); } - // In case you want the same + + error_filename += ".tex"; + std::ofstream error_table_file(error_filename.c_str()); + + convergence_table.write_tex(error_table_file); + + + // In case you want the same // caption for several columns, you // can merge some columns to a // super column by @@ -1731,9 +1790,9 @@ void HelmholtzProblem::run () filename += "-adaptive"; break; default: - Assert (false, ExcInternalError()); - }; - switch (dynamic_cast&>(*fe).get_degree()) + Assert (false, ExcNotImplemented()); + } + switch (fe->degree) { case 1: filename += "-q1"; @@ -1742,8 +1801,8 @@ void HelmholtzProblem::run () filename += "-q2"; break; default: - Assert (false, ExcInternalError()); - }; + Assert (false, ExcNotImplemented()); + } filename += ".tex"; std::ofstream table_file(filename.c_str()); @@ -1802,7 +1861,7 @@ int main () helmholtz_problem_2d.run (); std::cout << std::endl; - }; + } { std::cout << "Solving with Q1 elements, global refinement" << std::endl @@ -1814,7 +1873,7 @@ int main () helmholtz_problem_2d.run (); std::cout << std::endl; - }; + } { std::cout << "Solving with Q2 elements, global refinement" << std::endl @@ -1826,7 +1885,7 @@ int main () helmholtz_problem_2d.run (); std::cout << std::endl; - }; + } } catch (std::exception &exc) @@ -1851,7 +1910,7 @@ int main () << "----------------------------------------------------" << std::endl; return 1; - }; + } return 0; }