]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Fix TransfiniteInterpolationManifold on complex geometry. 5397/head
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Sat, 4 Nov 2017 14:21:02 +0000 (15:21 +0100)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Sat, 4 Nov 2017 17:31:05 +0000 (18:31 +0100)
include/deal.II/grid/manifold_lib.h
source/grid/manifold_lib.cc
tests/manifold/transfinite_manifold_05.cc [new file with mode: 0644]
tests/manifold/transfinite_manifold_05.output [new file with mode: 0644]

index bb808be3dea2aefb50647195b523cd31087871df..82087a89e150fad7d3b6d68e9fa94a6851d1ad0e 100644 (file)
@@ -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 <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;
 
   /**
index 5256a5ead9da32ae197dd5ec8c2ac01d55e806d1..b493f0e68f6358af1e907f6c05ee9ddaeae3e092 100644 (file)
@@ -1224,7 +1224,7 @@ TransfiniteInterpolationManifold<dim,spacedim>
 
 
 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
 {
@@ -1289,9 +1289,9 @@ TransfiniteInterpolationManifold<dim,spacedim>
   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;
@@ -1309,14 +1309,12 @@ TransfiniteInterpolationManifold<dim, spacedim>
          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]);
@@ -1338,11 +1336,43 @@ TransfiniteInterpolationManifold<dim, spacedim>
         {
           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();
 }
 
diff --git a/tests/manifold/transfinite_manifold_05.cc b/tests/manifold/transfinite_manifold_05.cc
new file mode 100644 (file)
index 0000000..1a0f992
--- /dev/null
@@ -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 <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;
+}
diff --git a/tests/manifold/transfinite_manifold_05.output b/tests/manifold/transfinite_manifold_05.output
new file mode 100644 (file)
index 0000000..e21e849
--- /dev/null
@@ -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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.