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 <tt>cell->index()</tt>.
- */
- std::array<unsigned int, 10>
+ * 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 <tt>cell->index()</tt>.
+ */
+ std::array<unsigned int, 20>
get_possible_cells_around_points(const ArrayView<const Point<spacedim>> &surrounding_points) const;
/**
template <int dim, int spacedim>
-std::array<unsigned int, 10>
+std::array<unsigned int, 20>
TransfiniteInterpolationManifold<dim,spacedim>
::get_possible_cells_around_points(const ArrayView<const Point<spacedim>> &points) const
{
AssertThrow(distances_and_cells.size() > 0,
(typename Mapping<dim,spacedim>::ExcTransformationFailed()));
std::sort(distances_and_cells.begin(), distances_and_cells.end());
- std::array<unsigned int,10> cells;
+ std::array<unsigned int,20> cells;
cells.fill(numbers::invalid_unsigned_int);
- for (unsigned int i=0; i<distances_and_cells.size() && i<10; ++i)
+ for (unsigned int i=0; i<distances_and_cells.size() && i<cells.size(); ++i)
cells[i] = distances_and_cells[i].second;
return cells;
ExcMessage("The chart points array view must be as large as the "
"surrounding points array view."));
- std::array<unsigned int,10> nearby_cells =
+ std::array<unsigned int,20> 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<nearby_cells.size(); ++c)
{
- AssertThrow(nearby_cells[c] != numbers::invalid_unsigned_int,
- (typename Mapping<dim,spacedim>::ExcTransformationFailed()));
typename Triangulation<dim,spacedim>::cell_iterator cell(triangulation,
level_coarse,
nearby_cells[c]);
{
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<c; ++b)
+ {
+ typename Triangulation<dim,spacedim>::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<GeometryInfo<dim>::vertices_per_cell; ++v)
+ message << cell->vertex(v) << " ";
+ message << std::endl;
+ message << "Transformation to chart coordinates: " << std::endl;
+ for (unsigned int i=0; i<surrounding_points.size(); ++i)
+ {
+ message << surrounding_points[i] << " -> "
+ << pull_back(cell, surrounding_points[i])
+ << std::endl;
+ }
+ }
+
+ AssertThrow(false,
+ (typename Mapping<dim,spacedim>::ExcTransformationFailed(message.str())));
+ }
}
- // a valid inversion should have returned a point above.
- AssertThrow(false,
- (typename Mapping<dim,spacedim>::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<dim,spacedim>::cell_iterator();
}
--- /dev/null
+//-------------------------------------------------------------------
+// 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 <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/manifold_lib.h>
+
+
+int main ()
+{
+ initlog();
+
+ const int dim = 3;
+ Triangulation<dim> tria1, tria2, tria;
+ GridGenerator::hyper_shell(tria1, Point<dim>(), 0.4, std::sqrt(dim), 6);
+ GridGenerator::hyper_ball(tria2, Point<dim>(), 0.4);
+ GridGenerator::merge_triangulations(tria1, tria2, tria);
+ tria.set_all_manifold_ids(0);
+ for (typename Triangulation<dim>::cell_iterator cell = tria.begin();
+ cell != tria.end(); ++cell)
+ {
+ for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+ {
+ bool face_at_sphere_boundary = true;
+ for (unsigned int v=0; v<GeometryInfo<dim-1>::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<dim> spherical_manifold;
+ tria.set_manifold(1, spherical_manifold);
+ static TransfiniteInterpolationManifold<dim> transfinite;
+ transfinite.initialize(tria);
+ tria.set_manifold(0, transfinite);
+
+ const std::array<Point<3>, 2> points({{Point<3>(0, 0.360566, 0), Point<3>(0, 0.321132, 0)}});
+ const std::array<double, 2> 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;
+}
--- /dev/null
+
+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