From 5bf5a9ec2288bb287845d8fa9fe99fd9c77ac05a Mon Sep 17 00:00:00 2001 From: Guido Kanschat Date: Thu, 3 Jul 2003 10:31:47 +0000 Subject: [PATCH] dof_renumbering.h moved git-svn-id: https://svn.dealii.org/trunk@7833 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/Makefile | 12 +- .../deal.II/source/dofs/dof_renumbering.cc | 2 +- deal.II/doc/news/c-4-0.html | 16 + deal.II/examples/step-16/step-16.cc | 349 +++++++----------- deal.II/examples/step-2/step-2.cc | 2 +- deal.II/examples/step-7/step-7.cc | 2 +- deal.II/lac/include/lac/solver_gmres.h | 189 +++++++++- tests/bits/anna_1.cc | 2 +- tests/bits/anna_2.cc | 2 +- tests/bits/anna_3.cc | 2 +- tests/deal.II/dof_renumbering.cc | 4 +- tests/deal.II/dof_test.cc | 2 +- tests/deal.II/intergrid_constraints.cc | 2 +- tests/deal.II/wave-test-3.cc | 2 +- 14 files changed, 359 insertions(+), 229 deletions(-) diff --git a/deal.II/Makefile b/deal.II/Makefile index 7be19f8c6b..4feff8d2f2 100644 --- a/deal.II/Makefile +++ b/deal.II/Makefile @@ -117,11 +117,13 @@ TODO: > $@ -TAGS: $(shell echo $D/base/source/*.cc $D/base/include/*/*.h \ - $D/lac/source/*.cc $D/lac/include/*/*.h \ - $D/deal.II/source/*/*.cc $D/deal.II/include/*/*.h \ - $D/tests/*/*.cc $D/tests/*/*.h ) - @cd $D/common ; etags $^ +TAGS: + @cd $D/common ; etags $D/base/source/*.cc $D/base/include/*/*.h \ + $D/lac/source/*.cc $D/lac/include/*/*.h \ + $D/deal.II/source/*/*.cc $D/deal.II/include/*/*.h \ + $D/tests/*/*.cc $D/tests/*/*.h \ + $D/examples/*/*.cc + clean: clean-contrib clean-base clean-lac clean-dealII clean-doc clean-examples clean-lib clean-tests diff --git a/deal.II/deal.II/source/dofs/dof_renumbering.cc b/deal.II/deal.II/source/dofs/dof_renumbering.cc index e392281427..3c72818669 100644 --- a/deal.II/deal.II/source/dofs/dof_renumbering.cc +++ b/deal.II/deal.II/source/dofs/dof_renumbering.cc @@ -26,7 +26,7 @@ #include #include #include -#include +#include #ifdef ENABLE_MULTIGRID #include diff --git a/deal.II/doc/news/c-4-0.html b/deal.II/doc/news/c-4-0.html index 9743a182d7..35d6d6da6a 100644 --- a/deal.II/doc/news/c-4-0.html +++ b/deal.II/doc/news/c-4-0.html @@ -28,6 +28,22 @@ contributor's names are abbreviated by WB (Wolfgang Bangerth), GK

+ +

Incompatibilities

+ +
    +
  1. Moved and changed: The header file + include/numerics/dof_renumbering.h has been moved to the + directory include/dofs, where it logically + belongs. Furthermore, the sorting parameter of the function DoFRenumbering::component_wise has changed its meaning. See + the reference documentation for details. +
    + (GK 2003/07/03) +

+ +

General

