From e21a8ec8f20212d8987a384e55eb4ccfd7867ed0 Mon Sep 17 00:00:00 2001 From: frohne Date: Wed, 9 Nov 2011 00:22:32 +0000 Subject: [PATCH] Change the commentaries git-svn-id: https://svn.dealii.org/trunk@24739 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-41/step-41.cc | 484 +++++----------------------- 1 file changed, 83 insertions(+), 401 deletions(-) diff --git a/deal.II/examples/step-41/step-41.cc b/deal.II/examples/step-41/step-41.cc index cb8c34fb5e..c44c5f68ba 100644 --- a/deal.II/examples/step-41/step-41.cc +++ b/deal.II/examples/step-41/step-41.cc @@ -12,11 +12,6 @@ // @sect3{Include files} - // The first few (many?) include - // files have already been used in - // the previous example, so we will - // not explain their meaning here - // again. #include #include #include @@ -48,44 +43,21 @@ #include #include - // This is new, however: in the previous - // example we got some unwanted output from - // the linear solvers. If we want to suppress - // it, we have to include this file and add a - // single line somewhere to the program (see - // the main() function below for that): -#include - - // The final step, as in previous - // programs, is to import all the - // deal.II class and function names - // into the global namespace: using namespace dealii; - // @sect3{The Step4 class template} - - // This is again the same - // Step4 class as in the - // previous example. The only - // difference is that we have now - // declared it as a class with a - // template parameter, and the - // template parameter is of course - // the spatial dimension in which we - // would like to solve the Laplace - // equation. Of course, several of - // the member variables depend on - // this dimension as well, in - // particular the Triangulation - // class, which has to represent - // quadrilaterals or hexahedra, - // respectively. Apart from this, - // everything is as before. + // @sect3{The Step41 class template} + + // This class supply all function and variables + // to an obstacle problem. The projection_active_set + // function and the ConstaintMatrix are important + // for the handling of the active set as we see + // later. + template -class Step4 +class Step41 { public: - Step4 (); + Step41 (); void run (); private: @@ -119,55 +91,6 @@ class Step4 // @sect3{Right hand side and boundary values} - // In the following, we declare two more - // classes denoting the right hand side and - // the non-homogeneous Dirichlet boundary - // values. Both are functions of a - // dim-dimensional space variable, so we - // declare them as templates as well. - // - // Each of these classes is derived from a - // common, abstract base class Function, - // which declares the common interface which - // all functions have to follow. In - // particular, concrete classes have to - // overload the value function, - // which takes a point in dim-dimensional - // space as parameters and shall return the - // value at that point as a - // double variable. - // - // The value function takes a - // second argument, which we have here named - // component: This is only meant - // for vector valued functions, where you may - // want to access a certain component of the - // vector at the point - // p. However, our functions are - // scalar, so we need not worry about this - // parameter and we will not use it in the - // implementation of the functions. Inside - // the library's header files, the Function - // base class's declaration of the - // value function has a default - // value of zero for the component, so we - // will access the value - // function of the right hand side with only - // one parameter, namely the point where we - // want to evaluate the function. A value for - // the component can then simply be omitted - // for scalar functions. - // - // Note that the C++ language forces - // us to declare and define a - // constructor to the following - // classes even though they are - // empty. This is due to the fact - // that the base class has no default - // constructor (i.e. one without - // arguments), even though it has a - // constructor which has default - // values for all arguments. template class RightHandSide : public Function { @@ -178,8 +101,6 @@ class RightHandSide : public Function const unsigned int component = 0) const; }; - - template class BoundaryValues : public Function { @@ -203,38 +124,8 @@ class Obstacle : public Function // For this example, we choose as right hand - // side function to function $4(x^4+y^4)$ in - // 2D, or $4(x^4+y^4+z^4)$ in 3D. We could - // write this distinction using an - // if-statement on the space dimension, but - // here is a simple way that also allows us - // to use the same function in 1D (or in 4D, - // if you should desire to do so), by using a - // short loop. Fortunately, the compiler - // knows the size of the loop at compile time - // (remember that at the time when you define - // the template, the compiler doesn't know - // the value of dim, but when it later - // encounters a statement or declaration - // RightHandSide@<2@>, it will take the - // template, replace all occurrences of dim - // by 2 and compile the resulting function); - // in other words, at the time of compiling - // this function, the number of times the - // body will be executed is known, and the - // compiler can optimize away the overhead - // needed for the loop and the result will be - // as fast as if we had used the formulas - // above right away. - // - // The last thing to note is that a - // Point@ denotes a point in - // dim-dimensionsal space, and its individual - // components (i.e. $x$, $y$, - // ... coordinates) can be accessed using the - // () operator (in fact, the [] operator will - // work just as well) with indices starting - // at zero as usual in C and C++. + // side function a constant force density + // like the gravitation attraction. template double RightHandSide::value (const Point &p, const unsigned int /*component*/) const @@ -245,13 +136,7 @@ double RightHandSide::value (const Point &p, } - // As boundary values, we choose x*x+y*y in - // 2D, and x*x+y*y+z*z in 3D. This happens to - // be equal to the square of the vector from - // the origin to the point at which we would - // like to evaluate the function, - // irrespective of the dimension. So that is - // what we return: + // As boundary values, we choose the zero. template double BoundaryValues::value (const Point &p, const unsigned int /*component*/) const @@ -261,6 +146,11 @@ double BoundaryValues::value (const Point &p, return return_value; } + + // The obstacle function describes a cascaded + // barrier. So if the gravitation attraction + // pulls the membrane down it blows over the + // steps. template double Obstacle::value (const Point &p, const unsigned int /*component*/) const @@ -281,86 +171,25 @@ double Obstacle::value (const Point &p, - // @sect3{Implementation of the Step4 class} - - // Next for the implementation of the class - // template that makes use of the functions - // above. As before, we will write everything - // as templates that have a formal parameter - // dim that we assume unknown at - // the time we define the template - // functions. Only later, the compiler will - // find a declaration of - // Step4@<2@> (in the - // main function, actually) and - // compile the entire class with - // dim replaced by 2, a process - // referred to as `instantiation of a - // template'. When doing so, it will also - // replace instances of - // RightHandSide@ by - // RightHandSide@<2@> and - // instantiate the latter class from the - // class template. - // - // In fact, the compiler will also find a - // declaration - // Step4@<3@> in - // main(). This will cause it to - // again go back to the general - // Step4@ - // template, replace all occurrences of - // dim, this time by 3, and - // compile the class a second time. Note that - // the two instantiations - // Step4@<2@> and - // Step4@<3@> are - // completely independent classes; their only - // common feature is that they are both - // instantiated from the same general - // template, but they are not convertible - // into each other, for example, and share no - // code (both instantiations are compiled - // completely independently). - - - // @sect4{Step4::Step4} - - // After this introduction, here is the - // constructor of the Step4 - // class. It specifies the desired polynomial - // degree of the finite elements and - // associates the DoFHandler to the - // triangulation just as in the previous - // example program, step-3: + // @sect3{Implementation of the Step41 class} + + + // @sect4{Step41::Step41} + template -Step4::Step4 () +Step41::Step41 () : fe (1), dof_handler (triangulation) {} - // @sect4{Step4::make_grid} - - // Grid creation is something inherently - // dimension dependent. However, as long as - // the domains are sufficiently similar in 2D - // or 3D, the library can abstract for - // you. In our case, we would like to again - // solve on the square $[-1,1]\times [-1,1]$ - // in 2D, or on the cube $[-1,1] \times - // [-1,1] \times [-1,1]$ in 3D; both can be - // termed GridGenerator::hyper_cube(), so we may - // use the same function in whatever - // dimension we are. Of course, the functions - // that create a hypercube in two and three - // dimensions are very much different, but - // that is something you need not care - // about. Let the library handle the - // difficult things. + // @sect4{Step41::make_grid} + + // We solve our obstacle problem on the square + // $[-1,1]\times [-1,1]$ in 2D. template -void Step4::make_grid () +void Step41::make_grid () { GridGenerator::hyper_cube (triangulation, -1, 1); n_refinements = 6; @@ -374,19 +203,10 @@ void Step4::make_grid () << std::endl; } - // @sect4{Step4::setup_system} - - // This function looks - // exactly like in the previous example, - // although it performs actions that in their - // details are quite different if - // dim happens to be 3. The only - // significant difference from a user's - // perspective is the number of cells - // resulting, which is much higher in three - // than in two space dimensions! + // @sect4{Step41::setup_system} + template -void Step4::setup_system () +void Step41::setup_system () { dof_handler.distribute_dofs (fe); @@ -411,78 +231,27 @@ void Step4::setup_system () } - // @sect4{Step4::assemble_system} + // @sect4{Step41::assemble_system} - // Unlike in the previous example, we - // would now like to use a - // non-constant right hand side - // function and non-zero boundary - // values. Both are tasks that are - // readily achieved with a only a few - // new lines of code in the - // assemblage of the matrix and right - // hand side. - // - // More interesting, though, is the - // way we assemble matrix and right - // hand side vector dimension - // independently: there is simply no - // difference to the - // two-dimensional case. Since the - // important objects used in this - // function (quadrature formula, - // FEValues) depend on the dimension - // by way of a template parameter as - // well, they can take care of - // setting up properly everything for - // the dimension for which this - // function is compiled. By declaring - // all classes which might depend on - // the dimension using a template - // parameter, the library can make - // nearly all work for you and you - // don't have to care about most - // things. + + // At once with assembling the system matrix and + // right-hand-side we apply the constraints + // to our system. The constraint consists not + // only of the zero Dirichlet boundary values, + // in addition they contain the obstacle values. + // The projection_active_set function are used + // to fill the ConstraintMatrix. template -void Step4::assemble_system () +void Step41::assemble_system () { QGauss quadrature_formula(2); - // We wanted to have a non-constant right - // hand side, so we use an object of the - // class declared above to generate the - // necessary data. Since this right hand - // side object is only used locally in the - // present function, we declare it here as - // a local variable: const RightHandSide right_hand_side; - // Compared to the previous example, in - // order to evaluate the non-constant right - // hand side function we now also need the - // quadrature points on the cell we are - // presently on (previously, we only - // required values and gradients of the - // shape function from the - // FEValues object, as well as - // the quadrature weights, - // FEValues::JxW() ). We can tell the - // FEValues object to do for - // us by also giving it the - // #update_quadrature_points - // flag: FEValues fe_values (fe, quadrature_formula, update_values | update_gradients | update_quadrature_points | update_JxW_values); - // We then again define a few - // abbreviations. The values of these - // variables of course depend on the - // dimension which we are presently - // using. However, the FE and Quadrature - // classes do all the necessary work for - // you and you don't have to care about the - // dimension dependent parts: const unsigned int dofs_per_cell = fe.dofs_per_cell; const unsigned int n_q_points = quadrature_formula.size(); @@ -491,19 +260,6 @@ void Step4::assemble_system () std::vector local_dof_indices (dofs_per_cell); - // Next, we again have to loop over all - // cells and assemble local contributions. - // Note, that a cell is a quadrilateral in - // two space dimensions, but a hexahedron - // in 3D. In fact, the - // active_cell_iterator data - // type is something different, depending - // on the dimension we are in, but to the - // outside world they look alike and you - // will probably never see a difference - // although the classes that this typedef - // stands for are in fact completely - // unrelated: typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); @@ -514,28 +270,6 @@ void Step4::assemble_system () cell_matrix = 0; cell_rhs = 0; - // Now we have to assemble the - // local matrix and right hand - // side. This is done exactly - // like in the previous - // example, but now we revert - // the order of the loops - // (which we can safely do - // since they are independent - // of each other) and merge the - // loops for the local matrix - // and the local vector as far - // as possible to make - // things a bit faster. - // - // Assembling the right hand side - // presents the only significant - // difference to how we did things in - // step-3: Instead of using a constant - // right hand side with value 1, we use - // the object representing the right - // hand side and evaluate it at the - // quadrature points: for (unsigned int q_point=0; q_point::assemble_system () right_hand_side.value (fe_values.quadrature_point (q_point)) * fe_values.JxW (q_point)); } - // As a final remark to these loops: - // when we assemble the local - // contributions into - // cell_matrix(i,j), we - // have to multiply the gradients of - // shape functions $i$ and $j$ at point - // q_point and multiply it with the - // scalar weights JxW. This is what - // actually happens: - // fe_values.shape_grad(i,q_point) - // returns a dim - // dimensional vector, represented by a - // Tensor@<1,dim@> object, - // and the operator* that multiplies it - // with the result of - // fe_values.shape_grad(j,q_point) - // makes sure that the dim - // components of the two vectors are - // properly contracted, and the result - // is a scalar floating point number - // that then is multiplied with the - // weights. Internally, this operator* - // makes sure that this happens - // correctly for all dim - // components of the vectors, whether - // dim be 2, 3, or any - // other space dimension; from a user's - // perspective, this is not something - // worth bothering with, however, - // making things a lot simpler if one - // wants to write code dimension - // independently. - - // With the local systems assembled, - // the transfer into the global matrix - // and right hand side is done exactly - // as before, but here we have again - // merged some loops for efficiency: - cell->get_dof_indices (local_dof_indices); + + cell->get_dof_indices (local_dof_indices); + + // This function apply the constraints + // to the system matrix and system rhs. + // The true parameter is set to make sure + // that the system rhs contains correct + // values in the rows with inhomogeneity + // constraints. constraints.distribute_local_to_global (cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs, true); } } - // @sect4{Step4::projection_active_set} + // @sect4{Step41::projection_active_set} - // Projection and updating of the active set - // for the dofs which penetrates the obstacle. + // Updating of the active set which means to + // set a inhomogeneity constraint in the + // ConstraintMatrix. At the same time we set + // the solution to the correct value - the obstacle value. + // To control the active set we use the vector + // active_set which contains a zero in a component + // that is not in the active set and elsewise a + // one. With the output file you can visualize it. template -void Step4::projection_active_set () +void Step41::projection_active_set () { const Obstacle obstacle; std::vector vertex_touched (triangulation.n_vertices(), @@ -624,10 +334,12 @@ void Step4::projection_active_set () double obstacle_value = obstacle.value (point); double solution_index_x = solution (index_x); - // to decide which dof belongs to the - // active-set. for that we scale the + // To decide which dof belongs to the + // active-set. For that we scale the // residual-vector with the cell-size and // the diag-entry of the mass-matrix. + + // TODO: I have to check the condition if ((resid_vector (index_x)*std::pow (2, 2*n_refinements)*diag_mass_matrix_vector (index_x) >= solution_index_x - obstacle_value)) { constraints.add_line (index_x); @@ -644,7 +356,7 @@ void Step4::projection_active_set () } std::cout<< "Number of Contact-Constaints: " << counter_contact_constraints <::projection_active_set () constraints.close (); } - // @sect4{Step4::solve} + // @sect4{Step41::solve} - // Solving the linear system of - // equations is something that looks - // almost identical in most - // programs. In particular, it is - // dimension independent, so this - // function is copied verbatim from the - // previous example. template -void Step4::solve () +void Step41::solve () { ReductionControl reduction_control (100, 1e-12, 1e-2); SolverCG solver (reduction_control); @@ -679,35 +384,13 @@ void Step4::solve () << std::endl; } - // @sect4{Step4::output_results} + // @sect4{Step41::output_results} - // This function also does what the - // respective one did in step-3. No changes - // here for dimension independence either. - // - // The only difference to the previous - // example is that we want to write output in - // VTK format, rather than for gnuplot. VTK - // format is currently the most widely used - // one and is supported by a number of - // visualization programs such as Visit and - // Paraview (for ways to obtain these - // programs see the ReadMe file of - // deal.II). To write data in this format, we - // simply replace the - // data_out.write_gnuplot call - // by data_out.write_vtk. - // - // Since the program will run both 2d and 3d - // versions of the laplace solver, we use the - // dimension in the filename to generate - // distinct filenames for each run (in a - // better program, one would check whether - // dim can have other values - // than 2 or 3, but we neglect this here for - // the sake of brevity). + // We use the vtk-format for the output. + // The file contains the displacement, + // the residual and active set vectors. template -void Step4::output_results (const std::string& title) const +void Step41::output_results (const std::string& title) const { DataOut data_out; @@ -718,23 +401,22 @@ void Step4::output_results (const std::string& title) const data_out.build_patches (); - std::ofstream output_vtk (dim == 2 ? - (title + ".vtk").c_str () : - (title + ".vtk").c_str ()); + std::ofstream output_vtk ((title + ".vtk").c_str ()); data_out.write_vtk (output_vtk); } - // @sect4{Step4::run} + // @sect4{Step41::run} // This is the function which has the - // top-level control over - // everything. Apart from one line of - // additional output, it is the same - // as for the previous example. + // top-level control over everything. + // Here the active set method is implemented. + + // TODO: I have to compare it with the algorithm + // in the Wohlmuth-paper template -void Step4::run () +void Step41::run () { std::cout << "Solving problem in " << dim << " space dimensions." << std::endl; @@ -812,7 +494,7 @@ void Step4::run () // looks mostly like in step-3, but if you // look at the code below, note how we first // create a variable of type - // Step4@<2@> (forcing + // Step41@<2@> (forcing // the compiler to compile the class template // with dim replaced by // 2) and run a 2d simulation, @@ -887,7 +569,7 @@ int main (int argc, char *argv[]) Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv); - Step4<2> laplace_problem_2d; + Step41<2> laplace_problem_2d; laplace_problem_2d.run (); return 0; -- 2.39.5