From: wolf
Date: Wed, 22 Dec 1999 21:28:11 +0000 (+0000)
Subject: Add step-4.
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=48c2798843f91aed31a115610e8f278ed175944c;p=dealii-svn.git
Add step-4.
git-svn-id: https://svn.dealii.org/trunk@2112 0785d39b-7218-0410-832d-ea1e28bc413d
---
diff --git a/deal.II/deal.II/Attic/examples/step-by-step/step-4/Makefile b/deal.II/deal.II/Attic/examples/step-by-step/step-4/Makefile
new file mode 100644
index 0000000000..9212069420
--- /dev/null
+++ b/deal.II/deal.II/Attic/examples/step-by-step/step-4/Makefile
@@ -0,0 +1,120 @@
+# $Id$
+# Copyright W. Bangerth, University of Heidelberg, 1998
+
+# Template for makefiles for the examples subdirectory. In principle,
+# everything should be done automatically if you set the target file
+# here correctly:
+target = step-4
+
+# 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
+
+# If you want your program to be linked with extra object or library
+# files, specify them here:
+user-libs =
+
+# To run the program, use "make run"; to give parameters to the program,
+# give the parameters to the following variable:
+run-parameters =
+
+# To execute additional action apart from running the program, fill
+# in this list:
+additional-run-action =
+
+# To specify which files are to be deleted by "make clean" (apart from
+# the usual ones: object files, executables, backups, etc), fill in the
+# following list
+delete-files = *gpl *inp *history
+
+
+
+
+###############################################################################
+# Internals
+
+#deal include base path
+D = ../../../..
+
+include $D/common/Make.global_options
+
+# get lists of files we need
+cc-files = $(filter-out *%, $(shell echo *.cc))
+o-files = $(cc-files:.cc=.o)
+go-files = $(cc-files:.cc=.go)
+h-files = $(filter-out *%, $(shell echo *.h))
+lib-h-files = $(filter-out *%, $(shell echo ../../include/*/*.h))
+
+# list of libraries needed to link with
+libs = ./Obj.a -ldeal_II_2d -llac -lbase
+libs.g = ./Obj.g.a -ldeal_II_3d.g -llac.g -lbase.g
+
+
+# check whether we use debug mode or not
+ifeq ($(debug-mode),on)
+libraries = $(libs.g)
+flags = $(CXXFLAGS.g)
+endif
+
+ifeq ($(debug-mode),off)
+libraries = $(libs)
+flags = $(CXXFLAGS)
+endif
+
+
+
+# make rule for the target
+$(target) : $(libraries)
+ @echo ============================ Linking $@
+ @$(CXX) $(flags) -o $@ $^ $(user-libs)
+
+# rule how to run the program
+run: $(target)
+ $(target) $(run-parameters)
+ $(additional-run-action)
+
+
+# rule to make object files
+%.go : %.cc
+ @echo ============================ Compiling with debugging information: $<
+ @echo $(CXX) ... -c $< -o $@
+ @$(CXX) $(CXXFLAGS.g) -c $< -o $@
+%.o : %.cc
+ @echo ============================ Compiling with optimization: $<
+ @echo $(CXX) ... -c $< -o $@
+ @$(CXX) $(CXXFLAGS) -c $< -o $@
+
+
+# rules which files the libraries depend upon
+Obj.a: ./Obj.a($(o-files))
+Obj.g.a: ./Obj.g.a($(go-files))
+
+
+clean:
+ -rm -f *.o *.go *~ Makefile.dep Obj.a Obj.g.a $(target) $(delete-files)
+
+
+
+.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.
+#
+#use perl to generate rules for the .go files as well
+#as to make rules not for tria.o and the like, but
+#rather for libnumerics.a(tria.o)
+Makefile.dep: $(cc-files) $(h-files) $(lib-h-files)
+ @echo ============================ Remaking Makefile
+ @perl $D/common/scripts/Make_dep.pl ./Obj $(INCLUDE) $(cc-files) \
+ > Makefile.dep
+
+
+include Makefile.dep
+
diff --git a/deal.II/deal.II/Attic/examples/step-by-step/step-4/step-4.cc b/deal.II/deal.II/Attic/examples/step-by-step/step-4/step-4.cc
new file mode 100644
index 0000000000..901bfa6e10
--- /dev/null
+++ b/deal.II/deal.II/Attic/examples/step-by-step/step-4/step-4.cc
@@ -0,0 +1,546 @@
+/* $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
+
+ // 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 line somewhere to the
+ // program; in this program, it was
+ // added to the main function.
+#include
+
+
+
+
+
+
+
+
+
+template
+class LaplaceProblem
+{
+ public:
+ LaplaceProblem ();
+ void run ();
+
+ private:
+ void make_grid_and_dofs ();
+ void assemble_system ();
+ void solve ();
+ void output_results ();
+
+ Triangulation triangulation;
+ FEQ1 fe;
+ DoFHandler dof_handler;
+
+ SparseMatrixStruct sparsity_pattern;
+ SparseMatrix system_matrix;
+
+ Vector solution;
+ Vector system_rhs;
+};
+
+
+ // In the following, we declare two
+ // more classes, which will represent
+ // the functions of the
+ // dim-dimensional space denoting the
+ // right hand side and the
+ // non-homogeneous Dirichlet boundary
+ // values.
+ //
+ // 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.
+template
+class RightHandSide : public Function
+{
+ public:
+ virtual double value (const Point &p,
+ const unsigned int component) const;
+};
+
+
+
+template
+class BoundaryValues : public Function
+{
+ public:
+ virtual double value (const Point &p,
+ const unsigned int component) const;
+};
+
+
+
+
+ // We wanted the right hand side
+ // function to be 4*(x**4+y**4) in
+ // 2D, or 4*(x**4+y**4+z**4) in
+ // 3D. Unfortunately, this is not as
+ // elegantly feasible dimension
+ // independently as much of the rest
+ // of this program, so we have to do
+ // it using a small
+ // loop. Fortunately, the compiler
+ // knows the size of the loop at
+ // compile time, i.e. the number of
+ // times the body will be executed,
+ // so it 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.
+ //
+ // Note that the different
+ // coordinates of the point are
+ // accessed using the () operator.
+template
+double RightHandSide::value (const Point &p,
+ const unsigned int) const
+{
+ double return_value = 0;
+ for (unsigned int i=0; i
+double BoundaryValues::value (const Point &p,
+ const unsigned int) const
+{
+ return p.square();
+};
+
+
+
+
+ // This is the constructor of the
+ // LaplaceProblem class. It
+ // associates the DoFHandler to the
+ // triangulation just as in the
+ // previous example.
+template
+LaplaceProblem::LaplaceProblem () :
+ dof_handler (triangulation)
+{};
+
+
+
+ // 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]x[-1,1] in 2D, or on
+ // the cube [-1,1]x[-1,1]x[-1,1] in
+ // 3D; both can be termed
+ // ``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.
+ //
+ // Likewise, associating a degree of
+ // freedom with each vertex is
+ // something which certainly looks
+ // different in 2D and 3D, but that
+ // does not need to bother you. This
+ // function therefore looks exactly
+ // like in the previous example,
+ // although it performs actions that
+ // in their details are quite
+ // different. The only significant
+ // difference is the number of cells
+ // resulting, which is much higher in
+ // three than in two space
+ // dimensions!
+template
+void LaplaceProblem::make_grid_and_dofs ()
+{
+ GridGenerator::hyper_cube (triangulation, -1, 1);
+ triangulation.refine_global (4);
+
+ 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());
+};
+
+
+
+ // 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 they
+ // way we assemble matrix and right
+ // hand side vector dimension
+ // independently: there is simply no
+ // difference to the pure
+ // 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.
+template
+void LaplaceProblem::assemble_system ()
+{
+ QGauss3 quadrature_formula;
+
+ // 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 in this
+ // function, we only declare it
+ // here, rather than as a member
+ // variable of the LaplaceProblem
+ // class, or somewhere else.
+ 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,
+ // they were only needed on the
+ // unit cell, in order to compute
+ // the values and gradients of the
+ // shape function, which are
+ // defined on the unit cell
+ // however). We can tell the
+ // FEValues object to do for us by
+ // giving it the update_q_points
+ // flag:
+ FEValues fe_values (fe, quadrature_formula,
+ UpdateFlags(update_values |
+ update_gradients |
+ update_q_points |
+ update_JxW_values));
+
+ // Note that the following numbers
+ // 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.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);
+
+ // Note here, 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 they are totally
+ // unrelated.
+ DoFHandler::active_cell_iterator cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+ for (; cell!=endc; ++cell)
+ {
+ fe_values.reinit (cell);
+ cell_matrix.clear ();
+ cell_rhs.clear ();
+
+ // 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; this makes
+ // things a bit faster.
+ 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,
+ BoundaryValues(),
+ boundary_values);
+ MatrixTools::apply_boundary_values (boundary_values,
+ system_matrix,
+ solution,
+ system_rhs);
+};
+
+
+ // Solving the linear system of
+ // equation is something that looks
+ // almost identical in most
+ // programs. In particular, it is
+ // dimension independent, so this
+ // function is mostly copied from the
+ // previous example.
+template
+void LaplaceProblem::solve ()
+{
+ SolverControl solver_control (1000, 1e-12);
+ PrimitiveVectorMemory<> vector_memory;
+ SolverCG<> cg (solver_control, vector_memory);
+ cg.solve (system_matrix, solution, system_rhs,
+ PreconditionIdentity());
+
+ // We have made one addition,
+ // though: since we suppress output
+ // from the linear solvers, we have
+ // to print the number of
+ // iterations by hand.
+ cout << " " << solver_control.last_step()
+ << " CG iterations needed to obtain convergence."
+ << endl;
+};
+
+
+
+ // This function also does what the
+ // respective one did in the previous
+ // example. No changes here for
+ // dimension independentce either.
+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 ();
+
+ // Only difference to the previous
+ // example: write output in GMV
+ // format, rather than for gnuplot.
+ ofstream output ("solution.gmv");
+ data_out.write_gmv (output);
+};
+
+
+
+ // 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.
+template
+void LaplaceProblem::run ()
+{
+ cout << "Solving problem in " << dim << " space dimensions." << endl;
+
+ make_grid_and_dofs();
+ assemble_system ();
+ solve ();
+ output_results ();
+};
+
+
+
+ // And this is the main function. It
+ // also looks mostly like in the
+ // previous example:
+int main ()
+{
+ // In the previous example, we had
+ // the output from the linear
+ // solvers about the starting
+ // residual and the number of the
+ // iteration where convergence was
+ // detected. This can be suppressed
+ // like this:
+ deallog.depth_console (0);
+ // The rationale here is the
+ // following: the deallog
+ // (i.e. deal-log, not de-allog)
+ // variable represents a stream to
+ // which some parts of the library
+ // write output. It redirects this
+ // output to the console and if
+ // required to a file. The output
+ // is nested in a way that each
+ // function can use a prefix string
+ // (separated by colons) for each
+ // line of output; if it calls
+ // another function, that may also
+ // use its prefix which is then
+ // printed after the one of the
+ // calling function. Since output
+ // from functions which are nested
+ // deep below is usually not as
+ // important as top-level output,
+ // you can give the deallog
+ // variable a maximal depth of
+ // nested output for output to
+ // console and file. The depth zero
+ // which we gave here means that no
+ // output is written.
+
+ // After having done this
+ // administrative stuff, we can go
+ // on just as before: define one of
+ // these top-level objects and
+ // transfer control to it:
+// LaplaceProblem<2> laplace_problem_2d;
+// laplace_problem_2d.run ();
+
+ LaplaceProblem<3> laplace_problem_3d;
+ laplace_problem_3d.run ();
+ return 0;
+};
diff --git a/deal.II/doc/tutorial/chapter-2.step-by-step/navbar.html b/deal.II/doc/tutorial/chapter-2.step-by-step/navbar.html
index 7307d29ed0..fbb14abf65 100644
--- a/deal.II/doc/tutorial/chapter-2.step-by-step/navbar.html
+++ b/deal.II/doc/tutorial/chapter-2.step-by-step/navbar.html
@@ -42,6 +42,12 @@
Step 3
+
+
+
+ Step 4
+
+
Back to the tutorial index
diff --git a/deal.II/doc/tutorial/chapter-2.step-by-step/step-4.data/grid-3d.jpg b/deal.II/doc/tutorial/chapter-2.step-by-step/step-4.data/grid-3d.jpg
new file mode 100644
index 0000000000..da8b54fe38
Binary files /dev/null and b/deal.II/doc/tutorial/chapter-2.step-by-step/step-4.data/grid-3d.jpg differ
diff --git a/deal.II/doc/tutorial/chapter-2.step-by-step/step-4.data/solution-3d.jpg b/deal.II/doc/tutorial/chapter-2.step-by-step/step-4.data/solution-3d.jpg
new file mode 100644
index 0000000000..d7db56da2b
Binary files /dev/null and b/deal.II/doc/tutorial/chapter-2.step-by-step/step-4.data/solution-3d.jpg differ
diff --git a/deal.II/doc/tutorial/chapter-2.step-by-step/step-4.intro b/deal.II/doc/tutorial/chapter-2.step-by-step/step-4.intro
new file mode 100644
index 0000000000..1518606e43
--- /dev/null
+++ b/deal.II/doc/tutorial/chapter-2.step-by-step/step-4.intro
@@ -0,0 +1,3 @@
+
+
Introduction
+
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
new file mode 100644
index 0000000000..4c30e00b53
--- /dev/null
+++ b/deal.II/doc/tutorial/chapter-2.step-by-step/step-4.results
@@ -0,0 +1,23 @@
+
+Results
+
+The output of the program looks as follows:
+
+
+...
+
+
+
+
+
+
+
+
+
+ |
+
+
+ |
+
+
+
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 2e05de8a0d..42b75b9e5c 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
@@ -57,6 +57,14 @@ At present, the following programs exist:
problem. Object-orientation. Assembling matrices and
vectors. Boundary values.
+
+
Step 4
+ What's new: This example is programmed in a
+ way that it is independent of the dimension for which we want to
+ solve Laplace's equation; we will solve the equation in 2D and
+ 3D, although the program is exactly the same. Non-constant right
+ hand side function. Non-homogeneous boundary values.
+
diff --git a/deal.II/examples/step-4/Makefile b/deal.II/examples/step-4/Makefile
new file mode 100644
index 0000000000..9212069420
--- /dev/null
+++ b/deal.II/examples/step-4/Makefile
@@ -0,0 +1,120 @@
+# $Id$
+# Copyright W. Bangerth, University of Heidelberg, 1998
+
+# Template for makefiles for the examples subdirectory. In principle,
+# everything should be done automatically if you set the target file
+# here correctly:
+target = step-4
+
+# 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
+
+# If you want your program to be linked with extra object or library
+# files, specify them here:
+user-libs =
+
+# To run the program, use "make run"; to give parameters to the program,
+# give the parameters to the following variable:
+run-parameters =
+
+# To execute additional action apart from running the program, fill
+# in this list:
+additional-run-action =
+
+# To specify which files are to be deleted by "make clean" (apart from
+# the usual ones: object files, executables, backups, etc), fill in the
+# following list
+delete-files = *gpl *inp *history
+
+
+
+
+###############################################################################
+# Internals
+
+#deal include base path
+D = ../../../..
+
+include $D/common/Make.global_options
+
+# get lists of files we need
+cc-files = $(filter-out *%, $(shell echo *.cc))
+o-files = $(cc-files:.cc=.o)
+go-files = $(cc-files:.cc=.go)
+h-files = $(filter-out *%, $(shell echo *.h))
+lib-h-files = $(filter-out *%, $(shell echo ../../include/*/*.h))
+
+# list of libraries needed to link with
+libs = ./Obj.a -ldeal_II_2d -llac -lbase
+libs.g = ./Obj.g.a -ldeal_II_3d.g -llac.g -lbase.g
+
+
+# check whether we use debug mode or not
+ifeq ($(debug-mode),on)
+libraries = $(libs.g)
+flags = $(CXXFLAGS.g)
+endif
+
+ifeq ($(debug-mode),off)
+libraries = $(libs)
+flags = $(CXXFLAGS)
+endif
+
+
+
+# make rule for the target
+$(target) : $(libraries)
+ @echo ============================ Linking $@
+ @$(CXX) $(flags) -o $@ $^ $(user-libs)
+
+# rule how to run the program
+run: $(target)
+ $(target) $(run-parameters)
+ $(additional-run-action)
+
+
+# rule to make object files
+%.go : %.cc
+ @echo ============================ Compiling with debugging information: $<
+ @echo $(CXX) ... -c $< -o $@
+ @$(CXX) $(CXXFLAGS.g) -c $< -o $@
+%.o : %.cc
+ @echo ============================ Compiling with optimization: $<
+ @echo $(CXX) ... -c $< -o $@
+ @$(CXX) $(CXXFLAGS) -c $< -o $@
+
+
+# rules which files the libraries depend upon
+Obj.a: ./Obj.a($(o-files))
+Obj.g.a: ./Obj.g.a($(go-files))
+
+
+clean:
+ -rm -f *.o *.go *~ Makefile.dep Obj.a Obj.g.a $(target) $(delete-files)
+
+
+
+.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.
+#
+#use perl to generate rules for the .go files as well
+#as to make rules not for tria.o and the like, but
+#rather for libnumerics.a(tria.o)
+Makefile.dep: $(cc-files) $(h-files) $(lib-h-files)
+ @echo ============================ Remaking Makefile
+ @perl $D/common/scripts/Make_dep.pl ./Obj $(INCLUDE) $(cc-files) \
+ > Makefile.dep
+
+
+include Makefile.dep
+
diff --git a/deal.II/examples/step-4/step-4.cc b/deal.II/examples/step-4/step-4.cc
new file mode 100644
index 0000000000..901bfa6e10
--- /dev/null
+++ b/deal.II/examples/step-4/step-4.cc
@@ -0,0 +1,546 @@
+/* $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
+
+ // 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 line somewhere to the
+ // program; in this program, it was
+ // added to the main function.
+#include
+
+
+
+
+
+
+
+
+
+template
+class LaplaceProblem
+{
+ public:
+ LaplaceProblem ();
+ void run ();
+
+ private:
+ void make_grid_and_dofs ();
+ void assemble_system ();
+ void solve ();
+ void output_results ();
+
+ Triangulation triangulation;
+ FEQ1 fe;
+ DoFHandler dof_handler;
+
+ SparseMatrixStruct sparsity_pattern;
+ SparseMatrix system_matrix;
+
+ Vector solution;
+ Vector system_rhs;
+};
+
+
+ // In the following, we declare two
+ // more classes, which will represent
+ // the functions of the
+ // dim-dimensional space denoting the
+ // right hand side and the
+ // non-homogeneous Dirichlet boundary
+ // values.
+ //
+ // 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.
+template
+class RightHandSide : public Function
+{
+ public:
+ virtual double value (const Point &p,
+ const unsigned int component) const;
+};
+
+
+
+template
+class BoundaryValues : public Function
+{
+ public:
+ virtual double value (const Point &p,
+ const unsigned int component) const;
+};
+
+
+
+
+ // We wanted the right hand side
+ // function to be 4*(x**4+y**4) in
+ // 2D, or 4*(x**4+y**4+z**4) in
+ // 3D. Unfortunately, this is not as
+ // elegantly feasible dimension
+ // independently as much of the rest
+ // of this program, so we have to do
+ // it using a small
+ // loop. Fortunately, the compiler
+ // knows the size of the loop at
+ // compile time, i.e. the number of
+ // times the body will be executed,
+ // so it 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.
+ //
+ // Note that the different
+ // coordinates of the point are
+ // accessed using the () operator.
+template
+double RightHandSide::value (const Point &p,
+ const unsigned int) const
+{
+ double return_value = 0;
+ for (unsigned int i=0; i
+double BoundaryValues::value (const Point &p,
+ const unsigned int) const
+{
+ return p.square();
+};
+
+
+
+
+ // This is the constructor of the
+ // LaplaceProblem class. It
+ // associates the DoFHandler to the
+ // triangulation just as in the
+ // previous example.
+template
+LaplaceProblem::LaplaceProblem () :
+ dof_handler (triangulation)
+{};
+
+
+
+ // 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]x[-1,1] in 2D, or on
+ // the cube [-1,1]x[-1,1]x[-1,1] in
+ // 3D; both can be termed
+ // ``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.
+ //
+ // Likewise, associating a degree of
+ // freedom with each vertex is
+ // something which certainly looks
+ // different in 2D and 3D, but that
+ // does not need to bother you. This
+ // function therefore looks exactly
+ // like in the previous example,
+ // although it performs actions that
+ // in their details are quite
+ // different. The only significant
+ // difference is the number of cells
+ // resulting, which is much higher in
+ // three than in two space
+ // dimensions!
+template
+void LaplaceProblem::make_grid_and_dofs ()
+{
+ GridGenerator::hyper_cube (triangulation, -1, 1);
+ triangulation.refine_global (4);
+
+ 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());
+};
+
+
+
+ // 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 they
+ // way we assemble matrix and right
+ // hand side vector dimension
+ // independently: there is simply no
+ // difference to the pure
+ // 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.
+template
+void LaplaceProblem::assemble_system ()
+{
+ QGauss3 quadrature_formula;
+
+ // 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 in this
+ // function, we only declare it
+ // here, rather than as a member
+ // variable of the LaplaceProblem
+ // class, or somewhere else.
+ 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,
+ // they were only needed on the
+ // unit cell, in order to compute
+ // the values and gradients of the
+ // shape function, which are
+ // defined on the unit cell
+ // however). We can tell the
+ // FEValues object to do for us by
+ // giving it the update_q_points
+ // flag:
+ FEValues fe_values (fe, quadrature_formula,
+ UpdateFlags(update_values |
+ update_gradients |
+ update_q_points |
+ update_JxW_values));
+
+ // Note that the following numbers
+ // 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.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);
+
+ // Note here, 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 they are totally
+ // unrelated.
+ DoFHandler::active_cell_iterator cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+ for (; cell!=endc; ++cell)
+ {
+ fe_values.reinit (cell);
+ cell_matrix.clear ();
+ cell_rhs.clear ();
+
+ // 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; this makes
+ // things a bit faster.
+ 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,
+ BoundaryValues(),
+ boundary_values);
+ MatrixTools::apply_boundary_values (boundary_values,
+ system_matrix,
+ solution,
+ system_rhs);
+};
+
+
+ // Solving the linear system of
+ // equation is something that looks
+ // almost identical in most
+ // programs. In particular, it is
+ // dimension independent, so this
+ // function is mostly copied from the
+ // previous example.
+template
+void LaplaceProblem::solve ()
+{
+ SolverControl solver_control (1000, 1e-12);
+ PrimitiveVectorMemory<> vector_memory;
+ SolverCG<> cg (solver_control, vector_memory);
+ cg.solve (system_matrix, solution, system_rhs,
+ PreconditionIdentity());
+
+ // We have made one addition,
+ // though: since we suppress output
+ // from the linear solvers, we have
+ // to print the number of
+ // iterations by hand.
+ cout << " " << solver_control.last_step()
+ << " CG iterations needed to obtain convergence."
+ << endl;
+};
+
+
+
+ // This function also does what the
+ // respective one did in the previous
+ // example. No changes here for
+ // dimension independentce either.
+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 ();
+
+ // Only difference to the previous
+ // example: write output in GMV
+ // format, rather than for gnuplot.
+ ofstream output ("solution.gmv");
+ data_out.write_gmv (output);
+};
+
+
+
+ // 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.
+template
+void LaplaceProblem::run ()
+{
+ cout << "Solving problem in " << dim << " space dimensions." << endl;
+
+ make_grid_and_dofs();
+ assemble_system ();
+ solve ();
+ output_results ();
+};
+
+
+
+ // And this is the main function. It
+ // also looks mostly like in the
+ // previous example:
+int main ()
+{
+ // In the previous example, we had
+ // the output from the linear
+ // solvers about the starting
+ // residual and the number of the
+ // iteration where convergence was
+ // detected. This can be suppressed
+ // like this:
+ deallog.depth_console (0);
+ // The rationale here is the
+ // following: the deallog
+ // (i.e. deal-log, not de-allog)
+ // variable represents a stream to
+ // which some parts of the library
+ // write output. It redirects this
+ // output to the console and if
+ // required to a file. The output
+ // is nested in a way that each
+ // function can use a prefix string
+ // (separated by colons) for each
+ // line of output; if it calls
+ // another function, that may also
+ // use its prefix which is then
+ // printed after the one of the
+ // calling function. Since output
+ // from functions which are nested
+ // deep below is usually not as
+ // important as top-level output,
+ // you can give the deallog
+ // variable a maximal depth of
+ // nested output for output to
+ // console and file. The depth zero
+ // which we gave here means that no
+ // output is written.
+
+ // After having done this
+ // administrative stuff, we can go
+ // on just as before: define one of
+ // these top-level objects and
+ // transfer control to it:
+// LaplaceProblem<2> laplace_problem_2d;
+// laplace_problem_2d.run ();
+
+ LaplaceProblem<3> laplace_problem_3d;
+ laplace_problem_3d.run ();
+ return 0;
+};