]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Raise exceptions for invalid input.
authorDavid Wells <wellsd2@rpi.edu>
Mon, 9 May 2016 01:20:14 +0000 (21:20 -0400)
committerDavid Wells <wellsd2@rpi.edu>
Sat, 14 May 2016 19:39:40 +0000 (15:39 -0400)
subdivided_parallelepiped could generate cells with negative volume if
the provided edges were not ordered correctly. Try to help the user by
raising an exception if this is the case in subdivided_parallelepiped
instead of create_triangulation.

This commit also fixes tests that would fail as a result of removing the
orientation fix.

Partly reverts ccb56e7fb5.

include/deal.II/grid/grid_generator.h
source/grid/grid_generator.cc
tests/grid/grid_parallelepiped.cc
tests/grid/grid_parallelepiped_01.cc
tests/grid/grid_parallelepiped_02.cc
tests/grid/grid_parallelepiped_03.cc
tests/grid/grid_parallelepiped_05.cc
tests/grid/grid_parallelepiped_05.output

index c70c33ca42cfa50b8cc6ea2d0f13b21b5c35a0ca..cace929cdaf308e4d123dab0007864a28f336b2d 100644 (file)
@@ -1139,6 +1139,12 @@ namespace GridGenerator
                   << "The vector of repetitions  must have "
                   << arg1 <<" elements.");
 
+  /**
+   * Exception for input that is not properly oriented.
+   */
+  DeclExceptionMsg (ExcInvalidInputOrientation,
+                    "The input to this function is oriented in a way that will"
+                    " cause all cells to have negative measure.");
   ///@}
 }
 
