+<br>
+
+<i>This program was contributed by Wolfgang Bangerth and Luca Heltai, using
+data provided by D. Sarah Stamps.</i>
+
+@note This program elaborates on concepts of geometry and the classes that
+implement it. These classes are grouped into the documentation module on @ref
+manifold "Manifold description for triangulations". See there for additional
+information.
+
+
<a name="Intro"></a>
<h1>Introduction</h1>
+Partial differential equations for realistic problems are often posed on
+domains with complicated geometries. To provide just a few examples, consider
+these cases:
+- Among the two arguably most important industrial applications for the finite
+ element method, aerodynamics and more generally fluid dynamics is
+ one. Computer simulations today are used in the design of every airplane,
+ car, train and ship. The domain in which the partial differential equation
+ is posed is, in these cases, the air surrounding the plane with its wings,
+ flaps and engines; the air surrounding the car with its wheel, wheel wells,
+ mirrors and, in the case of race cars, all sorts of aerodynamic equipment;
+ the air surrounding the train with its wheels and gaps between cars. In the
+ case of ships, the domain is the water surrounding the ship with its rudders
+ and propellers.
+- The other of the two big applications of the finite element method is
+ structural engineering in which the domains are bridges, airplane necelles
+ and wings, and other solid bodies of often complicated shapes.
+- Finite element modeling is also often used to describe the generation and
+ propagation of earthquake waves. In these cases, one needs to accurately
+ represent the geometry of faults in the Earth crust. Since faults intersect,
+ dip at angles, and are often not completely straight, domains are frequently
+ very complex.
+One could cite many more examples of complicated geometries in which one wants
+to pose and solve a partial differential equation. What this shows is that the
+"real" world is much more complicated than what we have shown in almost all of
+the tutorial programs preceding this one.
+
+This program is therefore devoted to showing how one deals with complex
+geometries using a concrete application. In particular, what it shows is how
+we make a mesh fit the domain we want to solve on. On the other hand, what the
+program does not show is how to create a coarse for a domain. The process to
+arrive at a coarse mesh is called "mesh generation" and there are a number of
+high-quality programs that do this much better than we could ever
+implement. However, deal.II does have the ability to read in meshes in many
+formats generated by mesh generators and then make them fit a given shape,
+either by deforming a mesh or refining it a number of times until it fits. The
+deal.II Frequently Asked Questions page referenced from http://www.dealii.org/
+provides resources to mesh generators.
+
+
+<h2>Where geometry and meshes intersect</h2>
+
+Let us assume that you have a complex domain and that you already have a
+coarse mesh that somehow represents the general features of the domain. Then
+there are two situations in which it is necessary to describe to a deal.II
+program the details of your geometry:
+
+- Mesh refinement: Whenever a cell is refined, it is necessary to introduce
+ new vertices in the Triangulation. In the simplest case, one assumes that
+ the objects that make up the Triangulation are straight line segments, a
+ bi-linear surface or a tri-linear volume. The next vertex is then simply put
+ into the middle of the old ones. However, for curved boundaries or if we
+ want to solve a PDE on a curved, lower-dimensional manifold embedded in a
+ higher-dimensional space, this is insufficient since it will not respect the
+ actual geometry. We will therefore have to tell Triangulation where to put
+ new points.
+
+- Integration: When using higher order finite element methods, it is often
+ necessary to compute integrals using curved approximations of the boundary,
+ i.e., describe each edge or face of cells as curves, instead of straight
+ line segments or bilinear patches). The same is, of course, true when
+ integrating boundary terms (e.g., inhomogenous Neumann boundary
+ conditions). For the purpose of integration, the various Mapping classes
+ then provide the transformation from the reference cell to the actual cell.
+
+In both cases, we need a way to provide information about the geometry of the
+domain at the level of an individual cell, its faces and edges. This is where
+the Manifold class comes into play. Manifold is an abstract base class that
+only defines an interface by which the Triangulation and Mapping classes can
+query geometric information about the domain. Conceptually, Manifold sees the
+world in a way not dissimilar to how the mathematical subdiscipline geometry
+sees it: a domain is essentially just a collection of points that is somehow
+equipped with the notion of a distance between points so that we can obtain a
+point "in the middle" of some other points.
+
+deal.II provides a number of classes that implement the interface provided by
+Manifold for a variety of common geometries. On the other hand, in this
+program we will consider only a very common and much simpler case, namely the
+situation where (a part of) the domain we want to solve on can be described by
+transforming a much simpler domain (we will call this the "reference domain").
+In the language of mathematics, this means
+that the (part of the) domain is a <a
+href="http://en.wikipedia.org/wiki/Chart_%28topology%29">chart</a>. Charts are
+described by a smooth function that maps from the simpler domain to the chart
+(the "push-forward" function) and its inverse (the "pull-back" function). If
+the domain as a whole is not a chart (e.g., the surface of a sphere), then it
+can often be described as a collection of charts (e.g., the northern
+hemisphere and the southern hemisphere are each charts) and the domain can then
+be describe by an <a
+href="http://en.wikipedia.org/wiki/Atlas_%28topology%29">atlas</a>.
+
+If a domain can be decomposed into an atlas, all we need to do is provide the
+pull-back and push-forward functions for each of the charts. In deal.II, this
+means providing a class derived from ChartManifold, and this is precisely what
+we will do in this program.
+
+
+<h2>The example case</h2>
+
+To illustrate how one describes geometries using charts in deal.II, we will
+consider a case that originates in an application of the <a
+href="http://aspect.dealii.org">ASPECT mantle convection code</a>, using a
+data set provided by D. Sarah Stamps. In the concrete application, we were
+interested in describing flow in the Earth mantle under the <a
+href="http://en.wikipedia.org/wiki/East_African_rift">East African Rift</a>, a
+zone where two continental plates drift apart. Not to beat around the bush,
+the geometry we want to describe looks like this:
+
+<img src="http://www.dealii.org/images/steps/developer/step-53.topo.png" alt="">
+
+In particular, though you cannot see this here, the top surface is not
+just colored by the elevation but is, in fact, deformed to follow the
+correct topography.
+While the actual application is not relevant here, the geometry is. The domain
+we are interested in is a part of the Earth that ranges from the surface to a
+depth of 500km, from 26 to 35 degrees East of the Greenwich meridian, and from
+5 degrees North of the equator to 10 degrees South.
+
+This description of the geometry suggests to start with a box
+$\hat U=[26,35]\times[-10,5]\times[-500000,0]$ (measured in degrees,
+degrees, and meters) and to provide a map $\varphi$ so
+that $\varphi^{-1}(\hat U)=\Omega$ where $\Omega$ is the domain we
+seek. $(\Omega,\varphi)$ is then a chart, $\varphi$ the pull-back operator, and
+$\varphi^{-1}$ the push-forward operator. If we need a point $q$ that is the
+"average" of other points $q_i\in\Omega$, the ChartManifold class then first
+applies the pull-back to obtain $\hat q_i=\varphi(q_i)$, averages these to a
+point $\hat p$ and then computes $p=\varphi^{-1}(\hat p)$.
+
+Our goal here is therefore to implement a class that describes $\varphi$ and
+$\varphi^{-1}$. If Earth was a sphere, then this would not be difficult: if we
+denote by $(\hat \phi,\hat \theta,\hat d)$ the points of $\hat U$ (i.e.,
+longitude counted eastward, latitude counted northward, and elevation relative
+to zero depth), then
+@f[
+ \mathbf x = \varphi^{-1}(\hat \phi,\hat \theta,\hat d)
+ = (R+\hat d) (\cos\hat \phi\cos\hat \theta, \sin\hat \phi\cos\hat \theta, \sin\hat \theta)^T
+@f]
+provides coordinates in a Cartesian coordinate system, where $R$ is the radius
+of the sphere. However, the Earth is not a sphere:
+
+<ol>
+<li> It is flattened at the poles and larger at the equator: the semi-major axis
+ is approximately 22km longer than the semi-minor axis. We will account for
+ this using the <a href="http://en.wikipedia.org/wiki/WGS84">WGS 84</a>
+ reference standard for the Earth shape. The formula used in WGS 84 to obtain
+ a position in Cartesian coordinates from longitude, latitude, and elevation
+ is
+@f[
+ \mathbf x = \varphi_\text{WGS84}^{-1}(\phi,\theta,d)
+ = \left(
+ \begin{array}{c}
+ (\bar R(\theta)+d) \cos\phi\cos\theta, \\
+ (\bar R(\theta)+d) \sin\phi\cos\theta, \\
+ ((1-e^2)\bar R(\theta)+d) \sin\theta
+ \end{array}
+ \right),
+@f]
+ where $\bar R(\theta)=\frac{R}{\sqrt{1-(e \sin\theta)^2}}$, and radius and
+ ellipticity are given by $R=6378137\text{m}, e=0.081819190842622$. In this formula,
+ we assume that the arguments to sines and cosines are evaluated in degree, not
+ radians (though we will have to change this assumption in the code).
+
+<li> It has topography in the form of mountains and valleys. We will account for
+ this using topography data originally obtained from the US Geologic Survey
+ (USGS). Using this data set, we can look up elevations on a
+ latitude-longitude mesh laid over the surface of the Earth. Starting with
+ the box $\hat U=[26,35]\times[-10,5]\times[-500000,0]$, we will therefore
+ first stretch it in vertical direction before handing it off to the WGS 84
+ function: if $h(\hat\phi,\hat\theta)$ is the height at longitude $\hat\phi$
+ and latitude $\hat\theta$, then we define
+@f[
+ (\phi,\theta,d) =
+ \varphi_\text{topo}^{-1}(\hat\phi,\hat\theta,\hat d)
+ = \left(
+ \hat\phi,
+ \hat\theta,
+ \hat d + \frac{\hat d+500000}{500000}h(\hat\phi,\hat\theta)
+ \right).
+@f]
+ Using this function, the top surface of the box $\hat U$ is displaced to the
+ correct topography, the bottom surface remains where it was, and points in
+ between are linearly interpolated.
+</ol>
+
+Using these two functions, we can then define the entire push-forward function
+$\varphi^{-1}: \hat U \rightarrow \Omega$ as
+@f[
+ \mathbf x
+ =
+ \varphi^{-1}(\hat\phi,\hat\theta,\hat d)
+ =
+ \varphi_\text{WGS84}^{-1}(\varphi_\text{topo}^{-1}(\hat\phi,\hat\theta,\hat d)).
+@f]
+In addition, we will have to define the inverse of this function, the
+pull-back operation, which we can write as
+@f[
+ (\hat\phi,\hat\theta,\hat d)
+ =
+ \varphi(\mathbf x)
+ =
+ \varphi_\text{topo}(\varphi_\text{WGS84}(\mathbf x)).
+@f]
+We can obtain one of the components of this function by inverting the formula above:
+@f[
+ (\hat\phi,\hat\theta,\hat d) =
+ \varphi_\text{topo}(\phi,\theta,d)
+ = \left(
+ \phi,
+ \theta,
+ 500000\frac{d-h(\phi,\theta)}{500000+h(\phi,\theta)}
+ \right).
+@f]
+Computing $\varphi_\text{WGS84}(\mathbf x)$ is also possible though a lot more
+awkward. We won't show the formula here but instead only provide the implementation
+in the program.
+
+
+<h2>Implementation</h2>
+
+There are a number of issues we need to address in the program. At the largest scale,
+we need to write a class that implements the interface of ChartManifold. This involves
+a function <code>push_forward()</code> that takes a point
+in the reference domain $\hat U$ and transform it into real space using the function
+$\varphi^{-1}$ outlined above, and its inverse function <code>pull_back()</code>
+implementing $\varphi$. We will do so in the <code>AfricaGeometry</code> class below
+that looks, in essence, like this:
+@code
+ class AfricaGeometry : public ChartManifold<3,3>
+ {
+ public:
+ virtual
+ Point<3>
+ pull_back(const Point<3> &space_point) const;
+
+ virtual
+ Point<3>
+ push_forward(const Point<3> &chart_point) const;
+
+ private:
+ ... some member variables and other member functions...;
+ };
+@endcode
+
+The transformations above have two parts: the WGS 84 transformations and the topography
+transformation. Consequently, the <code>AfricaGeometry</code> class will have
+additional (non-virtual) member functions
+<code>AfricaGeometry::push_forward_wgs84()</code> and
+<code>AfricaGeometry::push_forward_topo()</code> that implement these two pieces, and
+corresponding pull back functions.
+
+The WGS 84 transformation functions are not particularly interesting (even though the
+formulas they implement are impressive). The more interesting part is the topography
+transformation. Recall that for this, we needed to evaluate the elevation function
+$h(\hat\phi,\hat\theta)$. There is of course no formula for this: Earth is what it is,
+the best one can do is look up the altitude from some table. This is, in fact what we
+will do.
+
+The data we use comes from the United States Geologic Survey (USGS) and was provided by
+D. Sarah Stamps who also wrote the initial version of the WGS 84 transformation functions.
+The topography data is stored in a file <code>topography.txt.gz</code> that, when unpackaed
+looks like this:
+@code
+6.983333 25.000000 700
+6.983333 25.016667 692
+6.983333 25.033333 701
+6.983333 25.050000 695
+6.983333 25.066667 710
+6.983333 25.083333 702
+...
+-11.983333 35.950000 707
+-11.983333 35.966667 687
+-11.983333 35.983333 659
+@endcode
+The data is formatted as <code>latitude longitude elevation</code> where the first two
+columns are provided in degrees North of the equator and degrees East of the Greenwich
+meridian. The final column is given in meters above the WGS 84 zero elevation.
+
+In the transformation functions, we need to evaluate $h(\hat\phi,\hat\theta)$ for a given
+longitude $\hat\phi$ and latitude $\hat\theta$. In general, this data point will not be
+available and we will have to interpolate between adjacent data points. Writing such an
+interpolation routine is not particularly difficult, but it is a bit tedious and error
+prone. Fortunately, we can somehow shoehorn this data set into an existing class:
+Functions::InterpolatedUniformGridData . Unfortunately, the class does not fit the bill
+quite exactly and so we need to work around it a bit. The problem comes from the way
+we initialize this class: in its simplest form, it takes a stream of values that it
+assumes form an equispaced mesh in the $x-y$ plane (or, here, the $\phi-\theta$ plane).
+Which is what they do here, sort of: they are ordered latitude first, longitude second;
+and more awkwardly, the first column starts at the largest values and counts down,
+rather than the usual other way around.
+
+Now, while tutorial programs are meant to illustrate how to code with deal.II, they do
+not necessarily have to satisfy the same quality standards as one would have to do
+with production codes. In a production code, we would write a function that reads the
+data and (i) automatically determines the extents of the first and second column,
+(ii) automatically determines the number of data points in each direction, (iii) does
+the interpolation regardless of the order in which data is arranged, if necessary
+by switching the order between reading and presenting it to the
+Functions::InterpolatedUniformGridData class.
+
+On the other hand, tutorial programs are best if they are short and demonstrate key
+points rather than dwell on unimportant aspects and, thereby, obscure what we really
+want to show. Consequently, we will allow ourselves a bit of leeway:
+- since this program is intended solely for a particular geometry around the area
+ of the East-African rift and since this is precisely the area described by the data
+ file, we will hardcode in the program that there are
+ $1139\times 660$ pieces of data;
+- we will hardcode the boundaries of the data
+ $[-11.98333^\circ,6.983333^\circ]\times[25^\circ,35.98333^\circ]$;
+- we will lie to the Functions::InterpolatedUniformGridData class: the class will
+ only see the data in the last column of this data file, and we will pretend that
+ the data is arranged in a way that there are 1139 data points in the first
+ coordinate direction that are arranged in <i>ascending</i> order but in an
+ interval $[-6.983333^\circ,11.98333^\circ]$ (not the negated bounds). Then,
+ when we need to look something up for a latitude $\hat\theta$, we can ask the
+ interpolating table class for a value at $-\hat\theta$. With this little
+ trick, we can avoid having to switch around the order of data as read from
+ file.
+
+All of this then calls for a class that essentially looks like this:
+@code
+ class AfricaTopography
+ {
+ public:
+ AfricaTopography ()
+ :
+ topography_data (...initialize somehow...)
+ {}
+
+ double value (const double lon, const double lat) const
+ {
+ return topography_data.value (Point<2>(-lat * 180/numbers::PI,
+ lon * 180/numbers::PI));
+ }
+
+ private:
+ const Functions::InterpolatedUniformGridData<2> topography_data;
+ };
+@endcode
+
+Note how the <code>value()</code> function negates the latitude. It also switches
+from the format $\phi,\theta$ that we use everywhere else to the latitude-longitude
+format used in the table. Finally, it takes its arguments in radians as that is what
+we do everywhere else in the program, but then converts them to the degree-based
+system used for table lookup. As you will see in the implementation below, the function
+has a few more (static) member functions that we will call in the initialization
+of the <code>topography_data</code> member variable: the class type of this variable
+has a constructor that allows us to set everything right at construction time,
+rather than having to fill data later on, but this constructor takes a number of
+objects that can't be constructed in-place (at least not in C++98). Consequently,
+the construction of each of the objects we want to pass in the initialization happens
+in a number of static member functions.
+
+Having discussed the general outline of how we want to implement things, let us go
+to the program and show how it is done in practice.
+
* ---------------------------------------------------------------------
*
- * Author: Wolfgang Bangerth, Texas A&M University, 2014
+ * Authors: Wolfgang Bangerth, Texas A&M University, 2014
+ * Luca Heltai, SISSA, 2014
+ * D. Sarah Stamps, MIT, 2014
*/
+// Let us start with the include files we need here. Obviously, we need the
+// ones that describe the triangulation (<code>tria.h</code>), and that allow
+// us to create and output triangulations (<code>grid_generator.h</code> and
+// <code>grid_out.h</code>). Furthermore, we need the header file that
+// declares the Manifold and ChartManifold classes that we will need to
+// describe the geometry (<code>manifold.h</code>). We will then also need
+// the GridTools::transform() function from the last of the following header
+// files; the purpose for this function will become discussed at the point
+// where we use it.
#include <deal.II/grid/tria.h>
-#include <deal.II/grid/tria_boundary.h>
-#include <deal.II/grid/tria_boundary_lib.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/manifold.h>
+#include <deal.II/grid/grid_tools.h>
-#include <deal.II/base/std_cxx11/array.h>
+// The remainder of the include files relate to reading the topography data.
+// As explained in the introduction, we will read it from a file and then
+// use the Functions::InterpolatedUniformGridData class that is declared in the
+// first of the following header files. Because the data is large, the file
+// we read from is stored as gzip compressed data and we make use of
+// some BOOST-provided functionality to read directly from gzipped data.
+#include <deal.II/base/function_lib.h>
-#include <fstream>
-
-
-using namespace dealii;
-
-
-template <int dim>
-class SphereGeometry : public Boundary<dim>
-{
-public:
- SphereGeometry (const Point<dim> ¢er);
- virtual
- Point<dim>
- get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const;
-
- virtual
- Point<dim>
- get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const;
-private:
+#include <boost/iostreams/filtering_stream.hpp>
+#include <boost/iostreams/filter/gzip.hpp>
+#include <boost/iostreams/device/file.hpp>
- std_cxx11::array<double,dim> pull_back (const Point<dim> &p) const;
- Point<dim> push_forward (const std_cxx11::array<double,dim> &preimage) const;
-
- template <int N>
- static std_cxx11::array<double,dim> average (const std_cxx11::array<double,dim> (&array)[N]);
-
- const Point<dim> center;
-};
-
-
-template <int dim>
-SphereGeometry<dim>::SphereGeometry (const Point<dim> ¢er)
- :
- center (center)
-{}
+#include <fstream>
-template <>
-std_cxx11::array<double,2>
-SphereGeometry<2>::pull_back (const Point<2> &p) const
+// The final part of the top matter is to open a namespace into which to put
+// everything, and then to import the dealii namespace into it.
+namespace Step53
{
- const Point<2> relative_point = p - center;
+ using namespace dealii;
- const double r = relative_point.norm();
- const double phi = std::atan2 (relative_point[1], relative_point[0]);
- std_cxx11::array<double,2> result;
- result[0] = r;
- result[1] = phi;
+ // @sect3{Describing topography: AfricaTopography}
+ //
+ // The first significant part of this program is the class that describes
+ // the topography $h(\hat phi,\hat \theta)$ as a function of longitude
+ // and latitude. As discussed in the introduction, we will make our life
+ // a bit easier here by not writing the class in the most general way
+ // possible but by only writing it for the particular purpose we are
+ // interested in here: interpolating data obtained from one very specific
+ // data file that contains information about a particular area of the
+ // world for which we know the extents.
+ //
+ // The general layout of the class has been discussed already above.
+ // Following is its declaration, including three static member functions
+ // that we will need in initializing the <code>topography_data</code>
+ // member variable.
+ class AfricaTopography
+ {
+ public:
+ AfricaTopography ();
+
+ double value (const double lon,
+ const double lat) const;
+
+ private:
+ const Functions::InterpolatedUniformGridData<2> topography_data;
+
+ static std_cxx1x::array<std::pair<double,double>,2> get_endpoints ();
+ static std_cxx1x::array<unsigned int,2> n_intervals ();
+ static std::vector<double> get_data ();
+ };
+
+
+ // Let us move to the implementation of the class. The interesting parts
+ // of the class are the constructor and the <code>value()</code> function.
+ // The former initializes the Functions::InterpolatedUniformGridData member
+ // variable and we will use the constructor that requires us to pass in
+ // the end points of the 2-dimensional data set we want to interpolate
+ // (which are here given by the intervals $[-6.983333, 11.98333]$,
+ // using the trick of switching end points discussed in the introduction,
+ // and $[25, 35.983333]$, both given in degrees), the number of intervals
+ // into which the data is split (1138 in latitude direction and 659 in
+ // longitude direction, for a total of $1139\times 660$ data points), and
+ // a Table object that contains the data. The data then of course has
+ // size $1139\times 660$ and we initialize it by providing an iterator
+ // to the first of the 751,740 elements of a std::vector object returned
+ // by the <code>get_data()</code> function below. Note that all of the
+ // member functions we call here are static because (i) they do not
+ // access any member variables of the class, and (ii) because they are
+ // called at a time when the object is not initialized fully anyway.
+ AfricaTopography::AfricaTopography ()
+ :
+ topography_data (get_endpoints(),
+ n_intervals(),
+ Table<2,double> (380, 220,
+ get_data().begin()))
+ {}
+
+
+ double
+ AfricaTopography::value (const double lon,
+ const double lat) const
+ {
+ return topography_data.value (Point<2>(-lat * 180/numbers::PI,
+ lon * 180/numbers::PI));
+ }
- return result;
-}
+ std_cxx1x::array<std::pair<double,double>,2>
+ AfricaTopography::get_endpoints ()
+ {
+ std_cxx1x::array<std::pair<double,double>,2> endpoints;
+ endpoints[0] = std::make_pair (-6.983333, 11.966667);
+ endpoints[1] = std::make_pair (25, 35.95);
+ return endpoints;
+ }
-template <>
-Point<2>
-SphereGeometry<2>::push_forward (const std_cxx11::array<double,2> &preimage) const
-{
- const Point<2> relative_point = preimage[0] * Point<2>(std::cos(preimage[1]), std::sin(preimage[1]));
- return relative_point + center;
-}
+ std_cxx1x::array<unsigned int,2>
+ AfricaTopography::n_intervals ()
+ {
+ std_cxx1x::array<unsigned int,2> endpoints;
+ endpoints[0] = 379;
+ endpoints[1] = 219;
+ return endpoints;
+ }
+ // The only other function of greater interest is the <code>get_data()</code>
+ // function. It returns a temporary vector that contains all 751,740 data
+ // points describing the altitude and is read from the file
+ // <code>topography.txt.gz</code>. Because the file is compressed by gzip,
+ // we cannot just read it through an object of type std::ifstream, but
+ // there are convenient methods in the BOOST library (see
+ // http://www.boost.org) that allows us to read from compressed files
+ // without first having to uncompress it on disk. The result is, basically,
+ // just another input stream that, for all practical purposes, looks just like
+ // the ones we always use.
+ //
+ // When reading the data, we read the three columns but throw ignore the
+ // first two. The datum in the last column is appended to an array that we
+ // the return and that will be copied into the table from which
+ // <code>topography_data</code> is initialized. Since the BOOST.iostreams
+ // library does not provide a very useful exception when the input file
+ // does not exist, is not readable, or does not contain the correct
+ // number of data lines, we catch all exceptions it may produce and
+ // create our own one. To this end, in the <code>catch</code>
+ // clause, we let the program run into an <code>AssertThrow(false, ...)</code>
+ // statement. Since the condition is always false, this always triggers an
+ // exception. In other words, this is equivalent to writing
+ // <code>throw ExcMessage("...")</code> but it also fills certain fields
+ // in the exception object that will later be printed on the screen
+ // identifying the function, file and line where the exception happened.
+ std::vector<double>
+ AfricaTopography::get_data ()
+ {
+ std::vector<double> data;
+
+ // create a stream where we read from gzipped data
+ boost::iostreams::filtering_istream in;
+ in.push(boost::iostreams::basic_gzip_decompressor<>());
+ in.push(boost::iostreams::file_source("topography.txt.gz"));
+
+ for (unsigned int line=0; line<751740; ++line)
+ {
+ try
+ {
+ double lat, lon, elevation;
+ in >> lat >> lon >> elevation;
+
+ data.push_back (elevation);
+ }
+ catch (...)
+ {
+ AssertThrow (false,
+ ExcMessage ("Could not read all 751,740 data points "
+ "from the file <topography.txt.gz>!"));
+ }
+ }
+
+ return data;
+ }
-template <>
-std_cxx11::array<double,3>
-SphereGeometry<3>::pull_back (const Point<3> &p) const
-{
- const Point<3> relative_point = p - center;
- const double r = relative_point.norm();
- const double phi = std::atan2 (relative_point[1], relative_point[0]);
- const double theta = std::atan2 (relative_point[2], std::sqrt (relative_point[0]*relative_point[0] +
- relative_point[1]*relative_point[1]));
+ // @sect3{Describing the geometry: AfricaGeometry}
+ //
+ // The following class is then the main one of this program. Its structure
+ // has been described in much detail in the introduction and does not need
+ // much introduction any more.
+ class AfricaGeometry : public ChartManifold<3,3>
+ {
+ public:
+ virtual
+ Point<3>
+ pull_back(const Point<3> &space_point) const;
- std_cxx11::array<double,3> result;
- result[0] = r;
- result[1] = phi;
- result[2] = theta;
+ virtual
+ Point<3>
+ push_forward(const Point<3> &chart_point) const;
- return result;
-}
+ private:
+ static const double R;
+ static const double ellipticity;
+ const AfricaTopography topography;
-template <>
-Point<3>
-SphereGeometry<3>::push_forward (const std_cxx11::array<double,3> &preimage) const
-{
- const Point<3> relative_point = (preimage[0] *
- Point<3>(std::cos(preimage[1]),
- std::sin(preimage[1]),
- std::cos(preimage[2])));
+ Point<3> push_forward_wgs84 (const Point<3> &phi_theta_d) const;
+ Point<3> pull_back_wgs84 (const Point<3> &x) const;
- return relative_point + center;
-}
+ Point<3> push_forward_topo (const Point<3> &phi_theta_d_hat) const;
+ Point<3> pull_back_topo (const Point<3> &phi_theta_d) const;
+ };
+ const double AfricaGeometry::R = 6378137;
+ const double AfricaGeometry::ellipticity = 8.1819190842622e-2;
-template <>
-template <int N>
-std_cxx11::array<double,2>
-SphereGeometry<2>::average (const std_cxx11::array<double,2> (&array)[N])
-{
- std_cxx11::array<double,2> result;
- // average the radii first. this is uncritical
+ // The implementation, as well, is pretty straightforward if you have
+ // read the introduction. In particular, both of the pull back and
+ // push forward functions are just concatenations of the respective
+ // functions of the WGS 84 and topography mappings:
+ Point<3>
+ AfricaGeometry::pull_back(const Point<3> &space_point) const
{
- result[0] = 0;
- for (unsigned int i=0; i<N; ++i)
- result[0] += array[i][0];
- result[0] /= N;
+ return pull_back_topo (pull_back_wgs84 (space_point));
}
- // now do the angle. there, we need to
- // be more careful because the average of 0.9*pi and -0.9*pi should not
- // be zero but pi. to this end, bring everything that is farther than
- // pi away from the first angle we want to average with, back to within pi
- // by adding/subtracting 2*pi
- //
- // we also want to make sure that we exclude from averaging all points
- // that lie at the origin since they have no angle at all
+ Point<3>
+ AfricaGeometry::push_forward(const Point<3> &chart_point) const
{
- bool origin_is_one_point = false;
-
- result[1] = 0;
- for (unsigned int i=0; i<N; ++i)
- if (array[i][0] > 1e-10)
- {
- const double angle = array[i][1];
- if (angle - array[0][1] > numbers::PI)
- result[1] += angle-2*numbers::PI;
- else if (angle - array[0][1] < -numbers::PI)
- result[1] += angle+2*numbers::PI;
- else
- result[1] += angle;
- }
- else
- origin_is_one_point = true;
-
- if (origin_is_one_point == false)
- result[1] /= N;
- else
- result[1] /= (N-1);
+ return push_forward_wgs84 (push_forward_topo (chart_point));
}
- return result;
-}
-
+ // The following two functions then define the forward and inverse
+ // transformations that correspond to the WGS 84 reference shape of
+ // Earth. The forward transform follows the formula shown in the
+ // introduction. The inverse transform is significantly more complicated
+ // and is, at the very least, not intuitive. It also suffers from the
+ // fact that it returns an angle that at the end of the function we
+ // need to clip back into the interval $[0,2\pi]$ if it should have
+ // escaped from there.
+ Point<3>
+ AfricaGeometry::push_forward_wgs84(const Point<3> &phi_theta_d) const
+ {
+ const double phi = phi_theta_d[0];
+ const double theta = phi_theta_d[1];
+ const double d = phi_theta_d[2];
-template <>
-template <int N>
-std_cxx11::array<double,3>
-SphereGeometry<3>::average (const std_cxx11::array<double,3> (&array)[N])
-{
- std_cxx11::array<double,3> result;
+ const double R_bar = R / std::sqrt(1 - (ellipticity * ellipticity *
+ std::sin(theta) * std::sin(theta)));
- // average the radii first. this is uncritical
- {
- result[0] = 0;
- for (unsigned int i=0; i<N; ++i)
- result[0] += array[i][0];
- result[0] /= N;
+ return Point<3> ((R_bar + d) * std::cos(phi) * std::cos(theta),
+ (R_bar + d) * std::sin(phi) * std::cos(theta),
+ ((1 - ellipticity * ellipticity) * R_bar + d) * std::sin(theta));
}
- // now do the angle along the equatorial direction. do the same as we did
- // in the 2d case
+ Point<3>
+ AfricaGeometry::pull_back_wgs84(const Point<3> &x) const
{
- bool origin_is_one_point = false;
-
- result[1] = 0;
- for (unsigned int i=0; i<N; ++i)
- if (array[i][0] > 1e-10)
- {
- const double angle = array[i][1];
- if (angle - array[0][1] > numbers::PI)
- result[1] += angle-2*numbers::PI;
- else if (angle - array[0][1] < -numbers::PI)
- result[1] += angle+2*numbers::PI;
- else
- result[1] += angle;
- }
- else
- origin_is_one_point = true;
-
- if (origin_is_one_point == false)
- result[1] /= N;
+ const double b = std::sqrt(R * R * (1 - ellipticity * ellipticity));
+ const double ep = std::sqrt((R * R - b * b) / (b * b));
+ const double p = std::sqrt(x(0) * x(0) + x(1) * x(1));
+ const double th = std::atan2(R * x(2), b * p);
+ const double phi = std::atan2(x(1), x(0));
+ const double theta = std::atan2(x(2) + ep * ep * b * std::pow(std::sin(th),3),
+ (p - (ellipticity * ellipticity * R * std::pow(std::cos(th),3))));
+ const double R_bar = R / (std::sqrt(1 - ellipticity * ellipticity * std::sin(theta) * std::sin(theta)));
+ const double R_plus_d = p / std::cos(theta);
+
+ Point<3> phi_theta_d;
+ if (phi < 0)
+ phi_theta_d[0] = phi + 2*numbers::PI;
+ else if (phi > 2*numbers::PI)
+ phi_theta_d[0] = phi - 2*numbers::PI;
else
- result[1] /= (N-1);
+ phi_theta_d[0] = phi;
+ phi_theta_d[1] = theta;
+ phi_theta_d[2] = R_plus_d - R_bar;
+ return phi_theta_d;
}
- // finally for the polar angle. the difficulty here is that we have, for
- // example, two points at $(r,\phi,\theta)$ and $(r,\phi+\pi,\pm \pi/2)$, then we want
- // to average these to $(r,\ast,\pi)$ where the equatorial angle does not matter
+
+ // In contrast, the topography transformations follow exactly the
+ // description in the introduction. There is not consequently not
+ // much to add:
+ Point<3>
+ AfricaGeometry::push_forward_topo(const Point<3> &phi_theta_d_hat) const
{
- //??? not sure how exactly to do this
+ const double d_hat = phi_theta_d_hat[2];
+ const double h = topography.value(phi_theta_d_hat[0],
+ phi_theta_d_hat[1]);
+ const double d = d_hat + (d_hat + 500000)/500000*h;
+ const Point<3> phi_theta_d (phi_theta_d_hat[0],
+ phi_theta_d_hat[1],
+ d);
+ return phi_theta_d;
}
+ Point<3>
+ AfricaGeometry::pull_back_topo(const Point<3> &phi_theta_d) const
+ {
+ const double d = phi_theta_d[2];
+ const double h = topography.value(phi_theta_d[0],
+ phi_theta_d[1]);
+ const double d_hat = 500000 * (d-h)/(500000+h);
+ const Point<3> phi_theta_d_hat (phi_theta_d[0],
+ phi_theta_d[1],
+ d_hat);
+ return phi_theta_d_hat;
+ }
return result;
}
for (unsigned int i=0; i<dim; ++i)
center[i] = .25;
const double radius=center.norm();
+=======
+>>>>>>> Replace the original attempt at providing step-53 by something real.
- SphereGeometry<dim> geometry(center);
- Triangulation<dim> triangulation;
- GridGenerator::hyper_cube (triangulation);
- triangulation.refine_global(1);
+ // @sect3{Creating the mesh}
+ //
+ // Having so described the properties of the geometry, not it is
+ // time to deal with the mesh used to discretize it. To this end,
+ // we create objects for the geometry and triangulation, and then
+ // proceed to create a $1\times 2\times 1$ rectangular mesh that
+ // corresponds to the reference domain
+ // $\hat U=[26,35]\times[-10,5]\times[-500000,0]$. We choose
+ // this number of subdivisions because it leads to cells that
+ // are roughly like cubes instead of stretched in one direction or
+ // another.
+ //
+ // Of course, we are not actually interested in meshing the
+ // reference domain. We are interested in meshing the real domain.
+ // Consequently, we will use the GridTools::transform() function
+ // that simply moves every point of a triangulation according to
+ // a given transformation. The transformation function it wants is
+ // a function that takes as its single argument a point in the reference
+ // domain and returns the corresponding location in the domain that we
+ // want to map to. This is, of course, exactly the push forward
+ // function of the geometry we use. However,
+ // <code>AfricaGeometry::push_forward()</code> requires two arguments:
+ // the <code>AfricaGeometry</code> object to work with via its implicit
+ // <code>this</code> pointer, and the point. We bind the first of these
+ // to the geometry object we have created at the top of the function
+ // and leave the second one open, obtaining the desired object to
+ // do the transformation.
+ void run ()
+ {
+ AfricaGeometry geometry;
+ Triangulation<3> triangulation;
- for (typename Triangulation<dim>::active_cell_iterator cell=triangulation.begin_active();
- cell!=triangulation.end(); ++cell)
{
- if (cell->center().distance(center)< radius)
- cell->set_manifold_id(1);
-
- for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
- if (cell->face(f)->center().distance(center)< radius)
- cell->face(f)->set_all_manifold_ids(1);
+ const Point<3> corner_points[2] = { Point<3>(26*numbers::PI/180,
+ -10*numbers::PI/180,
+ -500000),
+ Point<3>(35*numbers::PI/180,
+ 5*numbers::PI/180,
+ 0)
+ };
+ std::vector<unsigned int> subdivisions(3);
+ subdivisions[0] = 1;
+ subdivisions[1] = 2;
+ subdivisions[2] = 1;
+ GridGenerator::subdivided_hyper_rectangle (triangulation, subdivisions,
+ corner_points[0], corner_points[1],
+ true);
+
+ GridTools::transform (std_cxx1x::bind(&AfricaGeometry::push_forward,
+ std_cxx1x::cref(geometry),
+ std_cxx1x::_1),
+ triangulation);
}
- triangulation.set_manifold(1,geometry);
- triangulation.refine_global(1);
-
- const std::string filename = "mesh-" + Utilities::int_to_string(dim) + "d.vtk";
- std::ofstream out (filename.c_str());
- GridOut grid_out;
- grid_out.write_vtk (triangulation, out);
+ // The next step is to explain to the triangulation to use our geometry
+ // object whenever a new point is needed upon refining the mesh. We do
+ // this by telling the triangulation to use our geometry for everythin
+ // that has manifold indicator zero, and then proceed to mark all cells
+ // and their bounding faces and edges with manifold indicator zero. This
+ // ensures that the triangulation consults our geometry object everytime
+ // a new vertex is needed. Since manifold indicators are inherited from
+ // mother to children, this also happens after several recursive
+ // refinement steps.
+ triangulation.set_manifold(0, geometry);
+ for (Triangulation<3>::active_cell_iterator cell=triangulation.begin_active();
+ cell!=triangulation.end(); ++cell)
+ cell->set_all_manifold_ids(0);
+
+ // The last step is to refine the mesh beyond its initial $1\times 2\times 1$
+ // coarse mesh. We could just refine globally a number of times, but since for
+ // the purpose of this tutorial program we're really only interested in what
+ // is happening close to the surface, we just refine 6 times all of the cells
+ // that have a face at a boundary with indicator 5. Looking this up in the
+ // documentation of the GridGenerator::subdivided_hyper_rectangle() function
+ // we have used above reveals that boundary indicator 5 corresponds to the top
+ // surface of the domain (and this is what the last <code>true</code> argument
+ // in the call to GridGenerator::subdivided_hyper_rectangle() above meant: to
+ // "color" the boundaries by assigning each boundary a unique boundary indicator).
+ for (unsigned int i=0; i<6; ++i)
+ {
+ for (Triangulation<3>::active_cell_iterator cell=triangulation.begin_active();
+ cell!=triangulation.end(); ++cell)
+ for (unsigned int f=0; f<GeometryInfo<3>::faces_per_cell; ++f)
+ if (cell->face(f)->boundary_indicator() == 5)
+ {
+ cell->set_refine_flag ();
+ break;
+ }
+ triangulation.execute_coarsening_and_refinement();
+
+ std::cout << "Refinement step " << i+1 << ": "
+ << triangulation.n_active_cells() << " cells, "
+ << GridTools::minimal_cell_diameter (triangulation)/1000
+ << "km minimal cell diameter"
+ << std::endl;
+ }
+
+ // Having done this all, we can now output the mesh into a file of its own:
+ const std::string filename = "mesh.vtu";
+ std::ofstream out (filename.c_str());
+ GridOut grid_out;
+ grid_out.write_vtu (triangulation, out);
+ }
}
-
// @sect3{The main function}
-// Finally, the main function. There isn't much to do here, only to call the
-// subfunctions.
+// Finally, the main function, which follows the same scheme used in all
+// tutorial programs starting with step-6. There isn't much to do here, only
+// to call the single <code>run()</code> function.
int main ()
{
- make_grid<2> ();
-// make_grid<3> ();
+ try
+ {
+ Step53::run ();
+ }
+ catch (std::exception &exc)
+ {
+ std::cerr << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ std::cerr << "Exception on processing: " << std::endl
+ << exc.what() << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+
+ return 1;
+ }
+ catch (...)
+ {
+ std::cerr << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ std::cerr << "Unknown exception!" << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ }
}
+
projects / dealii.git / commitdiff