From: Wolfgang Bangerth Date: Mon, 25 Aug 2014 22:26:40 +0000 (-0500) Subject: Replace the original attempt at providing step-53 by something real. X-Git-Tag: v8.2.0-rc1~169^2~5 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c76982038e5d528bf1b596f7329fd708038dfe68;p=dealii.git Replace the original attempt at providing step-53 by something real. In particular, provide a step-53 that shows how to describe geometries using charts and use a section of Earth under Eastern Africa as an example. This has about the right level of complexity, and produces awesome pictures. This commit was generated from a number of smaller one: More text. Implement the Africa model for the WGS 84 geometry (actual topography still missing). Provide the actual data set. Provide a way of reading it when not compressed. Results and finish the program. Much still to be done. More text. Finish documenting everything. Remove a picture we now no longer need. --- diff --git a/doc/doxygen/headers/manifold.h b/doc/doxygen/headers/manifold.h index 84cfbe3b01..36666f111b 100644 --- a/doc/doxygen/headers/manifold.h +++ b/doc/doxygen/headers/manifold.h @@ -75,9 +75,10 @@ * is the ChartManifold class, and this is the class that users will * likely overload for complex geometries. * - * While this process is non trivial in most cases of interest, for - * most of the trivial geometries, like cylinders, spheres or shells, - * we provide reasonable implementations. + * While this process is non trivial in most cases of interest, for most of + * the trivial geometries, like cylinders, spheres or shells, deal.II provides + * reasonable implementations. More complicated examples can be described + * using the techniques shown in step-53. * * The boundary of a Triangulation is a special case of Manifold, for * which additional information can be useful in user codes, such as diff --git a/examples/step-53/doc/intro.dox b/examples/step-53/doc/intro.dox index 1518606e43..a3da00ad63 100644 --- a/examples/step-53/doc/intro.dox +++ b/examples/step-53/doc/intro.dox @@ -1,3 +1,368 @@ +
+ +This program was contributed by Wolfgang Bangerth and Luca Heltai, using +data provided by D. Sarah Stamps. + +@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. + +

Introduction

+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. + + +

Where geometry and meshes intersect

+ +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 chart. 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 atlas. + +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. + + +

The example case

+ +To illustrate how one describes geometries using charts in deal.II, we will +consider a case that originates in an application of the ASPECT mantle convection code, 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 East African Rift, a +zone where two continental plates drift apart. Not to beat around the bush, +the geometry we want to describe looks like this: + + + +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: + +
    +
  1. 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 WGS 84 + 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). + +
  2. 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. +
+ +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. + + +

Implementation

+ +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 push_forward() 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 pull_back() +implementing $\varphi$. We will do so in the AfricaGeometry 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 AfricaGeometry class will have +additional (non-virtual) member functions +AfricaGeometry::push_forward_wgs84() and +AfricaGeometry::push_forward_topo() 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 topography.txt.gz 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 latitude longitude elevation 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 ascending 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 value() 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 topography_data 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. + diff --git a/examples/step-53/doc/results.dox b/examples/step-53/doc/results.dox index b5eaba9377..feb821428c 100644 --- a/examples/step-53/doc/results.dox +++ b/examples/step-53/doc/results.dox @@ -1,2 +1,101 @@

Results

