From 61d22ffd9b7148528451ec3648731dbea5c06418 Mon Sep 17 00:00:00 2001
From: guido
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
+
+
+ -
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