From 34cf8d3cfa6e1ee5772ffedc0ce894edc1a82f96 Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Wed, 8 May 2019 23:01:58 -0400 Subject: [PATCH] Modernize step-53 --- examples/step-53/step-53.cc | 73 +++++++++++-------------------------- 1 file changed, 21 insertions(+), 52 deletions(-) diff --git a/examples/step-53/step-53.cc b/examples/step-53/step-53.cc index 3e4eb436a5..602d8f6f0e 100644 --- a/examples/step-53/step-53.cc +++ b/examples/step-53/step-53.cc @@ -88,9 +88,7 @@ namespace Step53 private: const Functions::InterpolatedUniformGridData<2> topography_data; - static std::array, 2> get_endpoints(); - static std::array n_intervals(); - static std::vector get_data(); + static std::vector get_data(); }; @@ -112,8 +110,9 @@ namespace Step53 // 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())) {} @@ -125,24 +124,6 @@ namespace Step53 } - std::array, 2> AfricaTopography::get_endpoints() - { - std::array, 2> endpoints; - endpoints[0] = std::make_pair(-6.983333, 11.966667); - endpoints[1] = std::make_pair(25, 35.95); - return endpoints; - } - - - std::array AfricaTopography::n_intervals() - { - std::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 83,600 data // points describing the altitude and is read from the file @@ -317,20 +298,18 @@ namespace Step53 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}; } @@ -354,13 +333,8 @@ namespace Step53 // 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. + // 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; @@ -377,10 +351,11 @@ namespace Step53 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 @@ -393,10 +368,7 @@ namespace Step53 // 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 @@ -412,10 +384,7 @@ namespace Step53 // 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) { -- 2.39.5