diff --git a/deal.II/examples/step-16/step-16.cc b/deal.II/examples/step-16/step-16.cc index 0fc88de0b4..a4c753b7bc 100644 --- a/deal.II/examples/step-16/step-16.cc +++ b/deal.II/examples/step-16/step-16.cc @@ -80,6 +80,10 @@ class LaplaceProblem private: void setup_system (); void assemble_system (); + // We add this function for + // assembling the multilevel + // matrices. + void assemble_multigrid (); void solve (); void output_results (const unsigned int cycle) const; @@ -90,7 +94,12 @@ class LaplaceProblem SparsityPattern sparsity_pattern; SparseMatrix system_matrix; + // Here are the new objects for + // handling level matrices. MGLevelObject mg_sparsity; + // We use number type float to + // save memory. It's only a + // preconditioner! MGLevelObject > mg_matrices; Vector solution; @@ -140,41 +149,33 @@ void LaplaceProblem::setup_system () const unsigned int nlevels = triangulation.n_levels(); mg_sparsity.resize(0, nlevels-1); mg_matrices.resize(0, nlevels-1); - + + // Now, we have to build a matrix + // on each level. Technically, we + // could use the matrix initialized + // above on the finest + // level. Beware that this is not + // true anymore with local + // refinement! for (unsigned int level=0;level void LaplaceProblem::assemble_system () { @@ -270,6 +271,95 @@ void LaplaceProblem::assemble_system () } + // Here is another assemble + // function. The integration core is + // the same as above. Only the loop + // goes over all existing cells now + // and the results must be entered + // into the right matrix. + + // Since we only do multi-level + // preconditioning, no right-hand + // side is assembled here. +template +void LaplaceProblem::assemble_multigrid () +{ + QGauss2 quadrature_formula; + + FEValues fe_values (fe, quadrature_formula, + UpdateFlags(update_gradients | + update_q_points | + update_JxW_values)); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.n_quadrature_points; + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + // + typename DoFHandler::cell_iterator cell = mg_dof_handler.begin(), + endc = mg_dof_handler.end(); + for (; cell!=endc; ++cell) + { + // Remember the level of the + // current cell. + const unsigned int level = cell->level() + cell_matrix.clear (); + + // Compute the values specified + // by update flags above. + fe_values.reinit (cell); + + // This is exactly the + // integration loop of the cell + // matrix above. + for (unsigned int q_point=0; q_pointget_mg_dof_indices (local_dof_indices); + for (unsigned int i=0; i boundary_values; + VectorTools::interpolate_boundary_values (mg_dof_handler, + 0, + ZeroFunction(), + boundary_values); + MatrixTools::apply_boundary_values (boundary_values, + system_matrix, + solution, + system_rhs); +} + + // The solution process again looks // mostly like in the previous @@ -332,184 +422,47 @@ void LaplaceProblem::solve () - // Writing output to a file is mostly - // the same as for the previous - // example, but here we will show how - // to modify some output options and - // how to construct a different - // filename for each refinement - // cycle. + // Here is the data output, which is + // a simplified version of step 5. We + // do a standard gnuplot output for + // each grid produced in the + // refinement process. template void LaplaceProblem::output_results (const unsigned int cycle) const { + // Construct and initialize a DataOut object DataOut data_out; data_out.attach_dof_handler (mg_dof_handler); data_out.add_data_vector (solution, "solution"); - data_out.build_patches (); - // For this example, we would like - // to write the output directly to - // a file in Encapsulated - // Postscript (EPS) format. The - // library supports this, but - // things may be a bit more - // difficult sometimes, since EPS - // is a printing format, unlike - // most other supported formats - // which serve as input for - // graphical tools. Therefore, you - // can't scale or rotate the image - // after it has been written to - // disk, and you have to decide - // about the viewpoint or the - // scaling in advance. - // - // The defaults in the library are - // usually quite reasonable, and - // regarding viewpoint and scaling - // they coincide with the defaults - // of Gnuplot. However, since this - // is a tutorial, we will - // demonstrate how to change - // them. For this, we first have to - // generate an object describing - // the flags for EPS output: - DataOutBase::EpsFlags eps_flags; - // They are initialized with the - // default values, so we only have - // to change those that we don't - // like. For example, we would like - // to scale the z-axis differently - // (stretch each data point in - // z-direction by a factor of four): - eps_flags.z_scaling = 4; - // Then we would also like to alter - // the viewpoint from which we look - // at the solution surface. The - // default is at an angle of 60 - // degrees down from the vertical - // axis, and 30 degrees rotated - // against it in mathematical - // positive sense. We raise our - // viewpoint a bit and look more - // along the y-axis: - eps_flags.azimut_angle = 40; - eps_flags.turn_angle = 10; - // That shall suffice. There are - // more flags, for example whether - // to draw the mesh lines, which - // data vectors to use for - // colorization of the interior of - // the cells, and so on. You may - // want to take a look at the - // documentation of the EpsFlags - // structure to get an overview of - // what is possible. - // - // The only thing still to be done, - // is to tell the output object to - // use these flags: - data_out.set_flags (eps_flags); - // The above way to modify flags - // requires recompilation each time - // we would like to use different - // flags. This is inconvenient, and - // we will see more advanced ways - // in following examples where the - // output flags are determined at - // run time using an input file. - - // Finally, we need the filename to - // which the results are to be - // written. We would like to have - // it of the form - // ``solution-N.eps'', where N is - // the number of the refinement - // cycle. Thus, we have to convert - // an integer to a part of a - // string; this can be done using - // the ``sprintf'' function, but in - // C++ there is a more elegant way: - // write everything into a special - // stream (just like writing into a - // file or to the screen) and - // retrieve what you wrote as a - // string. This applies the usual - // conversions from integer to - // strings, and one could as well - // give stream modifiers such as - // ``setw'', ``setprecision'', and - // so on. - // - // In C++, you can do this by using - // the so-called stringstream - // classes. As already discussed at - // the point of inclusion of the - // respective header file above, - // there is some historical - // confusion we have to work around - // here, since the class we'd like - // to use used to be called - // ``ostrstream'', but now is named - // ``ostringstream''. In the same - // way as done above in deciding - // which file to include, we here - // decide which class name to use: + // The following block generates + // the file name and opens the + // file. This looks awkward because + // of a change in the handling of + // string streams (See step 5 for explanation). + #ifdef HAVE_STD_STRINGSTREAM std::ostringstream filename; #else std::ostrstream filename; #endif - // Fortunately, the interface of - // the two classes which we might - // now be using, depending on which - // one is available, is close - // enough that we need to take care - // about the differences only once - // below, so we can use them in a - // rather straightforward way, even - // if they are not identical. - - // In order to now actually - // generate a filename, we fill the - // stringstream variable with the - // base of the filename, then the - // number part, and finally the - // suffix indicating the file type: filename << "solution-" << cycle - << ".eps"; - - // For the old string stream - // classes, we have to append the - // final '\0' that appears at the - // end of ``char *'' - // variables. This is done by the - // following construct: + << ".gnuplot"; + #ifndef HAVE_STD_STRINGSTREAM filename << std::ends; #endif - // We can get whatever we wrote to - // the stream using the ``str()'' - // function. If the new - // stringstream classes are used, - // then the result is a string - // which we have to convert to a - // char* using the ``c_str()'' - // function, otherwise the result - // is a char* right away. Use that - // as filename for the output - // stream: + #ifdef HAVE_STD_STRINGSTREAM std::ofstream output (filename.str().c_str()); #else std::ofstream output (filename.str()); #endif - // And then write the data to the - // file. - data_out.write_eps (output); + + data_out.write_gnuplot (output); } @@ -521,35 +474,11 @@ void LaplaceProblem::run () { std::cout << "Cycle " << cycle << ':' << std::endl; - // If this is the first round, - // then we have no grid yet, - // and we will create it - // here. In previous examples, - // we have already used some of - // the functions from the - // GridGenerator class. Here we - // would like to read a grid - // from a file where the cells - // are stored and which may - // originate from someone else, - // or may be the product of a - // mesh generator tool. - // - // In order to read a grid from - // a file, we generate an - // object of data type GridIn - // and associate the - // triangulation to it (i.e. we - // tell it to fill our - // triangulation object when we - // ask it to read the - // file). Then we open the - // respective file and - // initialize the triangulation - // with the data in the file: if (cycle == 0) { - // Generate grid here! + // Generate a simple hypercube grid. + GridGenerator::hyper_cube(tr); + } // If this is not the first // cycle, then simply refine diff --git a/deal.II/examples/step-2/step-2.cc b/deal.II/examples/step-2/step-2.cc index 6b29f90afa..bfef563132 100644 --- a/deal.II/examples/step-2/step-2.cc +++ b/deal.II/examples/step-2/step-2.cc @@ -56,7 +56,7 @@ // We will want to use a special // algorithm to renumber degrees of // freedom. It is declared here: -#include +#include // This is needed for C++ output: #include diff --git a/deal.II/examples/step-7/step-7.cc b/deal.II/examples/step-7/step-7.cc index d0462ca73d..34174a25fb 100644 --- a/deal.II/examples/step-7/step-7.cc +++ b/deal.II/examples/step-7/step-7.cc @@ -47,7 +47,7 @@ // Cuthill-McKee algorithm. The // necessary functions are declared // in the following file: -#include +#include // Then we will show a little trick // how we can make sure that objects // are not deleted while they are diff --git a/deal.II/lac/include/lac/solver_gmres.h b/deal.II/lac/include/lac/solver_gmres.h index 88a7b882db..3f189b14d1 100644 --- a/deal.II/lac/include/lac/solver_gmres.h +++ b/deal.II/lac/include/lac/solver_gmres.h @@ -14,7 +14,6 @@ #define __deal2__solver_gmres_h -/*---------------------------- solver_pgmres.h ---------------------------*/ #include #include @@ -289,6 +288,85 @@ class SolverGMRES : public Solver }; +/** + * Flexible GMRES. + * + * This version of the GMRES method allows for the use of a different + * preconditioner in each iteration step. Therefore, it is also more + * robust with respect to inaccurate evaluation of the + * preconditioner. An important application is also the use of a + * Krylov space method inside the preconditioner. + * + * FGMRES needs two vectors in each iteration steps yielding a total + * of @p{2*AdditionalData::max_basis_size+1} auxiliary vectors. + */ +template > +class SolverFGMRES : public Solver +{ + public: + /** + * Standardized data struct to + * pipe additional data to the + * solver. + */ + struct AdditionalData + { + /** + * Constructor. By default, + * set the number of + * temporary vectors to 30, + * preconditioning from left + * and the residual of the + * stopping criterion to the + * default residual + * (cf. class documentation). + */ + AdditionalData(const unsigned int max_basis_size = 30, + bool use_default_residual = true) + : + max_basis_size(max_basis_size) + {}; + + /** + * Maximum number of + * tmp vectors. + */ + unsigned int max_basis_size; + }; + + /** + * Constructor. + */ + SolverFGMRES (SolverControl &cn, + VectorMemory &mem, + const AdditionalData &data=AdditionalData()); + + /** + * Solve the linear system $Ax=b$ + * for x. + */ + template + void + solve (const MATRIX &A, + VECTOR &x, + const VECTOR &b, + const PRECONDITIONER &precondition); + + private: + /** + * Additional flags. + */ + AdditionalData additional_data; + /** + * Projected system matrix + */ + FullMatrix H; + /** + * Auxiliary matrix for inverting @p{H} + */ + FullMatrix H1; +}; + /* --------------------- Inline and template functions ------------------- */ @@ -708,7 +786,112 @@ SolverGMRES::criterion () } -/*---------------------------- solver_pgmres.h ---------------------------*/ +//----------------------------------------------------------------------// + +template +SolverFGMRES::SolverFGMRES (SolverControl &cn, + VectorMemory &mem, + const AdditionalData &data) + : + Solver (cn, mem), + additional_data(data) +{} + + +template +template +void +SolverFGMRES::solve ( + const MATRIX& A, + VECTOR& x, + const VECTOR& b, + const PRECONDITIONER& precondition) +{ + deallog.push("FGMRES"); + + SolverControl::State iteration_state = SolverControl::iterate; + + const unsigned int basis_size = additional_data.max_basis_size; + + // Generate an object where basis + // vectors are stored. + typename SolverGMRES::TmpVectors v (basis_size, this->memory); + typename SolverGMRES::TmpVectors z (basis_size, this->memory); + + // number of the present iteration; this + // number is not reset to zero upon a + // restart + unsigned int accumulated_iterations = 0; + + // matrix used for the orthogonalization + // process later + H.reinit(basis_size+1, basis_size); + + // Vectors for projected system + Vector projected_rhs; + Vector y; + + // Iteration starts here + + do + { + VECTOR* aux = this->memory.alloc(); + aux->reinit(x); + A.vmult(*aux, x); + aux->sadd(-1., 1., b); + + double beta = aux->l2_norm(); + if (this->control().check(accumulated_iterations,beta) + == SolverControl::success) + break; + + H.reinit(basis_size+1, basis_size); + double a = beta; + + for (unsigned int j=0;jadd(-H(i,j), v[i]); + } + H(j+1,j) = a = aux->l2_norm(); + + // Compute projected solution + + if (j>0) + { + H1.reinit(j+1,j); + projected_rhs.reinit(j+1); + y.reinit(j); + projected_rhs(0) = beta; + H1.fill(H); + + double res = H1.least_squares(y, projected_rhs); + iteration_state = this->control().check(++accumulated_iterations, res); + if (iteration_state != SolverControl::iterate) + break; + } + } + // Update solution vector + for (unsigned int j=0;jmemory.free(aux); + } while (iteration_state == SolverControl::iterate); + + deallog.pop(); + // in case of failure: throw + // exception + if (this->control().last_check() != SolverControl::success) + throw SolverControl::NoConvergence (this->control().last_step(), + this->control().last_value()); +} #endif -/*---------------------------- solver_pgmres.h ---------------------------*/ diff --git a/tests/bits/anna_1.cc b/tests/bits/anna_1.cc index 6d64cda088..98e2f417c7 100644 --- a/tests/bits/anna_1.cc +++ b/tests/bits/anna_1.cc @@ -33,7 +33,7 @@ #include #include #include -#include +#include #include #include diff --git a/tests/bits/anna_2.cc b/tests/bits/anna_2.cc index 96f55918e3..0cd178873c 100644 --- a/tests/bits/anna_2.cc +++ b/tests/bits/anna_2.cc @@ -37,7 +37,7 @@ #include #include #include -#include +#include #include #include diff --git a/tests/bits/anna_3.cc b/tests/bits/anna_3.cc index 59b688b0bf..c9baa8864f 100644 --- a/tests/bits/anna_3.cc +++ b/tests/bits/anna_3.cc @@ -32,7 +32,7 @@ #include #include #include -#include +#include #include #include diff --git a/tests/deal.II/dof_renumbering.cc b/tests/deal.II/dof_renumbering.cc index 93568b02df..72c55b5dcf 100644 --- a/tests/deal.II/dof_renumbering.cc +++ b/tests/deal.II/dof_renumbering.cc @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2000, 2001 by the deal.II authors +// Copyright (C) 2000, 2001, 2003 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -24,12 +24,12 @@ #include #include #include +#include #include #include #include #include #include -#include #include diff --git a/tests/deal.II/dof_test.cc b/tests/deal.II/dof_test.cc index c28798dc7f..01f4c4202e 100644 --- a/tests/deal.II/dof_test.cc +++ b/tests/deal.II/dof_test.cc @@ -22,7 +22,7 @@ #include #include #include -#include +#include #include #include diff --git a/tests/deal.II/intergrid_constraints.cc b/tests/deal.II/intergrid_constraints.cc index 65cc7eb669..21841617e8 100644 --- a/tests/deal.II/intergrid_constraints.cc +++ b/tests/deal.II/intergrid_constraints.cc @@ -9,7 +9,7 @@ #include #include #include -#include +#include #include #include #include diff --git a/tests/deal.II/wave-test-3.cc b/tests/deal.II/wave-test-3.cc index 15eed7412a..a9428921d9 100644 --- a/tests/deal.II/wave-test-3.cc +++ b/tests/deal.II/wave-test-3.cc @@ -4468,7 +4468,7 @@ void SweepInfo::write_summary (const std::list*> &eval_list, #include #include #include -#include +#include #include -- 2.39.5