]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
starting step-49: more complex meshes
authorheister <heister@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 1 Feb 2013 23:19:03 +0000 (23:19 +0000)
committerheister <heister@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 1 Feb 2013 23:19:03 +0000 (23:19 +0000)
git-svn-id: https://svn.dealii.org/trunk@28211 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-49/Makefile [new file with mode: 0644]
deal.II/examples/step-49/doc/builds-on [new file with mode: 0644]
deal.II/examples/step-49/doc/intro.dox [new file with mode: 0644]
deal.II/examples/step-49/doc/kind [new file with mode: 0644]
deal.II/examples/step-49/doc/results.dox [new file with mode: 0644]
deal.II/examples/step-49/doc/step-1.grid-1.png [new file with mode: 0644]
deal.II/examples/step-49/doc/step-1.grid-2.png [new file with mode: 0644]
deal.II/examples/step-49/doc/tooltip [new file with mode: 0644]
deal.II/examples/step-49/step-49.cc [new file with mode: 0644]

diff --git a/deal.II/examples/step-49/Makefile b/deal.II/examples/step-49/Makefile
new file mode 100644 (file)
index 0000000..41c0fb6
--- /dev/null
@@ -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========= $(<F)  ->  $@"
+       @$(CXX) $(CXXFLAGS.g) -c $< -o $@
+./%.$(OBJEXT) :
+       @echo "==============optimized===== $(<F)  ->  $@"
+       @$(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 (file)
index 0000000..8b13789
--- /dev/null
@@ -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 (file)
index 0000000..c3629d2
--- /dev/null
@@ -0,0 +1,157 @@
+<a name="Intro"></a>
+<h1>Introduction</h1>
+
+<h3> About the tutorial </h3>
+
+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:
+<ol>
+  <li> <b>Introduction:</b> This is a discussion of what the program
+       does, including the mathematical model, and
+       what programming techniques are new compared to previous
+       tutorial programs.
+  <li> <b>The commented program:</b> 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 <i>what</i> the program wants to achieve
+       (a goal discussed in the introduction) before you can
+       understand <i>how</i> it intends to get there.
+  <li> <b>Results:</b> 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.
+  <li> <b>The plain program:</b> 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.
+</ol>
+
+The tutorials are not only meant to be static documentation, but you
+should play with them. To this end, go to the
+<code>examples/step-1</code> 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, <code>make run</code> 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 <code>examples/step-1/step-1.cc</code>
+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.
+
+
+<h3> What this program does </h3>
+
+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.
+
+
+<h3> About scientific computing in general </h3>
+
+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 &mdash; 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 <a
+  href="http://dealii.sourceforge.net/index.php/Deal.II_Questions_and_Answers">deal.II
+  Frequently Asked Questions</a>: 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 <a href="http://software-carpentry.org/">Software Carpentry project</a>
+  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 (file)
index 0000000..15a13db
--- /dev/null
@@ -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 (file)
index 0000000..595a388
--- /dev/null
@@ -0,0 +1,79 @@
+<h1>Results</h1>
+
+The program has, after having been run, produced two grids, which look
+like this:
+
+<TABLE WIDTH="60%" ALIGN="center">
+  <tr>
+    <td ALIGN="center">
+      @image html step-1.grid-1.png
+    </td>
+
+    <td ALIGN="center">
+      @image html step-1.grid-2.png
+    </td>
+  </tr>
+</table>
+
+The left one, well, is not very exciting. The right one is &mdash; at least
+&mdash; 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.
+
+
+<h3> Possible extensions </h3>
+
+<h4> Different adaptive refinement strategies </h4>
+
+This program obviously does not have a whole lot of functionality, but
+in particular the <code>second_grid</code> 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 <code>TriaAccessor::center</code>
+function that we call by dereferencing the <code>cell</code> iterator
+returns a Point<2> object; subscripting <code>[0]</code> would give
+the $x$-coordinate, subscripting <code>[1]</code> the
+$y$-coordinate). By looking at the functions that TriaAccessor
+provides, you can also use more complicated criteria for refinement.
+
+<h4> Different geometries </h4>
+
+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.
+
+<h4> Comments about programming and debugging </h4>
+
+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 <a
+href="http://www.dealii.org/">deal.II webpage</a> also provides a good number
+of hints on debugging deal.II programs.
diff --git a/deal.II/examples/step-49/doc/step-1.grid-1.png b/deal.II/examples/step-49/doc/step-1.grid-1.png
new file mode 100644 (file)
index 0000000..e7edc25
Binary files /dev/null and b/deal.II/examples/step-49/doc/step-1.grid-1.png differ
diff --git a/deal.II/examples/step-49/doc/step-1.grid-2.png b/deal.II/examples/step-49/doc/step-1.grid-2.png
new file mode 100644 (file)
index 0000000..8dc1d31
Binary files /dev/null and b/deal.II/examples/step-49/doc/step-1.grid-2.png differ
diff --git a/deal.II/examples/step-49/doc/tooltip b/deal.II/examples/step-49/doc/tooltip
new file mode 100644 (file)
index 0000000..71927bd
--- /dev/null
@@ -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 (file)
index 0000000..2da4f71
--- /dev/null
@@ -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 <deal.II/grid/tria.h>
+// We need the following two includes for loops over cells and/or faces:
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+// Here are some functions to generate standard grids:
+#include <deal.II/grid/grid_generator.h>
+// We would like to use boundaries which are not straight lines, so we import
+// some classes which predefine some boundary descriptions:
+#include <deal.II/grid/tria_boundary_lib.h>
+// Output of grids in various graphics formats:
+#include <deal.II/grid/grid_out.h>
+
+// This is needed for C++ output:
+#include <fstream>
+// And this for the declarations of the `sqrt' and `fabs' functions:
+#include <cmath>
+
+// The final step in importing deal.II is this: All deal.II functions and
+// classes are in a namespace <code>dealii</code>, 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 <code>dealii::</code>, 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 <code>cell</code> and <code>endc</code> for the iterator
+      // pointing to the present cell and to the <code>one-past-the-end</code>
+      // 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 <code>&lt;2&gt;</code> to <code>&lt;3&gt;</code>, 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 ();
+}

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.