+Running the program produces a mesh file mesh.vtu that we can +visualize with any of the usual visualization programs that can read the VTU +file format. If one just looks at the mesh itself, it is actually very difficult +to see anything that doesn't just look like a perfectly round piece of a +sphere (though if one modified the program so that it does produce a sphere and +looked at them at the same time, the difference between the overall sphere and +WGS 84 shape is quite apparent). Apparently, Earth is actually quite a flat place. +Of course we already know this from satellite pictures. +However, we can tease out something more by +coloring cells by their volume. This both produces slight variations in hue +along the top surface and something for the visualization programs to apply +their shading algorithms to (because the top surfaces of the cells are now no +longer just tangential to a sphere but tilted): + + + +Yet, at least as far as visualizations are concerned, this is still not too +impressive. Rather, let us visualize things in a way so that we show the +actual elevation along the top surface. In other words, we want a picture like +this, with an incredible amount of detail: + + + +This image was produced with three small modifications: +
    +
  1. An addition seventh mesh refinement towards the top surface, + +
  2. The addition of the following piece function that, given a point + x computes the elevation by converting the point to + reference WGS 84 coordinates and only keeping the depth variable (the + function is, consequently, a simplified version of the + AfricaGeometry::pull_back_wgs84() version: + +@code +#include +#include +#include +#include + + +double get_elevation (const Point<3> &x) + { + const double R = 6378137; + const double ellipticity = 8.1819190842622e-2; + + 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 theta = std::atan2((x(2) + ep * ep * b * std::sin(th) * std::sin(th) * std::sin(th)), + (p - (ellipticity * ellipticity * R * (std::cos(th) * std::cos(th) * std::cos(th))))); + 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); + + return R_plus_d - R_bar; + } +@endcode + +
  3. Adding the following piece to the bottom of the run() function: + +@code + FE_Q<3> fe(1); + DoFHandler<3> dof_handler (triangulation); + dof_handler.distribute_dofs(fe); + + Vector elevation (dof_handler.n_dofs()); + { + std::map boundary_values; + VectorTools::interpolate_boundary_values(dof_handler, + 5, + ScalarFunctionFromFunctionObject<3>(get_elevation), + boundary_values); + for (std::map::const_iterator p = boundary_values.begin(); + p!=boundary_values.end(); ++p) + elevation[p->first] = p->second; + } + + DataOut<3> data_out; + data_out.attach_dof_handler(dof_handler); + data_out.add_data_vector (elevation, "elevation"); + data_out.build_patches(); + + std::ofstream out ("data.vtu"); + data_out.write_vtu (out); +@endcode +
+This last piece of code first creates a $Q_1$ finite element space on the mesh. +It then (ab)uses VectorTools::interpolate_boundary_values() to evaluate the +elevation function for every node at the top boundary (the one with boundary +indicator 5). We here wrap the call to get_elevation() with the +ScalarFunctionFromFunctionObject class to make a regular C++ function look +like an object of a class derived from the Function class that we want +to use in VectorTools::interpolate_boundary_values(). Having so gotten a list +of degrees of freedom located at the top boundary and corresponding elevation +values, we just go down this list and set these elevations in the +elevation vector (leaving all interior degrees of freedom at +their original zero value). This vector is then output using DataOut as +usual and can be visualized as shown above. + diff --git a/examples/step-53/doc/tooltip b/examples/step-53/doc/tooltip index 5e77ecb288..5cd4028333 100644 --- a/examples/step-53/doc/tooltip +++ b/examples/step-53/doc/tooltip @@ -1 +1 @@ -Geometry: Dealing with curved boundaries. +Geometry: Dealing with deformed domains. diff --git a/examples/step-53/step-53.cc b/examples/step-53/step-53.cc index 27ff0f0d8b..24d960ddea 100644 --- a/examples/step-53/step-53.cc +++ b/examples/step-53/step-53.cc @@ -15,216 +15,319 @@ * --------------------------------------------------------------------- * - * 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 (tria.h), and that allow +// us to create and output triangulations (grid_generator.h and +// grid_out.h). Furthermore, we need the header file that +// declares the Manifold and ChartManifold classes that we will need to +// describe the geometry (manifold.h). 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 -#include -#include #include #include +#include +#include -#include +// 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 -#include - - -using namespace dealii; - - -template -class SphereGeometry : public Boundary -{ -public: - SphereGeometry (const Point ¢er); - virtual - Point - get_new_point_on_line (const typename Triangulation::line_iterator &line) const; - - virtual - Point - get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const; -private: +#include +#include +#include - std_cxx11::array pull_back (const Point &p) const; - Point push_forward (const std_cxx11::array &preimage) const; - - template - static std_cxx11::array average (const std_cxx11::array (&array)[N]); - - const Point center; -}; - - -template -SphereGeometry::SphereGeometry (const Point ¢er) - : - center (center) -{} +#include -template <> -std_cxx11::array -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 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 topography_data + // 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,2> get_endpoints (); + static std_cxx1x::array n_intervals (); + static std::vector get_data (); + }; + + + // Let us move to the implementation of the class. The interesting parts + // of the class are the constructor and the value() 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 get_data() 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,2> + AfricaTopography::get_endpoints () + { + std_cxx1x::array,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 &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 + AfricaTopography::n_intervals () + { + std_cxx1x::array endpoints; + endpoints[0] = 379; + endpoints[1] = 219; + return endpoints; + } + // The only other function of greater interest is the get_data() + // function. It returns a temporary vector that contains all 751,740 data + // points describing the altitude and is read from the file + // topography.txt.gz. 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 + // topography_data 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 catch + // clause, we let the program run into an AssertThrow(false, ...) + // statement. Since the condition is always false, this always triggers an + // exception. In other words, this is equivalent to writing + // throw ExcMessage("...") 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 + AfricaTopography::get_data () + { + std::vector 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 !")); + } + } + + return data; + } -template <> -std_cxx11::array -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 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 &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 -std_cxx11::array -SphereGeometry<2>::average (const std_cxx11::array (&array)[N]) -{ - std_cxx11::array 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 + AfricaGeometry::push_forward(const Point<3> &chart_point) const { - bool origin_is_one_point = false; - - result[1] = 0; - for (unsigned int i=0; i 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 -std_cxx11::array -SphereGeometry<3>::average (const std_cxx11::array (&array)[N]) -{ - std_cxx11::array 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 ((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 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; } @@ -269,42 +372,151 @@ void make_grid () for (unsigned int i=0; i>>>>>> Replace the original attempt at providing step-53 by something real. - SphereGeometry geometry(center); - Triangulation 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, + // AfricaGeometry::push_forward() requires two arguments: + // the AfricaGeometry object to work with via its implicit + // this 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::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::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 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 true 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::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 run() 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; + } } + diff --git a/examples/step-53/topography.txt.gz b/examples/step-53/topography.txt.gz new file mode 100644 index 0000000000..223a85b7de Binary files /dev/null and b/examples/step-53/topography.txt.gz differ