From: David Wells Date: Mon, 9 May 2016 01:20:14 +0000 (-0400) Subject: Raise exceptions for invalid input. X-Git-Tag: v8.5.0-rc1~1033^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=1c209aff4467e191f2fed9985a64cde2eb59cb21;p=dealii.git Raise exceptions for invalid input. 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. --- diff --git a/include/deal.II/grid/grid_generator.h b/include/deal.II/grid/grid_generator.h index c70c33ca42..cace929cda 100644 --- a/include/deal.II/grid/grid_generator.h +++ b/include/deal.II/grid/grid_generator.h @@ -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."); ///@} } diff --git a/source/grid/grid_generator.cc b/source/grid/grid_generator.cc index 341b82f339..20354c17ee 100644 --- a/source/grid/grid_generator.cc +++ b/source/grid/grid_generator.cc @@ -856,9 +856,6 @@ namespace GridGenerator const std::vector &subdivisions, const bool colorize) { - // The provided edges and subdivisions may not be valid, so possibly - // recompute both of them. - std_cxx11::array, dim> compute_edges = edges; std::vector 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 (0.0, 0.5); - corners[1] = Point (0.5, 0.0); + corners[0] = Point (0.5, 0.0); + corners[1] = Point (0.0, 0.5); break; case 3: diff --git a/tests/grid/grid_parallelepiped_01.cc b/tests/grid/grid_parallelepiped_01.cc index de9c7740a0..da4467cc48 100644 --- a/tests/grid/grid_parallelepiped_01.cc +++ b/tests/grid/grid_parallelepiped_01.cc @@ -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); diff --git a/tests/grid/grid_parallelepiped_02.cc b/tests/grid/grid_parallelepiped_02.cc index 382c2868eb..6a11233ba7 100644 --- a/tests/grid/grid_parallelepiped_02.cc +++ b/tests/grid/grid_parallelepiped_02.cc @@ -59,8 +59,8 @@ void check_nd_parallelepiped_by_comparison (bool log) case 2: { - corners[0] = Point (0.25, 0.50); - corners[1] = Point (0.50, 0.25); + corners[0] = Point (0.50, 0.25); + corners[1] = Point (0.25, 0.50); break; } diff --git a/tests/grid/grid_parallelepiped_03.cc b/tests/grid/grid_parallelepiped_03.cc index c2b8a84d8c..25401aab5e 100644 --- a/tests/grid/grid_parallelepiped_03.cc +++ b/tests/grid/grid_parallelepiped_03.cc @@ -51,8 +51,8 @@ void check_subdivided_parallelepiped (bool colorize, bool log) break; case 2: - corners[0] = Point (0.25, 0.50); - corners[1] = Point (0.50, 0.25); + corners[0] = Point (0.50, 0.25); + corners[1] = Point (0.25, 0.50); break; case 3: diff --git a/tests/grid/grid_parallelepiped_05.cc b/tests/grid/grid_parallelepiped_05.cc index 4ed8dc346f..c070110d22 100644 --- a/tests/grid/grid_parallelepiped_05.cc +++ b/tests/grid/grid_parallelepiped_05.cc @@ -62,8 +62,8 @@ void check (bool subdivide) break; case 2: - edges[0] = point(0.15, 0.50, 0.03); - edges[1] = point(0.70, 0.25, -0.01); + edges[0] = point(0.70, 0.25, -0.01); + edges[1] = point(0.15, 0.50, 0.03); break; case 3: diff --git a/tests/grid/grid_parallelepiped_05.output b/tests/grid/grid_parallelepiped_05.output index 0b04927230..41ea526d75 100644 --- a/tests/grid/grid_parallelepiped_05.output +++ b/tests/grid/grid_parallelepiped_05.output @@ -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