From 5184b0ca1720493433c45e6e0692287218874e44 Mon Sep 17 00:00:00 2001 From: Timo Heister Date: Fri, 1 Feb 2013 23:19:03 +0000 Subject: [PATCH] starting step-49: more complex meshes git-svn-id: https://svn.dealii.org/trunk@28211 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-49/Makefile | 144 ++++++++++++++++ deal.II/examples/step-49/doc/builds-on | 1 + deal.II/examples/step-49/doc/intro.dox | 157 +++++++++++++++++ deal.II/examples/step-49/doc/kind | 1 + deal.II/examples/step-49/doc/results.dox | 79 +++++++++ deal.II/examples/step-49/doc/tooltip | 1 + deal.II/examples/step-49/step-49.cc | 209 +++++++++++++++++++++++ 7 files changed, 592 insertions(+) create mode 100644 deal.II/examples/step-49/Makefile create mode 100644 deal.II/examples/step-49/doc/builds-on create mode 100644 deal.II/examples/step-49/doc/intro.dox create mode 100644 deal.II/examples/step-49/doc/kind create mode 100644 deal.II/examples/step-49/doc/results.dox create mode 100644 deal.II/examples/step-49/doc/tooltip create mode 100644 deal.II/examples/step-49/step-49.cc diff --git a/deal.II/examples/step-49/Makefile b/deal.II/examples/step-49/Makefile new file mode 100644 index 0000000000..41c0fb67bc --- /dev/null +++ b/deal.II/examples/step-49/Makefile @@ -0,0 +1,144 @@ +# $Id: Makefile 25724 2012-07-24 23:35:36Z bangerth $ + + +# 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 = step-49 + +# 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 *vtk *ucd *.d2 + + + + +# +# +# Usually, you will not need to change anything beyond this point. +# +# +# The next statement tells 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. deal.II has two +# libraries: one for the debug mode version of the +# application and one for optimized mode. +libs.g := $(lib-deal2.g) +libs.o := $(lib-deal2.o) + + +# We now use the variable defined above to 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 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. +all: $(target)$(EXEEXT) +$(target)$(EXEEXT) : $(libraries) + @echo ============================ Linking $@ + @$(CXX) -o $@ $^ $(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)$(EXEEXT) + @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========= $( $@" + @$(CXX) $(CXXFLAGS.g) -c $< -o $@ +./%.$(OBJEXT) : + @echo "==============optimized===== $( $@" + @$(CXX) $(CXXFLAGS.o) -c $< -o $@ + + +# The following statement tells make that the rules `run' and `clean' +# are not expected to produce files of the same name as Makefile rules +# usually do. +.PHONY: all run clean + + +# Finally there is a rule which you normally need not care much about: +# since the executable depends on some include files from the library, +# besides the C++ application file of course, it is necessary to +# re-generate the executable when one of the files it depends on has +# changed. The following rule creates a dependency file +# `Makefile.dep', which `make' uses to determine when to regenerate +# the executable. 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. +# +# If the creation of Makefile.dep fails, blow it away and fail +Makefile.dep: $(target).cc Makefile \ + $(shell echo $D/include/deal.II/*/*.h) + @echo ============================ Remaking $@ + @$D/common/scripts/make_dependencies $(INCLUDE) -B. $(target).cc \ + > $@ \ + || (rm -f $@ ; false) + @if test -s $@ ; then true ; else rm $@ ; false ; fi + +# To make the dependencies known to `make', we finally have to include +# them: +include Makefile.dep + + diff --git a/deal.II/examples/step-49/doc/builds-on b/deal.II/examples/step-49/doc/builds-on new file mode 100644 index 0000000000..8b13789179 --- /dev/null +++ b/deal.II/examples/step-49/doc/builds-on @@ -0,0 +1 @@ + diff --git a/deal.II/examples/step-49/doc/intro.dox b/deal.II/examples/step-49/doc/intro.dox new file mode 100644 index 0000000000..c3629d2973 --- /dev/null +++ b/deal.II/examples/step-49/doc/intro.dox @@ -0,0 +1,157 @@ + +

Introduction

+ +

About the tutorial

