From: Wolfgang Bangerth Date: Wed, 19 Oct 2005 04:52:14 +0000 (+0000) Subject: Add earliest, dysfunctional draft of a step-20. X-Git-Tag: v8.0.0~12983 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=684fd2eb4520d339446e5a6c6f6485788e19703d;p=dealii.git Add earliest, dysfunctional draft of a step-20. git-svn-id: https://svn.dealii.org/trunk@11620 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-20/Makefile b/deal.II/examples/step-20/Makefile new file mode 100644 index 0000000000..ba429f5d42 --- /dev/null +++ b/deal.II/examples/step-20/Makefile @@ -0,0 +1,157 @@ +# $Id$ + + +# For the small projects Makefile, you basically need to fill in only +# four fields. +# +# The first is the name of the application. It is assumed that the +# application name is the same as the base file name of the single C++ +# file from which the application is generated. +target = $(basename $(shell echo step-*.cc)) + +# The second field determines whether you want to run your program in +# debug or optimized mode. The latter is significantly faster, but no +# run-time checking of parameters and internal states is performed, so +# you should set this value to `on' while you develop your program, +# and to `off' when running production computations. +debug-mode = on + + +# As third field, we need to give the path to the top-level deal.II +# directory. You need to adjust this to your needs. Since this path is +# probably the most often needed one in the Makefile internals, it is +# designated by a single-character variable, since that can be +# reference using $D only, i.e. without the parentheses that are +# required for most other parameters, as e.g. in $(target). +D = ../../ + + +# The last field specifies the names of data and other files that +# shall be deleted when calling `make clean'. Object and backup files, +# executables and the like are removed anyway. Here, we give a list of +# files in the various output formats that deal.II supports. +clean-up-files = *gmv *gnuplot *gpl *eps *pov + + + + +# +# +# Usually, you will not need to change something beyond this point. +# +# +# The next statement tell the `make' program where to find the +# deal.II top level directory and to include the file with the global +# settings +include $D/common/Make.global_options + + +# Since the whole project consists of only one file, we need not +# consider difficult dependencies. We only have to declare the +# libraries which we want to link to the object file, and there need +# to be two sets of libraries: one for the debug mode version of the +# application and one for the optimized mode. Here we have selected +# the versions for 2d. Note that the order in which the libraries are +# given here is important and that your applications won't link +# properly if they are given in another order. +# +# You may need to augment the lists of libraries when compiling your +# program for other dimensions, or when using third party libraries +libs.g = $(lib-deal2-2d.g) \ + $(lib-deal2-3d.g) \ + $(lib-lac.g) \ + $(lib-base.g) +libs.o = $(lib-deal2-2d.o) \ + $(lib-deal2-3d.o) \ + $(lib-lac.o) \ + $(lib-base.o) + + +# We now use the variable defined above which switch between debug and +# optimized mode to select the set of libraries to link with. Included +# in the list of libraries is the name of the object file which we +# will produce from the single C++ file. Note that by default we use +# the extension .g.o for object files compiled in debug mode and .o for +# object files in optimized mode (or whatever the local default on your +# system is instead of .o). +ifeq ($(debug-mode),on) + libraries = $(target).g.$(OBJEXT) $(libs.g) +else + libraries = $(target).$(OBJEXT) $(libs.o) +endif + + +# Now comes the first production rule: how to link the single object +# file produced from the single C++ file into the executable. Since +# this is the first rule in the Makefile, it is the one `make' selects +# if you call it without arguments. +$(target) : $(libraries) + @echo ============================ Linking $@ + @$(CXX) -o $@$(EXEEXT) $^ $(LIBS) $(LDFLAGS) + + +# To make running the application somewhat independent of the actual +# program name, we usually declare a rule `run' which simply runs the +# program. You can then run it by typing `make run'. This is also +# useful if you want to call the executable with arguments which do +# not change frequently. You may then want to add them to the +# following rule: +run: $(target) + @echo ============================ Running $< + @./$(target)$(EXEEXT) + + +# As a last rule to the `make' program, we define what to do when +# cleaning up a directory. This usually involves deleting object files +# and other automatically created files such as the executable itself, +# backup files, and data files. Since the latter are not usually quite +# diverse, you needed to declare them at the top of this file. +clean: + -rm -f *.$(OBJEXT) *~ Makefile.dep $(target)$(EXEEXT) $(clean-up-files) + + +# Since we have not yet stated how to make an object file from a C++ +# file, we should do so now. Since the many flags passed to the +# compiler are usually not of much interest, we suppress the actual +# command line using the `at' sign in the first column of the rules +# and write the string indicating what we do instead. +./%.g.$(OBJEXT) : + @echo ==============debug========= $( Makefile.dep + @if test -s $@ ; then : else rm $@ ; fi + + +# To make the dependencies known to `make', we finally have to include +# them: +include Makefile.dep + + diff --git a/deal.II/examples/step-20/step-20.cc b/deal.II/examples/step-20/step-20.cc new file mode 100644 index 0000000000..af5db06e1a --- /dev/null +++ b/deal.II/examples/step-20/step-20.cc @@ -0,0 +1,748 @@ +/* $Id$ */ +/* Author: Wolfgang Bangerth, University of Heidelberg, 1999 */ + +/* $Id$ */ +/* Version: $Name$ */ +/* */ +/* Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005 by the deal.II authors */ +/* */ +/* This file is subject to QPL and may not be distributed */ +/* without copyright and license information. Please refer */ +/* to the file deal.II/doc/license.html for the text and */ +/* further information on this license. */ + + // 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 +#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 + + + + // This is again the same + // LaplaceProblem class as in the + // previous example. The only + // difference is that we have now + // declared it as a class with a + // template parameter, and the + // template parameter is of course + // the spatial dimension in which we + // would like to solve the Laplace + // equation. Of course, several of + // the member variables depend on + // this dimension as well, in + // particular the Triangulation + // class, which has to represent + // quadrilaterals or hexahedra, + // respectively. Apart from this, + // everything is as before. +template +class LaplaceProblem +{ + public: + LaplaceProblem (); + void run (); + + private: + void make_grid_and_dofs (); + void assemble_system (); + void solve (); + void output_results () const; + + Triangulation triangulation; + FESystem fe; + DoFHandler dof_handler; + + BlockSparsityPattern sparsity_pattern; + BlockSparseMatrix system_matrix; + + BlockVector solution; + BlockVector 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. + // + // The `value' function takes a + // second argument, which we have + // here named `component': This is + // only meant for vector valued + // functions, where you may want to + // access a certain component of the + // vector at the point `p'. However, + // our functions are scalar, so we + // need not worry about this + // parameter and we will not use it + // in the implementation of the + // functions. Note that in the base + // class (Function), the declaration + // of the `value' function has a + // default value of zero for the + // component, so we will access the + // `value' function of the right hand + // side with only one parameter, + // namely the point where we want to + // evaluate the function. + // + // Note that the C++ language forces + // us to declare and define a + // constructor to the following + // classes even though they are + // empty. This is due to the fact + // that the base class has no default + // constructor (i.e. one without + // arguments), even though it has a + // constructor which has default + // values for all arguments. +template +class RightHandSide : public Function +{ + public: + RightHandSide () : Function() {}; + + virtual double value (const Point &p, + const unsigned int component = 0) const; +}; + + + +template +class BoundaryValues : public Function +{ + public: + BoundaryValues () : Function() {}; + + virtual double value (const Point &p, + const unsigned int component = 0) 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 (i.e. `x', `y', ...) + // of the point are accessed using + // the () operator. +template +double RightHandSide::value (const Point &p, + const unsigned int) const +{ + double return_value = deal_II_numbers::PI * deal_II_numbers::PI * dim; + 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 specifies + // the desired polynomial degree of + // the finite elements and associates + // the DoFHandler to the + // triangulation just as in the + // previous example. +template +LaplaceProblem::LaplaceProblem () : + fe (FE_RaviartThomas(2),1,FE_DGQ(2),1), + 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, 0, 1); + triangulation.refine_global (0); + + std::cout << " Number of active cells: " + << triangulation.n_active_cells() + << std::endl + << " Total number of cells: " + << triangulation.n_cells() + << std::endl; + + dof_handler.distribute_dofs (fe); + DoFRenumbering::component_wise (dof_handler); + + std::vector dofs_per_component (dim+1); + DoFTools::count_dofs_per_component (dof_handler, dofs_per_component); + const unsigned int n_u = dofs_per_component[0], + n_p = dofs_per_component[dim]; + + std::cout << " Number of degrees of freedom: " + << dof_handler.n_dofs() + << " (" << n_u << '+' << n_p << ')' + << std::endl; + + sparsity_pattern.reinit (2,2); + sparsity_pattern.block(0,0).reinit (n_u, n_u, + dof_handler.max_couplings_between_dofs()); + sparsity_pattern.block(1,0).reinit (n_p, n_u, + dof_handler.max_couplings_between_dofs()); + sparsity_pattern.block(0,1).reinit (n_u, n_p, + dof_handler.max_couplings_between_dofs()); + sparsity_pattern.block(1,1).reinit (n_p, n_p, + dof_handler.max_couplings_between_dofs()); + sparsity_pattern.collect_sizes(); + DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern); + sparsity_pattern.compress(); + + system_matrix.reinit (sparsity_pattern); + + std::vector block_components (2); + block_components[0] = n_u; + block_components[1] = n_p; + solution.reinit (block_components); + system_rhs.reinit (block_components); +} + + +Tensor<1,2> extract_u (const FEValues<2> &fe_values, + const unsigned int j, + const unsigned int q) +{ + Tensor<1,2> tmp; + tmp[0] = fe_values.shape_value_component (j,q,0); + tmp[1] = fe_values.shape_value_component (j,q,1); + return tmp; +} + + + +Tensor<1,3> extract_u (const FEValues<3> &fe_values, + const unsigned int j, + const unsigned int q) +{ + Tensor<1,3> tmp; + tmp[0] = fe_values.shape_value_component (j,q,0); + tmp[1] = fe_values.shape_value_component (j,q,1); + tmp[2] = fe_values.shape_value_component (j,q,2); + return tmp; +} + + + + + +double extract_div_u (const FEValues<2> &fe_values, + const unsigned int j, + const unsigned int q) +{ + return fe_values.shape_grad_component (j,q,0)[0] + + fe_values.shape_grad_component (j,q,1)[1]; +} + + +double extract_div_u (const FEValues<3> &fe_values, + const unsigned int j, + const unsigned int q) +{ + return fe_values.shape_grad_component (j,q,0)[0] + + fe_values.shape_grad_component (j,q,1)[1] + + fe_values.shape_grad_component (j,q,2)[2]; +} + + +template +double extract_p (const FEValues &fe_values, + const unsigned int j, + const unsigned int q) +{ + return fe_values.shape_value_component (j,q,dim); +} + + + + // 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 the + // 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 () +{ + QGauss quadrature_formula(2); + + // 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 local_matrix (dofs_per_cell, dofs_per_cell); + Vector local_rhs (dofs_per_cell); + std::vector rhs_values (n_q_points); + + std::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. + typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + { + fe_values.reinit (cell); + local_matrix = 0; + local_rhs = 0; + + right_hand_side.value_list (fe_values.get_quadrature_points(), + rhs_values); + + for (unsigned int q=0; q phi_i_u = extract_u (fe_values, i, q); + const double div_phi_i_u = extract_div_u (fe_values, i, q); + const double phi_i_p = extract_p (fe_values, i, q); + + for (unsigned int j=0; j phi_j_u = extract_u (fe_values, j, q); + const double div_phi_j_u = extract_div_u (fe_values, j, q); + const double phi_j_p = extract_p (fe_values, j, q); + + local_matrix(i,j) += (phi_i_u * phi_j_u + - div_phi_i_u * phi_j_p + + phi_i_p * div_phi_j_u) + * fe_values.JxW(q); + } + + local_rhs(i) += phi_i_p * + rhs_values[q] * + fe_values.JxW(q); + } + + cell->get_dof_indices (local_dof_indices); + for (unsigned int i=0; i &A) + : + A (A), + tmp1 (A.block(0,0).m()), + tmp2 (A.block(0,0).m()) + {} + + void vmult (Vector &dst, + const Vector &src) const + { + A.block(0,1).vmult (tmp1, src); + + SolverControl solver_control (tmp1.size(), + 1e-8*tmp1.l2_norm()); + PrimitiveVectorMemory<> vector_memory; + SolverGMRES<> cg (solver_control, vector_memory); + + A.block(0,0).print_formatted(std::cout, 2, false, 6, " ", 81); + FullMatrix F(24,24); + F.copy_from (A.block(0,0)); + std::cout << F.norm2() << ' ' << F.relative_symmetry_norm2() << std::endl; + + abort (); + + PreconditionSSOR<> precondition; + precondition.initialize(A.block(0,0)); + cg.solve (A.block(0,0), tmp2, tmp1, precondition); + + std::cout << " " << solver_control.last_step() + << " inner iterations needed to obtain convergence." + << std::endl; + + A.block(1,0).vmult (dst, tmp2); + + dst *= -1; + } + + private: + const BlockSparseMatrix &A; + + mutable Vector tmp1, tmp2; +}; + + + + // 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 (system_matrix.block(0,0).m(), + 1e-6*system_rhs.block(1).l2_norm()); + PrimitiveVectorMemory<> vector_memory; + SolverGMRES<> cg (solver_control, vector_memory); + + cg.solve (SchurComplement(system_matrix), solution.block(1), + system_rhs.block(1), + 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. + std::cout << " " << solver_control.last_step() + << " CG mass matrix iterations needed to obtain convergence." + << std::endl; + } + { + Vector tmp (system_matrix.block(0,0).m()); + system_matrix.block(0,1).vmult (tmp, solution.block(1)); + + SolverControl solver_control (system_matrix.block(0,0).m(), + 1e-6*tmp.l2_norm()); + PrimitiveVectorMemory<> vector_memory; + SolverGMRES<> cg (solver_control, vector_memory); + + cg.solve (system_matrix.block(0,0), solution.block(0), + tmp, 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. + std::cout << " " << solver_control.last_step() + << " CG Schur complement iterations needed to obtain convergence." + << std::endl; + } +} + + + + // This function also does what the + // respective one did in the previous + // example. No changes here for + // dimension independence either. +template +void LaplaceProblem::output_results () const +{ + 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. We use the dimension in + // the filename to generate + // distinct filenames for each run + // (in a better program, one would + // check whether `dim' can have + // other values than 2 or 3, but we + // neglect this here for the sake + // of brevity). + std::ofstream output (dim == 2 ? + "solution-2d" : + "solution-3d"); + data_out.write_gnuplot (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 () +{ + std::cout << "Solving problem in " << dim << " space dimensions." << std::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. Actually, now is the point + // where we have to tell the + // compiler which dimension we + // would like to use; all functions + // up to now including the classes + // were only templates and nothing + // has been compiled by now, but by + // declaring the following objects, + // the compiler will start to + // compile all the functions at the + // top using the template parameter + // replaced with a concrete value. + // + // For demonstration, we will first + // let the whole thing run in 2D + // and then in 3D: + LaplaceProblem<2> laplace_problem_2d; + laplace_problem_2d.run (); + +// LaplaceProblem<3> laplace_problem_3d; +// laplace_problem_3d.run (); + + return 0; +}