From: wolf
+ Step 5 +
+Back to the tutorial index diff --git a/deal.II/doc/tutorial/chapter-2.step-by-step/step-4.results b/deal.II/doc/tutorial/chapter-2.step-by-step/step-4.results index 3201c10650..7514fac593 100644 --- a/deal.II/doc/tutorial/chapter-2.step-by-step/step-4.results +++ b/deal.II/doc/tutorial/chapter-2.step-by-step/step-4.results @@ -1,7 +1,10 @@
Solving problem in 2 space dimensions.
diff --git a/deal.II/doc/tutorial/chapter-2.step-by-step/step-5.intro b/deal.II/doc/tutorial/chapter-2.step-by-step/step-5.intro
new file mode 100644
index 0000000000..1518606e43
--- /dev/null
+++ b/deal.II/doc/tutorial/chapter-2.step-by-step/step-5.intro
@@ -0,0 +1,3 @@
+
+Introduction
+
diff --git a/deal.II/doc/tutorial/chapter-2.step-by-step/step-5.results b/deal.II/doc/tutorial/chapter-2.step-by-step/step-5.results
new file mode 100644
index 0000000000..747bdbbb3d
--- /dev/null
+++ b/deal.II/doc/tutorial/chapter-2.step-by-step/step-5.results
@@ -0,0 +1,3 @@
+
+Results
+
diff --git a/deal.II/doc/tutorial/chapter-2.step-by-step/toc.html b/deal.II/doc/tutorial/chapter-2.step-by-step/toc.html
index 32b036a221..0bec85471b 100644
--- a/deal.II/doc/tutorial/chapter-2.step-by-step/toc.html
+++ b/deal.II/doc/tutorial/chapter-2.step-by-step/toc.html
@@ -65,6 +65,13 @@ At present, the following programs exist:
3D, although the program is exactly the same. Non-constant right
hand side function. Non-homogeneous boundary values.
+
+ Step 5
+ What's new: Computations on successively
+ refined grids. Some optimizations. Non-constant coefficient in
+ the elliptic operator. Preconditioning of the CG solver for the
+ linear system of equations.
+
diff --git a/deal.II/examples/step-4/step-4.cc b/deal.II/examples/step-4/step-4.cc
index cf073c5509..5c95b88071 100644
--- a/deal.II/examples/step-4/step-4.cc
+++ b/deal.II/examples/step-4/step-4.cc
@@ -252,7 +252,7 @@ void LaplaceProblem::make_grid_and_dofs ()
cout << " Number of active cells: "
<< triangulation.n_active_cells()
<< endl
- << " Total number of cells: "
+ << " Total number of cells: "
<< triangulation.n_cells()
<< endl;
diff --git a/deal.II/examples/step-5/Makefile b/deal.II/examples/step-5/Makefile
new file mode 100644
index 0000000000..dc2798f839
--- /dev/null
+++ b/deal.II/examples/step-5/Makefile
@@ -0,0 +1,94 @@
+# $Id$
+# Copyright W. Bangerth, University of Heidelberg, 1999
+
+# Template for makefiles for the examples subdirectory. In principle,
+# everything should be done automatically if you set the target file
+# here correctly. We get deduce it from the files in the present
+# directory:
+target = $(basename $(shell echo step-*.cc))
+
+# All dependencies between files should be updated by the included
+# file Makefile.dep if necessary. Object files are compiled into
+# the archives ./Obj.a and ./Obj.g.a. By default, the debug version
+# is used to link. It you don't like that, change the following
+# variable to "off"
+debug-mode = on
+
+
+
+###############################################################################
+# Internals
+
+#deal include base path
+D = ../../../..
+
+include $D/common/Make.global_options
+
+# get lists of files we need
+
+
+# list of libraries needed to link with
+libs = -ldeal_II_2d -llac -lbase
+libs.g = -ldeal_II_2d.g -llac.g -lbase.g
+
+
+# check whether we use debug mode or not
+ifeq ($(debug-mode),on)
+ libraries = $(target).go $(libs.g)
+ flags = $(CXXFLAGS.g)
+else
+ libraries = $(target).go $(libs)
+ flags = $(CXXFLAGS.o)
+endif
+
+
+
+# make rule for the target. $^ is the object file $(target).g?o
+$(target) : $(libraries)
+ @echo ============================ Linking $@
+ @$(CXX) $(flags) -o $@ $^
+
+# rule how to run the program
+run: $(target)
+ @echo ============================ Running $@
+ @./$(target)
+
+
+# rule to make object files
+%.go : %.cc
+ @echo ============================ Compiling with debugging information: $<
+ @$(CXX) $(CXXFLAGS.g) -c $< -o $@
+%.o : %.cc
+ @echo ============================ Compiling with optimization: $<
+ @$(CXX) $(CXXFLAGS) -c $< -o $@
+
+
+clean:
+ -rm -f *.o *.go *~ Makefile.dep $(target) *gmv *gnuplot *gpl *eps
+
+
+
+.PHONY: clean
+
+
+# Rule to generate the dependency file. This file is
+# automagically remade whenever needed, i.e. whenever
+# one of the cc-/h-files changed. Make detects whether
+# to remake this file upon inclusion at the bottom
+# of this file.
+#
+# Since the script prefixes the output names by lib/g?o, we have to
+# strip that again (the script was written for the main libraries and
+# large projects where object files are put into subdirs)
+Makefile.dep: $(target).cc Makefile \
+ $(shell echo $D/base/include/base/*.h \
+ $D/lac/include/lac/*.h \
+ $D/deal.II/include/*/*.h)
+ @echo ============================ Remaking Makefile
+ @perl $D/common/scripts/make_dependencies.pl $(INCLUDE) $(target).cc \
+ | perl -pi -e 's!lib/g?o/!!g;' \
+ > Makefile.dep
+
+
+include Makefile.dep
+
diff --git a/deal.II/examples/step-5/step-5.cc b/deal.II/examples/step-5/step-5.cc
new file mode 100644
index 0000000000..16e2c5c52e
--- /dev/null
+++ b/deal.II/examples/step-5/step-5.cc
@@ -0,0 +1,504 @@
+/* $Id$ */
+
+ // 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
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+#include
+#include
+
+#include
+
+
+template
+class LaplaceProblem
+{
+ public:
+ LaplaceProblem ();
+ void run ();
+
+ private:
+ void make_grid_and_dofs (const unsigned int refinement);
+ void assemble_system ();
+ void solve ();
+ void output_results ();
+ void clear ();
+
+ Triangulation triangulation;
+ FEQ1 fe;
+ DoFHandler dof_handler;
+
+ SparseMatrixStruct sparsity_pattern;
+ SparseMatrix system_matrix;
+
+ Vector solution;
+ Vector system_rhs;
+};
+
+
+template
+class Coefficient : public Function
+{
+ public:
+ virtual double value (const Point &p,
+ const unsigned int component = 0) const;
+ virtual void value_list (const vector > &points,
+ vector &values,
+ const unsigned int component = 0) const;
+};
+
+
+
+template
+double Coefficient::value (const Point &p,
+ const unsigned int) const
+{
+ if (p.square() < 0.5*0.5)
+ return 10;
+ else
+ return 1;
+};
+
+
+
+template
+void Coefficient::value_list (const vector > &points,
+ vector &values,
+ const unsigned int component) const
+{
+ const unsigned int n_points = points.size();
+
+ Assert (values.size() == n_points,
+ ExcVectorHasWrongSize (values.size(), n_points));
+
+ Assert (component == 0,
+ ExcWrongComponent (component, 1));
+
+ for (unsigned int i=0; i
+LaplaceProblem::LaplaceProblem () :
+ dof_handler (triangulation)
+{};
+
+
+
+template
+void LaplaceProblem::make_grid_and_dofs (const unsigned int refinement)
+{
+ GridGenerator::hyper_cube (triangulation, -1, 1);
+ triangulation.refine_global (refinement);
+
+ cout << " Number of active cells: "
+ << triangulation.n_active_cells()
+ << endl
+ << " Total number of cells: "
+ << triangulation.n_cells()
+ << endl;
+
+ dof_handler.distribute_dofs (fe);
+
+ cout << " Number of degrees of freedom: "
+ << dof_handler.n_dofs()
+ << endl;
+
+ sparsity_pattern.reinit (dof_handler.n_dofs(),
+ dof_handler.n_dofs(),
+ dof_handler.max_couplings_between_dofs());
+ DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
+ sparsity_pattern.compress();
+
+ system_matrix.reinit (sparsity_pattern);
+
+ solution.reinit (dof_handler.n_dofs());
+ system_rhs.reinit (dof_handler.n_dofs());
+};
+
+
+
+ // As in the previous examples, this
+ // function is not changed much with
+ // regard to its functionality, but
+ // there are still some optimizations
+ // which we will show. For this, it
+ // is important to note that if
+ // efficient solvers are used (such
+ // as the preconditions CG method),
+ // assembling the matrix and right
+ // hand side can take a comparable
+ // time, and it is worth the effort
+ // to use one or two optimizations at
+ // some places.
+ //
+ // What we will show here is how we
+ // can avoid calls to the
+ // shape_value, shape_grad, and
+ // quadrature_point functions of the
+ // FEValues object, and in particular
+ // optimize away most of the virtual
+ // function calls of the Function
+ // object. The way to do so will be
+ // explained in the following, while
+ // those parts of this function that
+ // are not changed with respect to
+ // the previous example are not
+ // commented on.
+template
+void LaplaceProblem::assemble_system ()
+{
+ // This time, we will again use a
+ // constant right hand side
+ // function, but a variable
+ // coefficient. The following
+ // object will be used for this:
+ const Coefficient coefficient;
+
+ QGauss3 quadrature_formula;
+
+ FEValues fe_values (fe, quadrature_formula,
+ UpdateFlags(update_values |
+ 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);
+ Vector cell_rhs (dofs_per_cell);
+
+ vector local_dof_indices (dofs_per_cell);
+
+ // ...
+ vector coefficient_values (n_q_points);
+
+ DoFHandler::active_cell_iterator cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+ for (; cell!=endc; ++cell)
+ {
+ cell_matrix.clear ();
+ cell_rhs.clear ();
+
+ // As before, we want the
+ // FEValues object to compute
+ // the quantities which we told
+ // him to compute in the
+ // constructor using the update
+ // flags.
+ fe_values.reinit (cell);
+ // Now, these quantities are
+ // stored in arrays in the
+ // FEValues object. Usually,
+ // the internals of how and
+ // where they are stored is not
+ // something that the outside
+ // world should know, but since
+ // this is a time critical
+ // function we decided to
+ // publicize these arrays a
+ // little bit, and provide
+ // facilities to export the
+ // address where this data is
+ // stored.
+ //
+ // For example, the values of
+ // shape function j at
+ // quadrature point q is stored
+ // in a matrix, of which we can
+ // get the address as follows
+ // (note that this is a
+ // reference to the matrix,
+ // symbolized by the ampersand,
+ // and that it must be a
+ // constant reference, since
+ // only read-only access is
+ // granted):
+ const FullMatrix
+ & shape_values = fe_values.get_shape_values();
+ // Instead of writing
+ // fe_values.shape_value(j,q)
+ // we can now write
+ // shape_values(j,q), i.e. the
+ // function call needed
+ // previously for each access
+ // has been otimized away.
+ //
+ // There are alike functions
+ // for almost all data elements
+ // in the FEValues class. The
+ // gradient are accessed as
+ // follows:
+ const vector > >
+ & shape_grads = fe_values.get_shape_grads();
+ // The data type looks a bit
+ // unwieldy, since each entry
+ // in the matrix (j,q) now
+ // needs to be the gradient of
+ // the shape function, which is
+ // a vector.
+ //
+ // Similarly, access to the
+ // place where quadrature
+ // points and the determinants
+ // of the Jacobian matrices
+ // times the weights of the
+ // respective quadrature points
+ // are stored, can be obtained
+ // like this:
+ const vector
+ & JxW_values = fe_values.get_JxW_values();
+ const vector >
+ & q_points = fe_values.get_quadrature_points();
+ // Admittedly, the declarations
+ // above are not easily
+ // readable, but they can save
+ // many function calls in the
+ // inner loops and can thus
+ // make assemblage faster.
+ //
+ // An additional advantage is
+ // that the inner loops are
+ // simpler to read, since the
+ // fe_values object is no more
+ // explicitely needed to access
+ // the different fields (see
+ // below). Unfortunately,
+ // things became a bit
+ // inconsistent, since the
+ // shape values are accessed
+ // via the FullMatrix operator
+ // (), i.e. using parentheses,
+ // while all the other fields
+ // are accessed through vector
+ // operator [], i.e. using
+ // brackets. This is due to
+ // historical reasons and
+ // frequently leads to a bit of
+ // confusion, but since the
+ // places where this happens
+ // are few in well-written
+ // programs, this is not too
+ // big a problem.
+
+ // There is one more thing: in
+ // this example, we want to use
+ // a non-constant
+ // coefficient. In the previous
+ // example, we have called the
+ // ``value'' function of the
+ // right hand side object for
+ // each quadrature
+ // point. Unfortunately, that
+ // is a virtual function, so
+ // calling it is relatively
+ // expensive. Therefore, we use
+ // a function of the Function
+ // class which returns the
+ // values at all quadrature
+ // points at once; that
+ // function is still virtual,
+ // but it needs to be computed
+ // once per cell only, not once
+ // in the inner loop:
+ coefficient.value_list (q_points, coefficient_values);
+ // It should be noted that the
+ // creation of the
+ // coefficient_values object is
+ // done outside the loop over
+ // all cells to avoid memory
+ // allocation each time we
+ // visit a new cell. Contrary
+ // to this, the other variables
+ // above were created inside
+ // the loop, but they were only
+ // references to memory that
+ // has already been allocated
+ // (i.e. they are pointers to
+ // that memory) and therefore,
+ // no new memory needs to be
+ // allocated; in particular, by
+ // declaring the pointers as
+ // close to their use as
+ // possible, we give the
+ // compiler a better choice to
+ // optimize them away
+ // altogether, something which
+ // it definitely can't do with
+ // the coefficient_values
+ // object since it is too
+ // complicated, but mostly
+ // because it's address is
+ // passed to a virtual function
+ // which is not knows at
+ // compile time.
+
+ // Using the various
+ // abbreviations, the loops
+ // then look like this (the
+ // parentheses around the
+ // product of the two gradients
+ // are needed to indicate the
+ // dot product; we have to
+ // overrule associativity of
+ // the operator* here, since
+ // the compiler would otherwise
+ // complain about an undefined
+ // product of double*gradient
+ // since it parses
+ // left-to-right):
+ for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices);
+ for (unsigned int i=0; i boundary_values;
+ VectorTools::interpolate_boundary_values (dof_handler,
+ 0,
+ ZeroFunction(),
+ boundary_values);
+ MatrixTools::apply_boundary_values (boundary_values,
+ system_matrix,
+ solution,
+ system_rhs);
+};
+
+
+
+template
+void LaplaceProblem::solve ()
+{
+ SolverControl solver_control (1000, 1e-12);
+ PrimitiveVectorMemory<> vector_memory;
+ SolverCG<> cg (solver_control, vector_memory);
+
+ // ...
+ PreconditionRelaxation<>
+ preconditioner(system_matrix,
+ &SparseMatrix::template precondition_SSOR,
+ 1.2);
+
+ cg.solve (system_matrix, solution, system_rhs,
+ preconditioner);
+
+ cout << " " << solver_control.last_step()
+ << " CG iterations needed to obtain convergence."
+ << endl;
+};
+
+
+
+template
+void LaplaceProblem::output_results ()
+{
+ DataOut data_out;
+
+ data_out.attach_dof_handler (dof_handler);
+ data_out.add_data_vector (solution, "solution");
+
+ data_out.build_patches ();
+
+ ofstream output (dim == 2 ?
+ "solution-2d.gmv" :
+ "solution-3d.gmv");
+ // ...
+ data_out.write_gnuplot (output);
+};
+
+
+
+template
+void LaplaceProblem::clear ()
+{
+ system_rhs.reinit (0);
+ solution.reinit (0);
+ system_matrix.reinit ();
+ sparsity_pattern.reinit (0, 0, 0);
+ dof_handler.clear ();
+ triangulation.clear ();
+};
+
+
+
+template
+void LaplaceProblem::run ()
+{
+ cout << "Solving problem in " << dim << " space dimensions." << endl;
+
+ for (unsigned int refinement=0; refinement<7; ++refinement)
+ {
+ cout << "Refinement step: " << refinement << endl;
+
+ make_grid_and_dofs(refinement);
+ assemble_system ();
+ solve ();
+ output_results ();
+
+ clear ();
+ };
+};
+
+
+
+int main ()
+{
+ deallog.depth_console (0);
+
+ LaplaceProblem<2> laplace_problem_2d;
+ laplace_problem_2d.run ();
+
+ return 0;
+};