+ +Since this is the first tutorial program, let us comment first on how +this tutorial and the rest of the deal.II documentation is supposed to +work. The documentation for deal.II comes essentially at three +different levels: +- The tutorial: This is a collection of programs that shows how + deal.II is used in practice. It doesn't typically discuss individual + functions at the level of individual arguments, but rather wants to + give the big picture of how things work together. In other words, it + discusses "concepts": what are the building blocks of deal.II and + how are they used together in finite element programs. +- The manual: This is the documentation of every single class and + every single (member) function in deal.II. You get there if, for + example, you click on the "Main page" or "Classes" tab at the top of + this page. This is the place where you would look up what the second + argument of Triangulation::create_triangulation_compatibility means, + to give just one slightly obscure example. You need this level of + documentation for when you know what you want to do, but forgot how + exactly the function was named, what its arguments are, or what it + returns. Note that you also get into the manual whenever you read + through the tutorial and click on any of the class or function + names, i.e. the tutorial contains a great many links into the manual + for whenever you need a more detailed description of a function or + class. On the other hand, the manual is not a good place to learn + deal.II since it gives you a microscopic view of things without + telling you how a function might fit into the bigger picture. +- Modules: These are groups of classes and functions that work + together or have related functionality. If you click on the + "Modules" tab at the top of this page, you end up on a page that + lists a number of such groups. Each module discusses the underlying + principles of these classes; for example, the @ref Sparsity module + talks about all sorts of different issues related to storing + sparsity patterns of matrices. This is documentation at an + intermediate level: they give you an overview of what's there in a + particular area. For example when you wonder what finite element + classes exist, you would take a look at the @ref fe module. The + modules are, of course, also cross-linked to the manual (and, at + times, to the tutorial); if you click on a class name, say on + Triangulation, would will also at the very top right under the class + name get a link to the modules this class is a member of if you want + to learn more about its context. + +Let's come back to the tutorial, since you are looking at the first program +(or "step") of it. Each tutorial program is subdivided into the following +sections: +
    +
  1. Introduction: This is a discussion of what the program + does, including the mathematical model, and + what programming techniques are new compared to previous + tutorial programs. +
  2. The commented program: An extensively documented listing of the + source code. Here, we often document individual lines, or + blocks of code, and discuss what they do, how they do it, and + why. The comments frequently reference the introduction, + i.e. you have to understand what the program wants to achieve + (a goal discussed in the introduction) before you can + understand how it intends to get there. +
  3. Results: The output of the program, with comments and + interpretation. This section also frequently has a subsection + that gives suggestions on how to extend the program in various + direction; in the earlier programs, this is intended to give + you directions for little experiments designed to make your + familiar with deal.II, while in later programs it is more about + how to use more advanced numerical techniques. +
  4. The plain program: The source code stripped of + all comments. This is useful if you want to see the "big + picture" of the code, since the commented version of the + program has so much text in between that it is often difficult + to see the entire code of a single function on the screen at + once. +
+ +The tutorials are not only meant to be static documentation, but you +should play with them. To this end, go to the +examples/step-1 directory (or whatever the number of the +tutorial is that you're interested in) and type +@code + make + make run +@endcode +The first command compiles the sources into an executable, while the +second executes it (strictly speaking, make run will also +compile the code if the executable doesn't exist yet, so you could +have skipped the first command if you wanted). This is all that's +needed to run the code and produce the output that is discussed in the +"Results" section of the tutorial programs. + +When learning the library, you need to play with it and see what +happens. To this end, open the examples/step-1/step-1.cc +source file with your favorite editor and modify it in some way, save it and +run it as above. A few suggestions for possibly modifications are given at the +end of the results section of this program, where we also provide a few links +to other useful pieces of information. + + +

What this program does

+ +Let's come back to step-1, the current program. +In this first example, we don't actually do very much, but show two +techniques: what is the syntax to generate triangulation objects, and +some elements of simple loops over all cells. We create two grids, one +which is a regularly refined square (not very exciting, but a common +starting grid for some problems), and one more geometric attempt: a +ring-shaped domain, which is refined towards the inner edge. Through +this, you will get to know three things every finite element program +will have to have somewhere: An object of type Triangulation for the +mesh; a call to the GridGenerator functions to generate a mesh; and +loops over all cells that involve iterators (iterators are a +generalization of pointers and are frequently used in the C++ standard +library; in the context of deal.II, the @ref Iterators module talks +about them). + +The program is otherwise small enough that it doesn't need a whole lot +of introduction. Let us just continue with its commented source. + + +

