From 3af08eab4d9deaee91633fe9425babc2c4661ad8 Mon Sep 17 00:00:00 2001 From: Matthias Maier Date: Mon, 2 Nov 2020 17:16:11 -0600 Subject: [PATCH] Transfinite interpolation: add a test --- tests/manifold/transfinite_manifold_11.cc | 155 ++++++++++++++++++ tests/manifold/transfinite_manifold_11.output | 1 + 2 files changed, 156 insertions(+) create mode 100644 tests/manifold/transfinite_manifold_11.cc create mode 100644 tests/manifold/transfinite_manifold_11.output diff --git a/tests/manifold/transfinite_manifold_11.cc b/tests/manifold/transfinite_manifold_11.cc new file mode 100644 index 0000000000..e7dba1dcfc --- /dev/null +++ b/tests/manifold/transfinite_manifold_11.cc @@ -0,0 +1,155 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 - 2020 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +// This test verifies that the transfinite interpolation works on a torus + +#include +#include + +#include + +#include "../tests.h" + +using namespace dealii; + +template +class GradingManifold : public ChartManifold +{ +public: + GradingManifold(const Point center, + const double grading, + const double epsilon) + : center(center) + , grading(grading) + , epsilon(epsilon) + , polar_manifold(center) + {} + + virtual Point + pull_back(const Point &space_point) const final override + { + auto point = polar_manifold.pull_back(space_point); + Assert(point[0] >= 0., ExcInternalError()); + point[0] = std::pow(point[0] + epsilon, 1. / grading) - + std::pow(epsilon, 1. / grading) + 1.e-14; + const auto chart_point = polar_manifold.push_forward(point); + return chart_point; + } + + virtual Point + push_forward(const Point &chart_point) const final override + { + auto point = polar_manifold.pull_back(chart_point); + point[0] = std::pow(point[0] + std::pow(epsilon, 1. / grading), grading) - + epsilon + 1.e-14; + Assert(point[0] >= 0., ExcInternalError()); + return polar_manifold.push_forward(point); + } + + std::unique_ptr> + clone() const final override + { + return std::make_unique>(center, grading, epsilon); + } + +private: + const Point center; + const double grading; + const double epsilon; + + PolarManifold polar_manifold; +}; + + +int +main() +{ + initlog(); + deallog << std::setprecision(9); + + const std::vector> vertices{ + {0.5 * 2.5, -std::sqrt(3.) / 2. * 2.5}, // 0 + {-0.9 + 0.7, -0.025}, // 1 + {-0.9 + 0.7, 0.025}, // 2 + {0.5 * 2.5, std::sqrt(3.) / 2. * 2.5}, // 3 + {2.5, -0.5 * 2.5}, // 4 + {2.5, -0.416667}, // 5 + {2.5, 0.416667}, // 6 + {2.5, 0.5 * 2.5}, // 7 + }; + + std::vector> cells(3); + cells[0].vertices = {0, 4, 1, 5}; + cells[1].vertices = {1, 5, 2, 6}; + cells[2].vertices = {2, 6, 3, 7}; + + Triangulation<2> triangulation; + triangulation.create_triangulation(vertices, cells, SubCellData()); + + { + triangulation.set_all_manifold_ids(3); + const auto center_cell = ++triangulation.begin_active(); + const auto lower_radial = center_cell->face(2); + const auto upper_radial = center_cell->face(3); + lower_radial->set_manifold_id(1); + upper_radial->set_manifold_id(2); + } + + GradingManifold<2> lower_radial{Point<2>{-0.9 + 0.7, -0.025}, 3.0, 0.01}; + triangulation.set_manifold(1, lower_radial); + + GradingManifold<2> upper_radial{Point<2>{-0.9 + 0.7, 0.025}, 3.0, 0.01}; + triangulation.set_manifold(2, upper_radial); + + TransfiniteInterpolationManifold<2> transfinite; + transfinite.initialize(triangulation); + triangulation.set_manifold(3, transfinite); + + /* + * triangulation.refine_global(8) fails with: + * + * An error occurred in line [...] in function + * [...]compute_chart_points(const ArrayView> &, + * ArrayView>) The violated condition was: false Additional + * information: Looking at cell 1_0: with vertices: -0.2 -0.025 2.5 + * -0.416667 -0.2 0.025 2.5 0.416667 Transformation to chart + * coordinates: -0.1271 -0.0177875 -> 0.1875 0.25 -0.1271 -0.015564 -> 20 20 + * Looking at cell 0_0: with vertices: + * 1.25 -2.16506 2.5 -1.25 -0.2 -0.025 2.5 -0.416667 + * Transformation to chart coordinates: + * -0.1271 -0.0177875 -> 0.211251 1.01047 + * -0.1271 -0.015564 -> 0.214091 1.01181 + * Looking at cell 2_0: with vertices: + * -0.2 0.025 2.5 0.416667 1.25 2.16506 2.5 1.25 + * Transformation to chart coordinates: + * -0.1271 -0.0177875 -> 0.253794 -0.0323121 + * -0.1271 -0.015564 -> 0.251301 -0.0309146 + * + * The following is a simplified version that fails with a sligthly + * different error message. + */ + + std::vector> points{Point<2>{-0.1271, -0.0177875}, + Point<2>{-0.1271, -0.015564}}; + ArrayView> points_view{points.data(), 2}; + + std::vector weights{0.5, 0.5}; + ArrayView weights_view{weights.data(), 2}; + + transfinite.get_new_point(points_view, weights_view); + + return 0; +} diff --git a/tests/manifold/transfinite_manifold_11.output b/tests/manifold/transfinite_manifold_11.output new file mode 100644 index 0000000000..8b13789179 --- /dev/null +++ b/tests/manifold/transfinite_manifold_11.output @@ -0,0 +1 @@ + -- 2.39.5