From: Wolfgang Bangerth Date: Thu, 10 Sep 2015 20:45:36 +0000 (-0500) Subject: A first test program. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=250c96e02adc078f3285f26ffa6aacde4958fe85;p=code-gallery.git A first test program. --- diff --git a/test-program/CMakeLists.txt b/test-program/CMakeLists.txt new file mode 100644 index 0000000..0dfbf45 --- /dev/null +++ b/test-program/CMakeLists.txt @@ -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 index 0000000..f3c1eca --- /dev/null +++ b/test-program/doc/author @@ -0,0 +1 @@ +Wolfgang Bangerth diff --git a/test-program/doc/builds-on b/test-program/doc/builds-on new file mode 100644 index 0000000..8b13789 --- /dev/null +++ b/test-program/doc/builds-on @@ -0,0 +1 @@ + diff --git a/test-program/doc/tooltip b/test-program/doc/tooltip new file mode 100644 index 0000000..a0f6101 --- /dev/null +++ b/test-program/doc/tooltip @@ -0,0 +1 @@ +A test program diff --git a/test-program/step-1.cc b/test-program/step-1.cc new file mode 100644 index 0000000..42c437b --- /dev/null +++ b/test-program/step-1.cc @@ -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 +// 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 faces and cells which are not straight lines, +// or bi-linear quads, so we import some classes which predefine some +// manifold descriptions: +#include +// Output of grids in various graphics formats: +#include + +// This is needed for C++ output: +#include +#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); + 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 this wikipedia + // link 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 all + // 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 + // cell and endc 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 cell_iterator 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 <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); + + 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 (); +}