From: Martin Kronbichler Date: Sat, 4 Nov 2017 14:21:02 +0000 (+0100) Subject: Fix TransfiniteInterpolationManifold on complex geometry. X-Git-Tag: v9.0.0-rc1~819^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=def568144b8433a4864a6572ef6a6412fd8c0381;p=dealii.git Fix TransfiniteInterpolationManifold on complex geometry. --- diff --git a/include/deal.II/grid/manifold_lib.h b/include/deal.II/grid/manifold_lib.h index bb808be3de..82087a89e1 100644 --- a/include/deal.II/grid/manifold_lib.h +++ b/include/deal.II/grid/manifold_lib.h @@ -731,15 +731,16 @@ public: private: /** * Internal function to identify the most suitable cells (=charts) where the - * given surrounding points are located. We use a cheap algorithm to identify - * the cells and rank the cells by probability before we actually do the - * search inside the relevant cells. The cells are sorted by the distance of - * a Q1 approximation of the inverse mapping to the unit cell of the - * surrounding points. We expect at most 10 cells (it should be less than 8 - * candidates even in 3D, typically only two or three), so get an array with - * 10 entries of a the indices cell->index(). - */ - std::array + * given surrounding points are located. We use a cheap algorithm to + * identify the cells and rank the cells by probability before we actually + * do the search inside the relevant cells. The cells are sorted by the + * distance of a Q1 approximation of the inverse mapping to the unit cell of + * the surrounding points. We expect at most 20 cells (it should be up to 8 + * candidates on a 3D structured mesh and a bit more on unstructured ones, + * typically we only get two or three), so get an array with 20 entries of a + * the indices cell->index(). + */ + std::array get_possible_cells_around_points(const ArrayView> &surrounding_points) const; /** diff --git a/source/grid/manifold_lib.cc b/source/grid/manifold_lib.cc index 5256a5ead9..b493f0e68f 100644 --- a/source/grid/manifold_lib.cc +++ b/source/grid/manifold_lib.cc @@ -1224,7 +1224,7 @@ TransfiniteInterpolationManifold template -std::array +std::array TransfiniteInterpolationManifold ::get_possible_cells_around_points(const ArrayView> &points) const { @@ -1289,9 +1289,9 @@ TransfiniteInterpolationManifold AssertThrow(distances_and_cells.size() > 0, (typename Mapping::ExcTransformationFailed())); std::sort(distances_and_cells.begin(), distances_and_cells.end()); - std::array cells; + std::array cells; cells.fill(numbers::invalid_unsigned_int); - for (unsigned int i=0; i ExcMessage("The chart points array view must be as large as the " "surrounding points array view.")); - std::array nearby_cells = + std::array nearby_cells = get_possible_cells_around_points(surrounding_points); // check whether all points are inside the unit cell of the current chart for (unsigned int c=0; c::ExcTransformationFailed())); typename Triangulation::cell_iterator cell(triangulation, level_coarse, nearby_cells[c]); @@ -1338,11 +1336,43 @@ TransfiniteInterpolationManifold { return cell; } + + // if we did not find a point and this was the last valid cell (the next + // iterate being the end of the array or an unvalid tag), we must stop + if (c == nearby_cells.size()-1 || + nearby_cells[c+1] == numbers::invalid_unsigned_int) + { + // generate additional information to help debugging why we did not + // get a point + std::ostringstream message; + for (unsigned int b=0; b::cell_iterator cell(triangulation, + level_coarse, + nearby_cells[b]); + message << "Looking at cell " + << cell->id() << " with vertices: " << std::endl; + for (unsigned int v=0; v::vertices_per_cell; ++v) + message << cell->vertex(v) << " "; + message << std::endl; + message << "Transformation to chart coordinates: " << std::endl; + for (unsigned int i=0; i " + << pull_back(cell, surrounding_points[i]) + << std::endl; + } + } + + AssertThrow(false, + (typename Mapping::ExcTransformationFailed(message.str()))); + } } - // a valid inversion should have returned a point above. - AssertThrow(false, - (typename Mapping::ExcTransformationFailed())); + // a valid inversion should have returned a point above. an invalid + // inversion should have triggered the assertion, so we should never end up + // here + Assert(false, ExcInternalError()); return typename Triangulation::cell_iterator(); } diff --git a/tests/manifold/transfinite_manifold_05.cc b/tests/manifold/transfinite_manifold_05.cc new file mode 100644 index 0000000000..1a0f992c62 --- /dev/null +++ b/tests/manifold/transfinite_manifold_05.cc @@ -0,0 +1,65 @@ +//------------------------------------------------------------------- +// Copyright (C) 2017 by the deal.II authors. +// +// This file is subject to LGPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//------------------------------------------------------------------- + + +// Test that transfinite interpolation manifold works properly for creating a +// particular point on a somewhat more complicated geometry. We used to +// restrict the search too much in an initial version of the manifold + +#include "../tests.h" +#include +#include +#include +#include +#include +#include + + +int main () +{ + initlog(); + + const int dim = 3; + Triangulation tria1, tria2, tria; + GridGenerator::hyper_shell(tria1, Point(), 0.4, std::sqrt(dim), 6); + GridGenerator::hyper_ball(tria2, Point(), 0.4); + GridGenerator::merge_triangulations(tria1, tria2, tria); + tria.set_all_manifold_ids(0); + for (typename Triangulation::cell_iterator cell = tria.begin(); + cell != tria.end(); ++cell) + { + for (unsigned int f=0; f::faces_per_cell; ++f) + { + bool face_at_sphere_boundary = true; + for (unsigned int v=0; v::vertices_per_cell; ++v) + if (std::abs(cell->face(f)->vertex(v).norm()-0.4) > 1e-12) + face_at_sphere_boundary = false; + if (face_at_sphere_boundary) + cell->face(f)->set_all_manifold_ids(1); + } + } + static const SphericalManifold spherical_manifold; + tria.set_manifold(1, spherical_manifold); + static TransfiniteInterpolationManifold transfinite; + transfinite.initialize(tria); + tria.set_manifold(0, transfinite); + + const std::array, 2> points({{Point<3>(0, 0.360566, 0), Point<3>(0, 0.321132, 0)}}); + const std::array weights({{0.4, 0.6}}); + deallog << "Interpolate between points " << points[0] << " and " + << points[1] << " with weights " << weights[0] << " and " + << weights[1] << ": " << transfinite.get_new_point(make_array_view(points.begin(), + points.end()), + make_array_view(weights.begin(), + weights.end())) + << std::endl; + + return 0; +} diff --git a/tests/manifold/transfinite_manifold_05.output b/tests/manifold/transfinite_manifold_05.output new file mode 100644 index 0000000000..e21e84989e --- /dev/null +++ b/tests/manifold/transfinite_manifold_05.output @@ -0,0 +1,2 @@ + +DEAL::Interpolate between points 0.00000 0.360566 0.00000 and 0.00000 0.321132 0.00000 with weights 0.400000 and 0.600000: 0.00000 0.336906 0.00000