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: +
    +
    +...
    +
    +
    + + +

    + + + + + +
    +solution-3d + +grid-3d +
    +

    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; +};