index 341b82f339f91f4af383279b4a1de17ab792220c..20354c17ee1049b8ae78ae09d9f94e9ab7d9d83a 100644 (file)
@@ -856,9 +856,6 @@ namespace GridGenerator
                              const std::vector<unsigned int> &subdivisions,
                              const bool colorize)
   {
-    // The provided edges and subdivisions may not be valid, so possibly
-    // recompute both of them.
-    std_cxx11::array<Tensor<1, spacedim>, dim> compute_edges = edges;
     std::vector<unsigned int> compute_subdivisions = subdivisions;
     if (compute_subdivisions.size() == 0)
       {
@@ -876,23 +873,75 @@ namespace GridGenerator
       }
 
     /*
-     * Verify that the vectors are oriented in a counter clockwise direction in 2D.
+     * Verify that the edge points to the right in 1D, vectors are oriented in
+     * a counter clockwise direction in 2D, or form a right handed system in
+     * 3D.
      */
-    if (dim == 2)
+    bool twisted_data = false;
+    switch (dim)
+      {
+      case 1:
+      {
+        twisted_data = (edges[0][0] < 0);
+        break;
+      }
+      case 2:
       {
-        const double plane_normal = edges[0][0]*edges[1][1] - edges[0][1]*edges[1][0];
-        if (plane_normal < 0.0)
+        if (spacedim == 2) // this check does not make sense otherwise
           {
-            std::swap(compute_edges[0], compute_edges[1]);
-            std::swap(compute_subdivisions[0], compute_subdivisions[1]);
+            const double plane_normal = edges[0][0]*edges[1][1] - edges[0][1]*edges[1][0];
+            twisted_data = (plane_normal < 0.0);
           }
+        break;
       }
-
+      case 3:
+      {
+        // Check that the first two vectors are not linear combinations to
+        // avoid zero division later on.
+        Assert(std::abs(edges[0]*edges[1]
+                        /(edges[0].norm()*edges[1].norm())
+                        - 1.0) > 1.0e-15,
+               ExcMessage("Edges in subdivided_parallelepiped() must point in"
+                          " different directions."));
+        const Tensor<1, spacedim> plane_normal = cross_product_3d
+                                                 (edges[0], edges[1]);
+
+        /*
+         * Ensure that edges 1, 2, and 3 form a right-handed set of
+         * vectors. This works by applying the definition of the dot product
+         *
+         *     cos(theta) = dot(x, y)/(norm(x)*norm(y))
+         *
+         * and then, since the normal vector and third edge should both point
+         * away from the plane formed by the first two edges, the angle
+         * between them must be between 0 and pi/2; hence we just need
+         *
+         *     0 < dot(x, y).
+         */
+        twisted_data = (plane_normal*edges[2] < 0.0);
+        break;
+      }
+      default:
+        Assert(false, ExcInternalError());
+      }
+    (void)twisted_data; // make the static analyzer happy
+    Assert(!twisted_data,
+           ExcInvalidInputOrientation
+           ("The triangulation you are trying to create will consist of cells"
+            " with negative measures. This is usually the result of input data"
+            " that does not define a right-handed coordinate system. The usual"
+            " fix for this is to ensure that in 1D the given point is to the"
+            " right of the origin (or the given edge tensor is positive), in 2D"
+            " that the two edges (and their cross product) obey the right-hand"
+            " rule (which may usually be done by switching the order of the"
+            " points or edge tensors), or in 3D that the edges form a"
+            " right-handed coordinate system (which may also be accomplished by"
+            " switching the order of the first two points or edge tensors)."));
 
     // Check corners do not overlap (unique)
     for (unsigned int i=0; i<dim; ++i)
       for (unsigned int j=i+1; j<dim; ++j)
-        Assert ((compute_edges[i]!=compute_edges[j]),
+        Assert ((edges[i]!=edges[j]),
                 ExcMessage ("Degenerate edges of subdivided_parallelepiped encountered."));
 
     // Create a list of points
@@ -902,15 +951,15 @@ namespace GridGenerator
       {
       case 1:
         for (unsigned int x=0; x<=compute_subdivisions[0]; ++x)
-          points.push_back (origin + compute_edges[0]/compute_subdivisions[0]*x);
+          points.push_back (origin + edges[0]/compute_subdivisions[0]*x);
         break;
 
       case 2:
         for (unsigned int y=0; y<=compute_subdivisions[1]; ++y)
           for (unsigned int x=0; x<=compute_subdivisions[0]; ++x)
             points.push_back (origin
-                              + compute_edges[0]/compute_subdivisions[0]*x
-                              + compute_edges[1]/compute_subdivisions[1]*y);
+                              + edges[0]/compute_subdivisions[0]*x
+                              + edges[1]/compute_subdivisions[1]*y);
         break;
 
       case 3:
@@ -919,9 +968,9 @@ namespace GridGenerator
             for (unsigned int x=0; x<=compute_subdivisions[0]; ++x)
               points.push_back (
                 origin
-                + compute_edges[0]/compute_subdivisions[0]*x
-                + compute_edges[1]/compute_subdivisions[1]*y
-                + compute_edges[2]/compute_subdivisions[2]*z);
+                + edges[0]/compute_subdivisions[0]*x
+                + edges[1]/compute_subdivisions[1]*y
+                + edges[2]/compute_subdivisions[2]*z);
         break;
 
       default:
index 7c9d560cdc8767d0809815b5fe8113fc0c6d943c..6345627c66203db364807fc77c025af849258541 100644 (file)
@@ -51,8 +51,8 @@ void check_parallelepiped (bool colorize, bool log)
       break;
 
     case 2:
-      corners[0] = Point<dim> (0.0, 0.5);
-      corners[1] = Point<dim> (0.5, 0.0);
+      corners[0] = Point<dim> (0.5, 0.0);
+      corners[1] = Point<dim> (0.0, 0.5);
       break;
 
     case 3:
index de9c7740a05c618e23cb48090286e79f8698bbf5..da4467cc482213171727e7d236eeaebec13cc6d5 100644 (file)
@@ -71,8 +71,8 @@ void check_2d_parallelepiped_by_comparison (bool log)
   // build corners for this particular dim that are known to give the
   // same output order as parallelogram:
   Point<2> (corners) [2];
-  corners[0] = Point<2> (0.0, 0.5);
-  corners[1] = Point<2> (0.5, 0.0);
+  corners[0] = Point<2> (0.5, 0.0);
+  corners[1] = Point<2> (0.0, 0.5);
 
   Triangulation<2> triangulation_parallelepiped;
   GridGenerator::parallelepiped (triangulation_parallelepiped, corners, false);
index 382c2868ebd4bd73dbe726cdc1538b85e9d3028f..6a11233ba71c32db042e92e7e5331c0f4cdf1bad 100644 (file)
@@ -59,8 +59,8 @@ void check_nd_parallelepiped_by_comparison (bool log)
 
     case 2:
     {
-      corners[0] = Point<dim> (0.25, 0.50);
-      corners[1] = Point<dim> (0.50, 0.25);
+      corners[0] = Point<dim> (0.50, 0.25);
+      corners[1] = Point<dim> (0.25, 0.50);
       break;
     }
 
index c2b8a84d8c0da7141ec983396473790772a07948..25401aab5ef37679592703a49b511e96097b1158 100644 (file)
@@ -51,8 +51,8 @@ void check_subdivided_parallelepiped (bool colorize, bool log)
       break;
 
     case 2:
-      corners[0] = Point<dim> (0.25, 0.50);
-      corners[1] = Point<dim> (0.50, 0.25);
+      corners[0] = Point<dim> (0.50, 0.25);
+      corners[1] = Point<dim> (0.25, 0.50);
       break;
 
     case 3:
index 4ed8dc346f189f1fbeff353a37be3c693c8048e4..c070110d223b356ea7196c5776c283de96ce4d96 100644 (file)
@@ -62,8 +62,8 @@ void check (bool subdivide)
       break;
 
     case 2:
-      edges[0] = point<spacedim>(0.15, 0.50, 0.03);
-      edges[1] = point<spacedim>(0.70, 0.25, -0.01);
+      edges[0] = point<spacedim>(0.70, 0.25, -0.01);
+      edges[1] = point<spacedim>(0.15, 0.50, 0.03);
       break;
 
     case 3:
index 0b04927230aa5f105d8ca3d355041be79c533102..41ea526d7522b82b6f79dc82b28f31c48df81629 100644 (file)
@@ -16,32 +16,32 @@ DEAL::dim 1 spacedim 3
 
 DEAL::dim 2 spacedim 2
 0.100000 1.10000 0 0
-0.450000 1.22500 0 0
-0.600000 1.72500 0 0
-0.250000 1.60000 0 0
+0.800000 1.35000 0 0
+0.875000 1.60000 0 0
+0.175000 1.35000 0 0
 0.100000 1.10000 0 0
 
 
-0.450000 1.22500 0 0
-0.800000 1.35000 0 0
+0.175000 1.35000 0 0
+0.875000 1.60000 0 0
 0.950000 1.85000 0 0
-0.600000 1.72500 0 0
-0.450000 1.22500 0 0
+0.250000 1.60000 0 0
+0.175000 1.35000 0 0
 
 
 DEAL::dim 2 spacedim 3
 0.100000 1.10000 2.10000 0 0
-0.450000 1.22500 2.09500 0 0
-0.600000 1.72500 2.12500 0 0
-0.250000 1.60000 2.13000 0 0
+0.800000 1.35000 2.09000 0 0
+0.875000 1.60000 2.10500 0 0
+0.175000 1.35000 2.11500 0 0
 0.100000 1.10000 2.10000 0 0
 
 
-0.450000 1.22500 2.09500 0 0
-0.800000 1.35000 2.09000 0 0
+0.175000 1.35000 2.11500 0 0
+0.875000 1.60000 2.10500 0 0
 0.950000 1.85000 2.12000 0 0
-0.600000 1.72500 2.12500 0 0
-0.450000 1.22500 2.09500 0 0
+0.250000 1.60000 2.13000 0 0
+0.175000 1.35000 2.11500 0 0
 
 
 DEAL::dim 3 spacedim 3

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.