About scientific computing in general

+ +If you are reading through this tutorial program, chances are that you are +interested in continuing to use deal.II for your own projects. Thus, you are +about to embark on an exercise in programming using a large-scale scientific +computing library. Unless you are already an experienced user of large-scale +programming methods, this may be new territory for you — with all the +new rules that go along with it such as the fact that you will have to deal +with code written by others, that you may have to think about documenting your +own code because you may not remember what exactly it is doing a year down the +road (or because others will be using it as well), or coming up with ways to +test that your program is doing the right thing. None of this is something +that we typically train mathematicians, engineers, or scientists in but that +is important when you start writing software of more than a few hundred +lines. Remember: Producing software is not the same as just writing code. + +To make your life easier on this journey let us point to two resources that +are worthwhile browsing through before you start any large-scale programming: + +- The deal.II + Frequently Asked Questions: This page has a good number of questions + that pertain to particular aspects of deal.II, but also to more general + questions such as "How do I debug scientific computing codes?" or "Can I + train myself to write code that has fewer bugs?". + +- The Software Carpentry project + that provides introductions to many topics that are important to dealing + with software, such as version control, make files, testing, etc. It is + specifically written for scientists and engineers, not for computer + scientists, and has a focus on short, practical lessons. + +As a general recommendation: If you expect to spend more than a few days +writing software in the future, do yourself the favor of learning tools that +can make your life more productive, in particular debuggers and integrated +development environments. You will find that you will get the time spent +learning these tools back severalfold soon by being more productive! diff --git a/deal.II/examples/step-49/doc/kind b/deal.II/examples/step-49/doc/kind new file mode 100644 index 0000000000..15a13db451 --- /dev/null +++ b/deal.II/examples/step-49/doc/kind @@ -0,0 +1 @@ +basic diff --git a/deal.II/examples/step-49/doc/results.dox b/deal.II/examples/step-49/doc/results.dox new file mode 100644 index 0000000000..595a388fa6 --- /dev/null +++ b/deal.II/examples/step-49/doc/results.dox @@ -0,0 +1,79 @@ +

Results

+ +The program has, after having been run, produced two grids, which look +like this: + + + + + + + +
+ @image html step-1.grid-1.png + + @image html step-1.grid-2.png +
+ +The left one, well, is not very exciting. The right one is — at least +— unconventional. + +While the second mesh is entirely artificial and made-up, and +certainly not very practical in applications, to everyone's surprise it +has found its way into the literature: see the paper by M. Mu +titled "PDE.MART: A network-based problem-solving environment", ACM +Trans. Math. Software, vol. 31, pp. 508-531, 2005. Apparently it is +good for some things at least. + + +

Possible extensions

+ +

Different adaptive refinement strategies

+ +This program obviously does not have a whole lot of functionality, but +in particular the second_grid function has a bunch of +places where you can play with it. For example, you could modify the +criterion by which we decide which cells to refine. An example would +be to change the condition to this: +@code + for (; cell!=endc; ++cell) + if (cell->center()[1] > 0) + cell->set_refine_flag (); +@endcode +This would refine all cells for which the $y$-coordinate of the cell's +center is greater than zero (the TriaAccessor::center +function that we call by dereferencing the cell iterator +returns a Point<2> object; subscripting [0] would give +the $x$-coordinate, subscripting [1] the +$y$-coordinate). By looking at the functions that TriaAccessor +provides, you can also use more complicated criteria for refinement. + +

Different geometries

+ +Another possibility would be to generate meshes of entirely different +geometries altogether. While for complex geometries there is no way around +using meshes obtained from mesh generators, there is a good number of +geometries for which deal.II can create meshes using the functions in the +GridGenerator namespace. Take a look at what it provides and see how it could +be used in a program like this. + +

Comments about programming and debugging

