From f92a144b969033b57c6a62a6ae60511ba06281b8 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Thu, 15 Apr 2004 20:19:57 +0000 Subject: [PATCH] Some documentation. git-svn-id: https://svn.dealii.org/trunk@9017 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-15/step-15.cc | 290 +++++++++++++++++++++------- 1 file changed, 222 insertions(+), 68 deletions(-) diff --git a/deal.II/examples/step-15/step-15.cc b/deal.II/examples/step-15/step-15.cc index 816c578ca9..fe12d4764c 100644 --- a/deal.II/examples/step-15/step-15.cc +++ b/deal.II/examples/step-15/step-15.cc @@ -11,38 +11,42 @@ /* to the file deal.II/doc/license.html for the text and */ /* further information on this license. */ + // As usual, most of the headers here have + // already been used and discussed in + // previous examples: #include #include #include #include #include #include - // #include #include #include #include -#include #include +#include #include #include #include +#include #include +#include #include +#include #include #include #include #include -#include -#include -#include - -#include -#include - -#include + // This is probably the only new one - it + // declares the class that we use to transfer + // a solution on one grid to the one we + // obtain after refining/coarsening it: +#include + // And here comes the usual assortment of C++ + // header files: #include #include @@ -53,27 +57,164 @@ #endif + // The first thing we have here is a helper + // function that computes an even power |v|^n + // of a vector ``v'', by evaluating + // (v*v)^(n/2). We need this in the + // computations below where we do not want to + // dwell on the fact that the gradient of the + // solution is actually a scalar in the 1d + // situation we consider in this program (in + // 1d, the gradient is a vector with a single + // element, which is easily extracted). Small + // tricks like this make it significantly + // simpler to later extend a program so that + // it also runs in higher space dimensions. + // + // While the implementation of the function + // is obvious, note the assertion at the + // beginning of the function body, which + // makes sure that the exponent is indeed an + // even number (here, we use that ``n/2'' is + // computed in integer arithmetic, i.e. any + // remainder of the division is + // lost). ``ExcMessage'' is a pre-defined + // exception class that takes a string + // argument explaining what goes wrong. It is + // a simpler way to declare exceptions than + // the ones shown in step-9 and step-13/14 + // where we explicitly declared exception + // classes. However, by using a generic + // exception class, we lose the ability to + // attach additional information at run-time + // to the exception message, such as the + // value of the variable ``n''. By following + // the way explained in above example + // programs, adding this feature is simple, + // though. +template +inline +double gradient_power (const Tensor<1,dim> &v, + const unsigned int n) +{ + Assert ((n/2)*2 == n, ExcMessage ("Value of 'n' must be even")); + double p = 1; + for (unsigned int k=0; k +{ + public: + InitializationValues () : Function<1>() {}; + + virtual double value (const Point<1> &p, + const unsigned int component = 0) const; +}; + + + // So here comes the function that implements + // the function object. The ``base'' value is + // ``x^1/3'', while ``random'' is a random + // number between -1 and 1 (note that + // ``rand()'' returns a random integer value + // between zero and ``RAND_MAX''; to convert + // it to a floating point value between 0 and + // 2, we have to divide by ``RAND_MAX'' and + // multiply by two -- note that the first + // multiplication has to happen in floating + // point arithmetic, so that the division is + // done in non-truncating floating point mode + // as well; the final step is then to shift + // the interval [0,2] to [-1,1]). + // + // In a second step, we add the base value + // and a random value in [-0.1,0.1] together + // and return it, unless it is less than + // zero, in which case we take zero. +double InitializationValues::value (const Point<1> &p, + const unsigned int) const +{ + const double base = std::pow(p(0), 1./3.); + const double random = 2.*rand()/RAND_MAX-1; + return std::max (base+.1*random, 0.); +} + + + + // Next is the declaration of the main + // class. As in most of the previous example + // programs, the public interface of the + // class consists only of a constructor and a + // ``run'' function that does the actual + // work. The constructor takes an additional + // argument that indicates the number of the + // run we are presently performing. This + // value is only used at the very end when we + // generate graphical output with a filename + // that matches this number. + // + // The private section of the class has the + // usual assortment of functions setting up + // the computations, doing one nonlinear + // step, refineming the mesh, doing a line + // search for step length computations, + // etc. The ``energy'' function computes the + // value of the optimization functional on an + // arbitrary finite element function with + // nodal values given on the ``DoFHandler'' + // given as an argument. Since it does not + // depend on the state of this object, we + // declare this function as ``static''. + // + // The member variables of this class are + // what we have seen before, and the + // variables that characterize the linear + // system to be solved in the next nonlinear + // step, as well as the present approximation + // of the solution. template class MinimizationProblem { public: - MinimizationProblem (); - ~MinimizationProblem (); + MinimizationProblem (const unsigned int run_number); void run (); - void output_results (const unsigned int cycle) const; private: - void setup_system (); + void initialize (); + void setup_system_on_mesh (); void assemble_step (); double line_search (const Vector & update) const; - void do_step (); - void initialize (); + void output_results () const; void refine_grid (); + void do_step (); static double energy (const DoFHandler &dof_handler, const Vector &function); + + const unsigned int run_number; Triangulation triangulation; @@ -91,48 +232,38 @@ class MinimizationProblem -class InitializationValues : public Function<1> -{ - public: - InitializationValues () : Function<1>() {}; - - virtual double value (const Point<1> &p, - const unsigned int component = 0) const; -}; - - - -double InitializationValues::value (const Point<1> &p, - const unsigned int) const -{ - const double base = std::pow(p(0), 1./3.); - const double random = 2.*rand()/RAND_MAX-1; - if (base+.1*random < 0 ) - return 0; - else - return base+.1*random; -} - - - + // The constructor of this class is actually + // somewhat boring: template -MinimizationProblem::MinimizationProblem () +MinimizationProblem::MinimizationProblem (const unsigned int run_number) : + run_number (run_number), fe (1), dof_handler (triangulation) {} + // And so is the function that prepares the + // member variables of this class for + // assembling the linear system in each + // nonlinear step. This has all been shown + // before in previous example programs. Note, + // however, that all this works in 1d just as + // in any other space dimension, and would + // not require any changes if we were to use + // the program in another space dimension. + // + // Note that this function is only called + // when the mesh has been changed (or before + // the first nonlinear step). It only + // initializes the variables to their right + // sizes, but since these sizes don't change + // as long as we don't change the mesh, we + // can use them for more than just one + // nonlinear iteration without reinitializing + // them. template -MinimizationProblem::~MinimizationProblem () -{ - dof_handler.clear (); -} - - - -template -void MinimizationProblem::setup_system () +void MinimizationProblem::setup_system_on_mesh () { hanging_node_constraints.clear (); DoFTools::make_hanging_node_constraints (dof_handler, @@ -150,26 +281,45 @@ void MinimizationProblem::setup_system () } -template -double gradient_power (const Tensor<1,dim> &v, - const unsigned int n) -{ - Assert ((n/2)*2 == n, ExcMessage ("Value of 'n' must be even")); - double p = 1; - for (unsigned int k=0; k void MinimizationProblem::assemble_step () { + // The first two lines of the function + // clear the matrix and right hand side + // values of their prior content, which + // could possibly still be there from the + // previous nonlinear step. matrix.reinit (sparsity_pattern); residual.reinit (dof_handler.n_dofs()); - - QGauss3 quadrature_formula; + // Then we initialize a ``FEValues'' object + // with a 3-point Gauss quadrature + // formula. This object will be used to + // compute the values and gradients of the + // shape functions at the quadrature + // points, which we need to assemble the + // matrix and right hand side of the + // nonlinear step as outlined in the + // introduction to this example program. In + // order to compute values and gradients, + // we need to pass the ``update_values'' + // and ``update_gradients'' flags to the + // constructor, and the + // ``update_JxW_values'' flag for the + // Jacobian times the weight at a + // quadrature point. In addition, we need + // to have the coordinate values of each + // quadrature point in real space for the + // ``x-u^3'' terms; to get these from the + // ``FEValues'' object, we need to pass it + // the ``update_q_points'' flag. + QGauss3 quadrature_formula; FEValues fe_values (fe, quadrature_formula, UpdateFlags(update_values | update_gradients | @@ -351,6 +501,7 @@ void MinimizationProblem::do_step () } + template <> void MinimizationProblem<1>::initialize () { @@ -545,7 +696,8 @@ MinimizationProblem::energy (const DoFHandler &dof_handler, template -void MinimizationProblem::output_results (const unsigned int cycle) const +void +MinimizationProblem::output_results () const { DataOut data_out; data_out.attach_dof_handler (dof_handler); @@ -558,7 +710,7 @@ void MinimizationProblem::output_results (const unsigned int cycle) const std::ostrstream filename; #endif filename << "solution-" - << cycle + << run_number << ".gnuplot" << std::ends; #ifdef HAVE_STD_STRINGSTREAM @@ -583,7 +735,7 @@ void MinimizationProblem::run () while (true) { - setup_system (); + setup_system_on_mesh (); unsigned int iteration=0; for (; iteration<5; ++iteration) @@ -599,6 +751,9 @@ void MinimizationProblem::run () refine_grid (); } + + output_results (); + std::cout << std::endl; } @@ -614,9 +769,8 @@ int main () { std::cout << "Realization " << realization << ":" << std::endl; - MinimizationProblem<1> minimization_problem_1d; + MinimizationProblem<1> minimization_problem_1d (realization); minimization_problem_1d.run (); - minimization_problem_1d.output_results (realization); } } catch (std::exception &exc) -- 2.39.5