diff --git a/include/aspect/geometry_model/box.h b/include/aspect/geometry_model/box.h index a8f7931a8c1..de695292bf5 100644 --- a/include/aspect/geometry_model/box.h +++ b/include/aspect/geometry_model/box.h @@ -223,6 +223,11 @@ namespace aspect parse_parameters (ParameterHandler &prm) override; private: + /** + * A pointer to the initial topography model. + */ + InitialTopographyModel::Interface *topo_model; + /** * Extent of the box in x-, y-, and z-direction (in 3d). */ @@ -242,11 +247,6 @@ namespace aspect * The number of cells in each coordinate direction. */ std::array repetitions; - - /** - * A pointer to the initial topography model. - */ - InitialTopographyModel::Interface *topo_model; }; } } diff --git a/include/aspect/geometry_model/chunk.h b/include/aspect/geometry_model/chunk.h index b4f6ac21391..ed34f5801f3 100644 --- a/include/aspect/geometry_model/chunk.h +++ b/include/aspect/geometry_model/chunk.h @@ -131,6 +131,11 @@ namespace aspect clone() const override; private: + /** + * A pointer to the topography model. + */ + const InitialTopographyModel::Interface *topo; + /** * The minimum longitude of the domain. */ @@ -164,11 +169,6 @@ namespace aspect virtual Point push_forward_topo(const Point &chart_point) const; - - /** - * A pointer to the topography model. - */ - const InitialTopographyModel::Interface *topo; }; } diff --git a/include/aspect/geometry_model/ellipsoidal_chunk.h b/include/aspect/geometry_model/ellipsoidal_chunk.h index fda07203bf1..516162ddfad 100644 --- a/include/aspect/geometry_model/ellipsoidal_chunk.h +++ b/include/aspect/geometry_model/ellipsoidal_chunk.h @@ -46,29 +46,17 @@ namespace aspect /** * Constructor */ - EllipsoidalChunkGeometry(); + EllipsoidalChunkGeometry(const InitialTopographyModel::Interface &topography, + const double para_semi_major_axis_a, + const double para_eccentricity, + const double para_semi_minor_axis_b, + const double para_bottom_depth, + const std::vector> ¶_corners); /** * Copy constructor */ - EllipsoidalChunkGeometry(const EllipsoidalChunkGeometry &other); - - /** - * An initialization function necessary to make sure that the - * manifold has access to the topography plugins. - */ - void - initialize(const InitialTopographyModel::Interface *topography); - - /** - * Sets several parameters for the ellipsoidal manifold object. - */ - void - set_manifold_parameters(const double para_semi_major_axis_a, - const double para_eccentricity, - const double para_semi_minor_axis_b, - const double para_bottom_depth, - const std::vector> ¶_corners); + EllipsoidalChunkGeometry(const EllipsoidalChunkGeometry &other) = default; /** * The deal.ii pull back function in 3d. This function receives @@ -127,13 +115,16 @@ namespace aspect */ Point<3> pull_back_topography (const Point<3> &phi_theta_d) const; - - double semi_major_axis_a; - double eccentricity; - double semi_minor_axis_b; - double bottom_depth; - std::vector> corners; + /** + * A pointer to the topography model. + */ const InitialTopographyModel::Interface *topography; + + const double semi_major_axis_a; + const double eccentricity; + const double semi_minor_axis_b; + const double bottom_depth; + const std::vector> corners; }; } @@ -333,9 +324,17 @@ namespace aspect /** - * Construct manifold object Pointer to an object that describes the geometry. + * An object that describes the geometry. This pointer is + * initialized in the initialize() function, and serves as the manifold + * object that the triangulation is later given in create_coarse_mesh() + * where the triangulation clones it. + * + * The object is marked as 'const' to make it clear that it should not + * be modified once created. That is because the triangulation copies it, + * and modifying the current object will not have any impact on the + * manifold used by the triangulation. */ - internal::EllipsoidalChunkGeometry manifold; + std::unique_ptr> manifold; void set_boundary_ids(parallel::distributed::Triangulation &coarse_grid) const; diff --git a/source/geometry_model/box.cc b/source/geometry_model/box.cc index 4d796ee6b11..882efeb7e68 100644 --- a/source/geometry_model/box.cc +++ b/source/geometry_model/box.cc @@ -41,6 +41,7 @@ namespace aspect { // Get pointer to initial topography model topo_model = const_cast*>(&this->get_initial_topography_model()); + // Check that initial topography is required. // If so, connect the initial topography function // to the right signals: It should be applied after diff --git a/source/geometry_model/chunk.cc b/source/geometry_model/chunk.cc index 9b264e96843..60faf061aac 100644 --- a/source/geometry_model/chunk.cc +++ b/source/geometry_model/chunk.cc @@ -43,10 +43,10 @@ namespace aspect const double min_radius, const double max_depth) : + topo (&topo), point1_lon(min_longitude), inner_radius(min_radius), - max_depth(max_depth), - topo (&topo) + max_depth(max_depth) {} diff --git a/source/geometry_model/ellipsoidal_chunk.cc b/source/geometry_model/ellipsoidal_chunk.cc index 0b4e66999a4..94200b97708 100644 --- a/source/geometry_model/ellipsoidal_chunk.cc +++ b/source/geometry_model/ellipsoidal_chunk.cc @@ -66,43 +66,24 @@ namespace aspect // Constructor template - EllipsoidalChunkGeometry::EllipsoidalChunkGeometry() + EllipsoidalChunkGeometry::EllipsoidalChunkGeometry(const InitialTopographyModel::Interface &topo, + const double para_semi_major_axis_a, + const double para_eccentricity, + const double para_semi_minor_axis_b, + const double para_bottom_depth, + const std::vector> ¶_corners) : - semi_major_axis_a (-1), - eccentricity (-1), - semi_minor_axis_b (-1), - bottom_depth (-1), - topography(nullptr) - {} - - // Copy constructor - template - EllipsoidalChunkGeometry::EllipsoidalChunkGeometry(const EllipsoidalChunkGeometry &other) - : - ChartManifold(other) + topography (&topo), + semi_major_axis_a (para_semi_major_axis_a), + eccentricity (para_eccentricity), + semi_minor_axis_b (para_semi_minor_axis_b), + bottom_depth (para_bottom_depth), + corners (para_corners) { - this->initialize(other.topography); - this->set_manifold_parameters(other.semi_major_axis_a, other.eccentricity, other.semi_minor_axis_b, - other.bottom_depth, other.corners); + AssertThrow (dim == 3, ExcMessage("This manifold can currently only be used in 3d.")); } - template - void - EllipsoidalChunkGeometry::set_manifold_parameters(const double para_semi_major_axis_a, - const double para_eccentricity, - const double para_semi_minor_axis_b, - const double para_bottom_depth, - const std::vector> ¶_corners) - { - AssertThrow (dim == 3, ExcMessage("This manifold can currently only be used in 3d.")); - - semi_major_axis_a = para_semi_major_axis_a; - eccentricity = para_eccentricity; - semi_minor_axis_b = para_semi_minor_axis_b; - bottom_depth = para_bottom_depth; - corners=para_corners; - } template Point<3> @@ -225,23 +206,18 @@ namespace aspect { return std::make_unique(*this); } - - - - template - void - EllipsoidalChunkGeometry::initialize(const InitialTopographyModel::Interface *topography_) - { - AssertThrow (dim == 3, ExcMessage("This manifold can currently only be used in 3d.")); - topography = topography_; - } } template void EllipsoidalChunk::initialize() { - manifold.initialize(&(this->get_initial_topography_model())); + manifold = std::make_unique>(this->get_initial_topography_model(), + semi_major_axis_a, + eccentricity, + semi_minor_axis_b, + bottom_depth, + corners); } @@ -280,7 +256,7 @@ namespace aspect GridTools::transform ( [&](const Point &x) -> Point { - return manifold.push_forward(x); + return manifold->push_forward(x); }, coarse_grid); @@ -288,7 +264,7 @@ namespace aspect // We won't use it during regular operation, but we set manifold_ids for // all cells, faces and edges immediately before refinement and // clear it again afterwards - coarse_grid.set_manifold (15, manifold); + coarse_grid.set_manifold (15, *manifold); coarse_grid.signals.pre_refinement.connect ( [&] {set_manifold_ids(coarse_grid);}); @@ -627,21 +603,13 @@ namespace aspect prm.leave_subsection(); } prm.leave_subsection(); - - - // Construct manifold object Pointer to an object that describes the geometry. - manifold.set_manifold_parameters(semi_major_axis_a, - eccentricity, - semi_minor_axis_b, - bottom_depth, - corners); } template double EllipsoidalChunk::depth(const Point &position) const { - return std::max(std::min(-manifold.pull_back(position)[dim-1], maximal_depth()), 0.0); + return std::max(std::min(-manifold->pull_back(position)[dim-1], maximal_depth()), 0.0); } template @@ -666,7 +634,7 @@ namespace aspect double EllipsoidalChunk::get_radius(const Point &position) const { - const Point long_lat_depth = manifold.pull_back(position); + const Point long_lat_depth = manifold->pull_back(position); return semi_major_axis_a / (std::sqrt(1 - eccentricity * eccentricity * std::sin(long_lat_depth[1]) * std::sin(long_lat_depth[1]))); } @@ -732,7 +700,7 @@ namespace aspect (southLatitude + northLatitude) * 0.5 * constants::degree_to_radians, -depth); - return manifold.push_forward(p); + return manifold->push_forward(p); } @@ -745,7 +713,7 @@ namespace aspect ExcMessage("After displacement of the mesh, this function can no longer be used to determine whether a point lies in the domain or not.")); // dim = 3 - const Point ellipsoidal_point = manifold.pull_back(point); + const Point ellipsoidal_point = manifold->pull_back(point); // compare deflection from the ellipsoid surface if (ellipsoidal_point[dim-1] > 0.0+std::numeric_limits::epsilon()*bottom_depth || @@ -771,7 +739,7 @@ namespace aspect for (unsigned int d=0; d transformed_point = manifold.pull_back_ellipsoid(cartesian_point, semi_major_axis_a, eccentricity); + Point<3> transformed_point = manifold->pull_back_ellipsoid(cartesian_point, semi_major_axis_a, eccentricity); const double radius = semi_major_axis_a / (std::sqrt(1 - eccentricity * eccentricity * std::sin(transformed_point[1]) * std::sin(transformed_point[1]))); @@ -807,7 +775,7 @@ namespace aspect const double radius = semi_major_axis_a / (std::sqrt(1 - eccentricity * eccentricity * std::sin(position_point(1)) * std::sin(position_point(1)))); position_point(2) = position_tensor[0] - radius; - Point<3> transformed_point = manifold.push_forward_ellipsoid(position_point, semi_major_axis_a, eccentricity); + Point<3> transformed_point = manifold->push_forward_ellipsoid(position_point, semi_major_axis_a, eccentricity); return transformed_point; } @@ -827,7 +795,7 @@ template typename aspect::GeometryModel::internal::EllipsoidalChunkGeometry aspect::GeometryModel::EllipsoidalChunk::get_manifold() const { - return manifold; + return *manifold; } // explicit instantiations diff --git a/tests/ellipsoidal_chunk_geometry.cc b/tests/ellipsoidal_chunk_geometry.cc index 93a18215873..fa6838e8c1e 100644 --- a/tests/ellipsoidal_chunk_geometry.cc +++ b/tests/ellipsoidal_chunk_geometry.cc @@ -53,20 +53,9 @@ int f() const unsigned int dim=3; - InitialTopographyModel::ZeroTopography topography; - GeometryModel::internal::EllipsoidalChunkGeometry ellipsoidal_manifold; - ellipsoidal_manifold.initialize(&topography); - std::vector> corners(2,Point<2>(-15.0,-15.0)); corners[1] *= -1.0; - std::cout << "Simple sphere test" << std::endl; - ellipsoidal_manifold.set_manifold_parameters(6371000.0, - 0.0, - 6371000.0, - 2890000.0, - corners); - std::vector> test_points; test_points.push_back(Point<3> (6371000.0,0,0)); test_points.push_back(Point<3> (6171000.0,0,0)); @@ -78,33 +67,50 @@ int f() test_points.push_back(Point<3> (25000.0,25000.0,500000.0)); test_points.push_back(Point<3> (25000.0,25000.0,50000.0)); - for (unsigned int i=0; i topography; + { + std::cout << "Simple sphere test" << std::endl; + GeometryModel::internal::EllipsoidalChunkGeometry ellipsoidal_manifold(topography, + 6371000.0, + 0.0, + 6371000.0, + 2890000.0, + corners); + + for (unsigned int i=0; i ellipsoidal_manifold(topography, + semi_major_axis_a, + eccentricity, + semi_minor, + 2890000.0, + corners); + + for (unsigned int i=0; i ellipsoidal_manifold(topography, + semi_major_axis_a, + eccentricity, + semi_minor, + 500000.0, + corners); + + for (unsigned int i=0; i