+ +We close with a comment about modifying or writing programs with deal.II in +general. When you start working with tutorial programs or your own +applications, you will find that mistakes happen: your program will contain +code that either aborts the program right away or bugs that simply lead to +wrong results. In either case, you will find it extremely helpful to know how +to work with a debugger: you may get by for a while by just putting debug +output into your program, compiling it, and running it, but ultimately finding +bugs with a debugger is much faster, much more convenient, and more reliable +because you don't have to recompile the program all the time and because you +can inspect the values of variables and how they change. + +Rather than postponing learning how to use a debugger till you really can't +see any other way to find a bug, here's the one piece of +advice we will provide in this program: learn how to use a debugger as soon as +possible. It will be time well invested. The deal.II Frequently Asked +Questions (FAQ) page linked to from the top-level deal.II webpage also provides a good number +of hints on debugging deal.II programs. diff --git a/deal.II/examples/step-49/doc/tooltip b/deal.II/examples/step-49/doc/tooltip new file mode 100644 index 0000000000..71927bda1d --- /dev/null +++ b/deal.II/examples/step-49/doc/tooltip @@ -0,0 +1 @@ +Creating a grid. Refining it. Writing it to a file diff --git a/deal.II/examples/step-49/step-49.cc b/deal.II/examples/step-49/step-49.cc new file mode 100644 index 0000000000..2da4f71da2 --- /dev/null +++ b/deal.II/examples/step-49/step-49.cc @@ -0,0 +1,209 @@ +/* $Id: step-1.cc 27657 2012-11-21 13:19:08Z bangerth $ + * + * Copyright (C) 1999-2003, 2005-2007, 2009, 2011-2012 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 most fundamental class in the library is the Triangulation class, which +// is declared here: +#include +// We need the following two includes for loops over cells and/or faces: +#include +#include +// Here are some functions to generate standard grids: +#include +// We would like to use boundaries which are not straight lines, so we import +// some classes which predefine some boundary descriptions: +#include +// Output of grids in various graphics formats: +#include + +// This is needed for C++ output: +#include +// And this for the declarations of the `sqrt' and `fabs' functions: +#include + +// The final step in importing deal.II is this: All deal.II functions and +// classes are in a namespace dealii, to make sure they don't +// clash with symbols from other libraries you may want to use in conjunction +// with deal.II. One could use these functions and classes by prefixing every +// use of these names by dealii::, but that would quickly become +// cumbersome and annoying. Rather, we simply import the entire deal.II +// namespace for general use: +using namespace dealii; + +// @sect3{Creating the first mesh} + +// In the following, first function, we simply use the unit square as domain +// and produce a globally refined grid from it. +void first_grid () +{ + // The first thing to do is to define an object for a triangulation of a + // two-dimensional domain: + Triangulation<2> triangulation; + // Here and in many following cases, the string "<2>" after a class name + // indicates that this is an object that shall work in two space + // dimensions. Likewise, there are versions of the triangulation class that + // are working in one ("<1>") and three ("<3>") space dimensions. The way + // this works is through some template magic that we will investigate in + // some more detail in later example programs; there, we will also see how + // to write programs in an essentially dimension independent way. + + // Next, we want to fill the triangulation with a single cell for a square + // domain. The triangulation is the refined four times, to yield 4^4=256 + // cells in total: + GridGenerator::hyper_cube (triangulation); + triangulation.refine_global (4); + + // Now we want to write a graphical representation of the mesh to an output + // file. The GridOut class of deal.II can do that in a number of different + // output formats; here, we choose encapsulated postscript (eps) format: + std::ofstream out ("grid-1.eps"); + GridOut grid_out; + grid_out.write_eps (triangulation, out); +} + + + +// @sect3{Creating the second mesh} + +// The grid in the following, second function is slightly more complicated in +// that we use a ring domain and refine the result once globally. +void second_grid () +{ + // We start again by defining an object for a triangulation of a + // two-dimensional domain: + Triangulation<2> triangulation; + + // We then fill it with a ring domain. The center of the ring shall be the + // point (1,0), and inner and outer radius shall be 0.5 and 1. The number of + // circumferential cells could be adjusted automatically by this function, + // but we choose to set it explicitely to 10 as the last argument: + const Point<2> center (1,0); + const double inner_radius = 0.5, + outer_radius = 1.0; + GridGenerator::hyper_shell (triangulation, + center, inner_radius, outer_radius, + 10); + // By default, the triangulation assumes that all boundaries are straight + // and given by the cells of the coarse grid (which we just created). It + // uses this information when cells at the boundary are refined and new + // points need to be introduced on the boundary; if the boundary is assumed + // to be straight, then new points will simply be in the middle of the + // surrounding ones. + // + // Here, however, we would like to have a curved boundary. Fortunately, some + // good soul implemented an object which describes the boundary of a ring + // domain; it only needs the center of the ring and automatically figures + // out the inner and outer radius when needed. Note that we associate this + // boundary object with that part of the boundary that has the "boundary + // indicator" zero. By default (at least in 2d and 3d, the 1d case is + // slightly different), all boundary parts have this number, but you can + // change this number for some parts of the boundary. In that case, the + // curved boundary thus associated with number zero will not apply on those + // parts with a non-zero boundary indicator, but other boundary description + // objects can be associated with those non-zero indicators. If no boundary + // description is associated with a particular boundary indicator, a + // straight boundary is implied. + const HyperShellBoundary<2> boundary_description(center); + triangulation.set_boundary (0, boundary_description); + + // In order to demonstrate how to write a loop over all cells, we will + // refine the grid in five steps towards the inner circle of the domain: + for (unsigned int step=0; step<5; ++step) + { + // Next, we need an iterator which points to a cell and which we will + // move over all active cells one by one (active cells are those that + // are not further refined, and the only ones that can be marked for + // further refinement, obviously). By convention, we almost always use + // the names cell and endc for the iterator + // pointing to the present cell and to the one-past-the-end + // iterator: + Triangulation<2>::active_cell_iterator + cell = triangulation.begin_active(), + endc = triangulation.end(); + + // The loop over all cells is then rather trivial, and looks like any + // loop involving pointers instead of iterators: + for (; cell!=endc; ++cell) + // Next, we want to loop over all vertices of the cells. Since we are + // in 2d, we know that each cell has exactly four vertices. However, + // instead of penning down a 4 in the loop bound, we make a first + // attempt at writing it in a dimension-independent way by which we + // find out about the number of vertices of a cell. Using the + // GeometryInfo class, we will later have an easier time getting the + // program to also run in 3d: we only have to change all occurrences + // of <2> to <3>, and do not + // have to audit our code for the hidden appearance of magic numbers + // like a 4 that needs to be replaced by an 8: + for (unsigned int v=0; + v < GeometryInfo<2>::vertices_per_cell; + ++v) + { + // If this cell is at the inner boundary, then at least one of its + // vertices must sit on the inner ring and therefore have a radial + // distance from the center of exactly 0.5, up to floating point + // accuracy. Compute this distance, and if we have found a vertex + // with this property flag this cell for later refinement. We can + // then also break the loop over all vertices and move on to the + // next cell. + const double distance_from_center + = center.distance (cell->vertex(v)); + + if (std::fabs(distance_from_center - inner_radius) < 1e-10) + { + cell->set_refine_flag (); + break; + } + } + + // Now that we have marked all the cells that we want refined, we let + // the triangulation actually do this refinement. The function that does + // so owes its long name to the fact that one can also mark cells for + // coarsening, and the function does coarsening and refinement all at + // once: + triangulation.execute_coarsening_and_refinement (); + } + + + // Finally, after these five iterations of refinement, we want to again + // write the resulting mesh to a file, again in eps format. This works just + // as above: + std::ofstream out ("grid-2.eps"); + GridOut grid_out; + grid_out.write_eps (triangulation, out); + + + // At this point, all objects created in this function will be destroyed in + // reverse order. Unfortunately, we defined the boundary object after the + // triangulation, which still has a pointer to it and the library will + // produce an error if the boundary object is destroyed before the + // triangulation. We therefore have to release it, which can be done as + // follows. Note that this sets the boundary object used for part "0" of the + // boundary back to a default object, over which the triangulation has full + // control. + triangulation.set_boundary (0); + // An alternative to doing so, and one that is frequently more convenient, + // would have been to declare the boundary object before the triangulation + // object. In that case, the triangulation would have let lose of the + // boundary object upon its destruction, and everything would have been + // fine. +} + + + +// @sect3{The main function} + +// Finally, the main function. There isn't much to do here, only to call the +// two subfunctions, which produce the two grids. +int main () +{ + first_grid (); + second_grid (); +} -- 2.39.5