]> https://gitweb.dealii.org/ - code-gallery.git/commitdiff
A first test program. 1/head 4/head
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 10 Sep 2015 20:45:36 +0000 (15:45 -0500)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 10 Sep 2015 20:46:38 +0000 (15:46 -0500)
test-program/CMakeLists.txt [new file with mode: 0644]
test-program/doc/author [new file with mode: 0644]
test-program/doc/builds-on [new file with mode: 0644]
test-program/doc/tooltip [new file with mode: 0644]
test-program/step-1.cc [new file with mode: 0644]

diff --git a/test-program/CMakeLists.txt b/test-program/CMakeLists.txt
new file mode 100644 (file)
index 0000000..0dfbf45
--- /dev/null
@@ -0,0 +1,39 @@
+##
+#  CMake script for the step-1 tutorial program:
+##
+
+# Set the name of the project and target:
+SET(TARGET "step-1")
+
+# Declare all source files the target consists of. Here, this is only
+# the one step-X.cc file, but as you expand your project you may wish
+# to add other source files as well. If your project becomes much larger,
+# you may want to either replace the following statement by something like
+#    FILE(GLOB_RECURSE TARGET_SRC  "source/*.cc")
+#    FILE(GLOB_RECURSE TARGET_INC  "include/*.h")
+#    SET(TARGET_SRC ${TARGET_SRC}  ${TARGET_INC}) 
+# or switch altogether to the large project CMakeLists.txt file discussed
+# in the "CMake in user projects" page accessible from the "User info"
+# page of the documentation.
+SET(TARGET_SRC
+  ${TARGET}.cc
+  )
+
+# Usually, you will not need to modify anything beyond this point...
+
+CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8)
+
+FIND_PACKAGE(deal.II 8.4 QUIET
+  HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
+  )
+IF(NOT ${deal.II_FOUND})
+  MESSAGE(FATAL_ERROR "\n"
+    "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n"
+    "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n"
+    "or set an environment variable \"DEAL_II_DIR\" that contains this path."
+    )
+ENDIF()
+
+DEAL_II_INITIALIZE_CACHED_VARIABLES()
+PROJECT(${TARGET})
+DEAL_II_INVOKE_AUTOPILOT()
diff --git a/test-program/doc/author b/test-program/doc/author
new file mode 100644 (file)
index 0000000..f3c1eca
--- /dev/null
@@ -0,0 +1 @@
+Wolfgang Bangerth <bangerth@math.tamu.edu>
diff --git a/test-program/doc/builds-on b/test-program/doc/builds-on
new file mode 100644 (file)
index 0000000..8b13789
--- /dev/null
@@ -0,0 +1 @@
+
diff --git a/test-program/doc/tooltip b/test-program/doc/tooltip
new file mode 100644 (file)
index 0000000..a0f6101
--- /dev/null
@@ -0,0 +1 @@
+A test program
diff --git a/test-program/step-1.cc b/test-program/step-1.cc
new file mode 100644 (file)
index 0000000..42c437b
--- /dev/null
@@ -0,0 +1,264 @@
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 1999 - 2015 by the deal.II authors
+ *
+ * This file is part of the deal.II library.
+ *
+ * The deal.II library is free software; you can use it, redistribute
+ * it, and/or modify it under the terms of the GNU Lesser General
+ * Public License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ * The full text of the license can be found in the file LICENSE at
+ * the top level of the deal.II distribution.
+ *
+ * ---------------------------------------------------------------------
+
+ */
+
+// @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 faces and cells which are not straight lines,
+// or bi-linear quads, so we import some classes which predefine some
+// manifold descriptions:
+#include <deal.II/grid/manifold_lib.h>
+// Output of grids in various graphics formats:
+#include <deal.II/grid/grid_out.h>
+
+// This is needed for C++ output:
+#include <iostream>
+#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);
+  std::cout << "Grid written to grid-1.eps" << std::endl;
+}
+
+
+
+// @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 explicitly 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 lines, and all cells are bi-linear quads or tri-linear
+  // hexes, and that they are defined by the cells of the coarse grid
+  // (which we just created). Unless we do something special, when new
+  // points need to be introduced; the domain is assumed to be
+  // delineated by the straight lines of the coarse mesh, and new
+  // points will simply be in the middle of the surrounding ones.
+  // Here, however, we know that the domain is curved, and we would
+  // like to have the Triangulation place new points according to the
+  // underlying geometry. Fortunately, some good soul implemented an
+  // object which describes a spherical domain, of which the ring is a
+  // section; it only needs the center of the ring and automatically
+  // figures out how to instruct the Triangulation where to place the
+  // new points. The way this works in deal.II is that you tag parts
+  // of the triangulation you want to be curved with a number that is
+  // usually referred to as "manifold indicator" and then tell the
+  // triangulation to use a particular "manifold object" for all
+  // places with this manifold indicator. How exactly this works is
+  // not important at this point (you can read up on it in step-53 and
+  // @ref manifold). Here, for simplicity, we will choose the manifold
+  // id to be zero.  By default, all cells and faces of the
+  // Triangulation have their manifold_id set to
+  // numbers::invalid_manifold_id, which is the default if you want a
+  // manifold that produces straight edges, but you can change this
+  // number for individual cells and faces. In that case, the curved
+  // manifold thus associated with number zero will not apply to those
+  // parts with a non-zero manifold indicator, but other manifold
+  // description objects can be associated with those non-zero
+  // indicators. If no manifold description is associated with a
+  // particular manifold indicator, a manifold that produces straight
+  // edges is implied. (Manifold indicators are a slightly complicated
+  // topic; if you're confused about what exactly is happening here,
+  // you may want to look at the @ref GlossManifoldIndicator "glossary
+  // entry on this topic".)
+  triangulation.set_all_manifold_ids(0);
+  const SphericalManifold<2> manifold_description(center);
+  triangulation.set_manifold (0, manifold_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 that points to a cell and which we will
+      // move over all active cells one by one. In a sense, you can think of a
+      // triangulation as a collection of cells. If it was an array, you would
+      // just get a pointer that you move from one to the next. In
+      // triangulations, cells aren't stored as an array, so simple pointers
+      // do not work, but one can generalize pointers to iterators (see <a
+      // href="http://en.wikipedia.org/wiki/Iterator#C.2B.2B">this wikipedia
+      // link</a> for more information). We will then get an iterator to the
+      // first cell and iterate over all of the cells until we hit the last
+      // one.
+      //
+      // The second important piece is that we only need the active cells.
+      // Active cells are those that are not further refined, and the only
+      // ones that can be marked for further refinement, obviously. deal.II
+      // provides iterator categories that allow us to iterate over <i>all</i>
+      // cells (including the parent cells of active ones) or only over the
+      // active cells. Because we want the latter, we need to choose
+      // Triangulation::active_cell_iterator as data type.
+      //
+      // Finally, 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 "one-past-the-end" iterator. This is, in
+      // a sense a misnomer, because the object is not really a "cell": it is
+      // an iterator/pointer to a cell. We should really have started to call
+      // these objects <code>cell_iterator</code> when deal.II started in
+      // 1998, but it is what it is.
+      //
+      // After declaring the iterator variable, the loop over all cells is
+      // then rather trivial, and looks like any loop involving pointers
+      // instead of iterators:
+      Triangulation<2>::active_cell_iterator
+      cell = triangulation.begin_active(),
+      endc = triangulation.end();
+      for (; cell!=endc; ++cell)
+        {
+          // @note Writing a loop like this requires a lot of typing, but it
+          // is the only way of doing it in C++98 and C++03. However, if you
+          // have a C++11-compliant compiler, you can also use the C++11
+          // range-based for loop style that requires significantly less
+          // typing. Take a look at @ref CPP11 "the deal.II C++11 page" to see
+          // how this works.
+          //
+          // 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);
+
+  std::cout << "Grid written to grid-2.eps" << std::endl;
+
+  // At this point, all objects created in this function will be destroyed in
+  // reverse order. Unfortunately, we defined the manifold object after the
+  // triangulation, which still has a pointer to it and the library will
+  // produce an error if the manifold object is destroyed before the
+  // triangulation. We therefore have to release it, which can be done as
+  // follows. Note that this sets the manifold object used for part "0" of the
+  // domain back to a default object, over which the triangulation has full
+  // control.
+  triangulation.set_manifold (0);
+  // An alternative to doing so, and one that is frequently more convenient,
+  // would have been to declare the manifold object before the triangulation
+  // object. In that case, the triangulation would have let lose of the
+  // manifold 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.