private:
const Functions::InterpolatedUniformGridData<2> topography_data;
- static std::array<std::pair<double, double>, 2> get_endpoints();
- static std::array<unsigned int, 2> n_intervals();
- static std::vector<double> get_data();
+ static std::vector<double> get_data();
};
// 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(),
+ : topography_data({std::make_pair(-6.983333, 11.966667),
+ std::make_pair(25, 35.95)},
+ {379, 219},
Table<2, double>(380, 220, get_data().begin()))
{}
}
- std::array<std::pair<double, double>, 2> AfricaTopography::get_endpoints()
- {
- std::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;
- }
-
-
- std::array<unsigned int, 2> AfricaTopography::n_intervals()
- {
- std::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 83,600 data
// points describing the altitude and is read from the file
Point<3>
AfricaGeometry::push_forward_topo(const Point<3> &phi_theta_d_hat) const
{
- 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;
+ 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;
+ return {phi_theta_d_hat[0], phi_theta_d_hat[1], 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;
+ 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);
+ return {phi_theta_d[0], phi_theta_d[1], d_hat};
}
// 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.
+ // function of the geometry we use. We wrap it by a lambda function to
+ // obtaining the desired object required for the transformation.
void run()
{
AfricaGeometry geometry;
GridGenerator::subdivided_hyper_rectangle(
triangulation, subdivisions, corner_points[0], corner_points[1], true);
- GridTools::transform(std::bind(&AfricaGeometry::push_forward,
- std::cref(geometry),
- std::placeholders::_1),
- triangulation);
+ GridTools::transform(
+ [&geometry](const Point<3> &chart_point) {
+ return geometry.push_forward(chart_point);
+ },
+ triangulation);
}
// The next step is to explain to the triangulation to use our geometry
// 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)
+ for (const auto &cell : triangulation.active_cell_iterators())
cell->set_all_manifold_ids(0);
// The last step is to refine the mesh beyond its initial $1\times 2\times
// 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 (const auto &cell : triangulation.active_cell_iterators())
for (unsigned int f = 0; f < GeometryInfo<3>::faces_per_cell; ++f)
if (cell->face(f)->boundary_id() == 5)
{