From a7d426ece473d19c6ac33ae1447c8e9ed4d41463 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 4 Sep 2006 02:39:18 +0000 Subject: [PATCH] Add step-26 git-svn-id: https://svn.dealii.org/trunk@13829 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-26/Makefile | 154 ++++++ deal.II/examples/step-26/step-26.cc | 759 ++++++++++++++++++++++++++++ 2 files changed, 913 insertions(+) create mode 100644 deal.II/examples/step-26/Makefile create mode 100644 deal.II/examples/step-26/step-26.cc diff --git a/deal.II/examples/step-26/Makefile b/deal.II/examples/step-26/Makefile new file mode 100644 index 0000000000..7f74cb9645 --- /dev/null +++ b/deal.II/examples/step-26/Makefile @@ -0,0 +1,154 @@ +# $Id: Makefile 11779 2005-11-23 15:54:25Z wolf $ + + +# 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 = off + + +# 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-3d.g) \ + $(lib-lac.g) \ + $(lib-base.g) +libs.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========= $( $@ \ + || (rm -f $@ ; false) + @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-26/step-26.cc b/deal.II/examples/step-26/step-26.cc new file mode 100644 index 0000000000..0cffe5165b --- /dev/null +++ b/deal.II/examples/step-26/step-26.cc @@ -0,0 +1,759 @@ +/* $Id: step-4.cc 13345 2006-07-06 13:44:43Z guido $ */ +/* Author: Wolfgang Bangerth, University of Heidelberg, 1999 */ + +/* $Id: step-4.cc 13345 2006-07-06 13:44:43Z guido $ */ +/* Version: $Name$ */ +/* */ +/* Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006 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. */ + + // @sect3{Include files} + + // 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 + // single line somewhere to the program (see + // the main() function below for that): +#include + + +#include +#include +#include + + +class PointDefinedSurface : public StraightBoundary<3> +{ + public: + PointDefinedSurface (const std::string &filename); + + Point<3> closest_point (const Point<3> &p) const; + + /** + * Let the new point be the + * arithmetic mean of the two + * vertices of the line. + * + * Refer to the general + * documentation of this class + * and the documentation of the + * base class for more + * information. + */ + virtual Point<3> + get_new_point_on_line (const Triangulation<3>::line_iterator &line) const; + + /** + * Let the new point be the + * arithmetic mean of the four + * vertices of this quad and the + * four midpoints of the lines, + * which are already created at + * the time of calling this + * function. + * + * Refer to the general + * documentation of this class + * and the documentation of the + * base class for more + * information. + */ + virtual Point<3> + get_new_point_on_quad (const Triangulation<3>::quad_iterator &quad) const; + + /** + * Gives n=points.size() + * points that splits the + * p{StraightBoundary} line into + * p{n+1} partitions of equal + * lengths. + * + * Refer to the general + * documentation of this class + * and the documentation of the + * base class. + */ + virtual void + get_intermediate_points_on_line (const Triangulation<3>::line_iterator &line, + std::vector > &points) const; + + /** + * Gives n=points.size()=m*m + * points that splits the + * p{StraightBoundary} quad into + * (m+1)(m+1) subquads of equal + * size. + * + * Refer to the general + * documentation of this class + * and the documentation of the + * base class. + */ + virtual void + get_intermediate_points_on_quad (const Triangulation<3>::quad_iterator &quad, + std::vector > &points) const; + private: + std::vector > point_list; +}; + + +PointDefinedSurface::PointDefinedSurface (const std::string &filename) +{ + // first read in all the points + { + std::ifstream in (filename.c_str()); + AssertThrow (in, ExcIO()); + + while (in) + { + Point<3> p; + in >> p; + point_list.push_back (p); + } + + AssertThrow (point_list.size() > 1, ExcIO()); + } + + // next fit a linear model through the data + // cloud to rectify it in a local + // coordinate system + // + // the first step is to move the center of + // mass of the points to the origin + { + const Point<3> c_o_m = std::accumulate (point_list.begin(), + point_list.end(), + Point<3>()) / + point_list.size(); + for (unsigned int i=0; i gradient_direction + = Point<2>(a,b) / std::sqrt(a*a+b*b); + const Point<2> orthogonal_direction + = Point<2>(-b,a) / std::sqrt(a*a+b*b); + + const double stretch_factor = std::sqrt(1.+a*a+b*b); + + for (unsigned int i=0; i xy (point_list[i][0], + point_list[i][1]); + const double grad_distance = xy * gradient_direction; + const double orth_distance = xy * orthogonal_direction; + + // we then have to stretch the points + // in the gradient direction. the + // stretch factor is defined above + // (zero if the original plane was + // already the xy plane, infinity if + // it was vertical) + const Point<2> new_xy + = (grad_distance * stretch_factor * gradient_direction + + orth_distance * orthogonal_direction); + point_list[i][0] = new_xy[0]; + point_list[i][1] = new_xy[1]; + } + } +} + + +Point<3> +PointDefinedSurface::closest_point (const Point<3> &p) const +{ + double distance = p.distance (point_list[0]); + Point<3> point = point_list[0]; + + for (std::vector >::const_iterator i=point_list.begin(); + i != point_list.end(); ++i) + { + const double d = p.distance (*i); + if (d < distance) + { + distance = d; + point = *i; + } + } + + return point; +} + + +Point<3> +PointDefinedSurface:: +get_new_point_on_line (const Triangulation<3>::line_iterator &line) const +{ + return closest_point (StraightBoundary<3>::get_new_point_on_line (line)); +} + + + +Point<3> +PointDefinedSurface:: +get_new_point_on_quad (const Triangulation<3>::quad_iterator &quad) const +{ + return closest_point (StraightBoundary<3>::get_new_point_on_quad (quad)); +} + + + +void +PointDefinedSurface:: +get_intermediate_points_on_line (const Triangulation<3>::line_iterator &line, + std::vector > &points) const +{ + StraightBoundary<3>::get_intermediate_points_on_line (line, + points); + for (unsigned int i=0; i::quad_iterator &quad, + std::vector > &points) const +{ + StraightBoundary<3>::get_intermediate_points_on_quad (quad, + points); + for (unsigned int i=0; iLaplaceProblem class template} + + // 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; + FE_Q fe; + DoFHandler dof_handler; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; +}; + + + // @sect3{Right hand side and boundary values} + + + + +template +class BoundaryValues : public Function +{ + public: + BoundaryValues () : Function() {}; + + virtual double value (const Point &p, + const unsigned int component = 0) const; +}; + + + +template +double BoundaryValues::value (const Point &p, + const unsigned int /*component*/) const +{ + return std::max(p[dim-1], -5.); +} + + + + // @sect3{Implementation of the LaplaceProblem class} + + // Next for the implementation of the class + // template that makes use of the functions + // above. As before, we will write everything + // as templates that have a formal parameter + // dim that we assume unknown at the time + // we define the template functions. Only + // later, the compiler will find a + // declaration of LaplaceProblem@<2@> (in + // the main function, actually) and + // compile the entire class with dim + // replaced by 2, a process referred to as + // `instantiation of a template'. When doing + // so, it will also replace instances of + // RightHandSide@ by + // RightHandSide@<2@> and instantiate the + // latter class from the class template. + // + // In fact, the compiler will also find a + // declaration LaplaceProblem@<3@> in + // main(). This will cause it to again go + // back to the general + // LaplaceProblem@ template, replace + // all occurrences of dim, this time by + // 3, and compile the class a second + // time. Note that the two instantiations + // LaplaceProblem@<2@> and + // LaplaceProblem@<3@> are completely + // independent classes; their only common + // feature is that they are both instantiated + // from the same general template, but they + // are not convertible into each other, for + // example, and share no code (both + // instantiations are compiled completely + // independently). + + + // @sect4{LaplaceProblem::LaplaceProblem} + + // After this introduction, here 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 program, step-3: +template +LaplaceProblem::LaplaceProblem () : + fe (1), + dof_handler (triangulation) +{} + + + // @sect4{LaplaceProblem::make_grid_and_dofs} + + // 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 + // either. This function therefore looks + // exactly like in the previous example, + // although it performs actions that in their + // details are quite different if dim + // happens to be 3. The only significant + // difference from a user's perspective 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, -30, 30); + + for (unsigned int f=0; f::faces_per_cell; ++f) + if (triangulation.begin()->face(f)->center()[2] > 15) + { + triangulation.begin()->face(f)->set_boundary_indicator (1); + for (unsigned int i=0; i::lines_per_face; ++i) + triangulation.begin()->face(f)->line(i)->set_boundary_indicator (1); + break; + } + triangulation.set_boundary (1, pds); + + + for (unsigned int v=0; v::vertices_per_cell; ++v) + if (triangulation.begin()->vertex(v)[2] > 0) + triangulation.begin()->vertex(v) + = pds.closest_point (Point<3>(triangulation.begin()->vertex(v)[0], + triangulation.begin()->vertex(v)[1], + 0)); + + for (unsigned int i=0; i<7; ++i) + { + for (typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + cell != triangulation.end(); ++cell) + for (unsigned int f=0; f::faces_per_cell; ++f) + if (cell->face(f)->boundary_indicator() == 1) + cell->set_refine_flag (); + + triangulation.execute_coarsening_and_refinement (); + + std::cout << "Refinement cycle " << i << std::endl + << " Number of active cells: " + << triangulation.n_active_cells() + << std::endl + << " Total number of cells: " + << triangulation.n_cells() + << std::endl; + + } + + + dof_handler.distribute_dofs (fe); + + std::cout << " Number of degrees of freedom: " + << dof_handler.n_dofs() + << std::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()); +} + + + // @sect4{LaplaceProblem::assemble_system} + + // 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 + // 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 () +{ + MatrixTools::create_laplace_matrix (dof_handler, + QGauss(2), + system_matrix); + system_rhs = 0; + + std::map boundary_values; + VectorTools::interpolate_boundary_values (dof_handler, + 0, + BoundaryValues(), + boundary_values); + MatrixTools::apply_boundary_values (boundary_values, + system_matrix, + solution, + system_rhs); +} + + + // @sect4{LaplaceProblem::solve} + + // Solving the linear system of + // equations is something that looks + // almost identical in most + // programs. In particular, it is + // dimension independent, so this + // function is copied verbatim from the + // previous example. +template +void LaplaceProblem::solve () +{ + SolverControl solver_control (1000, 1e-12); + SolverCG<> cg (solver_control); + + PreconditionSSOR<> preconditioner; + preconditioner.initialize(system_matrix, 1.2); + + cg.solve (system_matrix, solution, system_rhs, + preconditioner); +} + + + // @sect4{LaplaceProblem::output_results} + + // This function also does what the + // respective one did in step-3. No changes + // here for dimension independence either. + // + // The only difference to the previous + // example is that we want to write output in + // GMV format, rather than for gnuplot (GMV + // is another graphics program that, contrary + // to gnuplot, shows data in nice colors, + // allows rotation of geometries with the + // mouse, and generates reasonable + // representations of 3d data; for ways to + // obtain it see the ReadMe file of + // deal.II). To write data in this format, we + // simply replace the + // data_out.write_gnuplot call by + // data_out.write_gmv. + // + // Since the program will run both 2d and 3d + // versions of the laplace solver, 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). +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 (); + + std::ofstream output (dim == 2 ? + "solution-2d.gmv" : + "solution-3d.gmv"); + data_out.write_gmv (output); +} + + + + // @sect4{LaplaceProblem::run} + + // 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 (); +} + + + // @sect3{The main function} + + // And this is the main function. It also + // looks mostly like in step-3, but if you + // look at the code below, note how we first + // create a variable of type + // LaplaceProblem@<2@> (forcing the + // compiler to compile the class template + // with dim replaced by 2) and run a + // 2d simulation, and then we do the whole + // thing over in 3d. + // + // In practice, this is probably not what you + // would do very frequently (you probably + // either want to solve a 2d problem, or one + // in 3d, but not both at the same + // time). However, it demonstrates the + // mechanism by which we can simply change + // which dimension we want in a single place, + // and thereby force the compiler to + // recompile the dimension independent class + // templates for the dimension we + // request. The emphasis here lies on the + // fact that we only need to change a single + // place. This makes it rather trivial to + // debug the program in 2d where computations + // are fast, and then switch a single place + // to a 3 to run the much more computing + // intensive program in 3d for `real' + // computations. + // + // Each of the two blocks is enclosed in + // braces to make sure that the + // laplace_problem_2d variable goes out + // of scope (and releases the memory it + // holds) before we move on to allocate + // memory for the 3d case. Without the + // additional braces, the + // laplace_problem_2d variable would only + // be destroyed at the end of the function, + // i.e. after running the 3d problem, and + // would needlessly hog memory while the 3d + // run could actually use it. + // + // Finally, the first line of the function is + // used to suppress some output. Remember + // that 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 through + // the deallog.depth_console(0) call. + // + // 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 so 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. By changing it you can get more + // information about the innards of the + // library. +int main () +{ + deallog.depth_console (0); + { + LaplaceProblem<3> laplace_problem_3d; + laplace_problem_3d.run (); + } + + return 0; +} -- 2.39.5