From ccb56e7fb59f4a4892f27fcff876846fa631b213 Mon Sep 17 00:00:00 2001
From: David Wells <wellsd2@rpi.edu>
Date: Sat, 7 May 2016 15:13:35 -0400
Subject: [PATCH] Fix a bug in the parallelepiped mesh generator.

This function generated cells with negative volume due to an orientation
issue. In particular, the function used to generate cells (see
tests/grid/grid_parallelepiped.output) with vertices (given here in
lexical order)

0 0
0 0.5
0.5 0.5
0.5 0

which results in a twisted cell, instead of

0 0
0.5 0
0.5 0.5
0 0.5

which is correct.
---
 doc/news/changes.h                       |   7 +
 source/grid/grid_generator.cc            |  73 ++++----
 tests/grid/grid_parallelepiped.output    |   4 +-
 tests/grid/grid_parallelepiped_03.output | 204 +++++++++++------------
 tests/grid/grid_parallelepiped_05.output |  16 +-
 5 files changed, 164 insertions(+), 140 deletions(-)

diff --git a/doc/news/changes.h b/doc/news/changes.h
index 7718a459df..c62ea5a10d 100644
--- a/doc/news/changes.h
+++ b/doc/news/changes.h
@@ -249,6 +249,13 @@ inconvenience this causes.
  (Joscha Gedicke, 2016/05/07)
  </li>
 
+ <li> Fixed: The function GridGenerator::subdivided_parallelepiped and its
+ variants could generate two dimensional meshes with cells that had negative
+ Jacobians. This has been fixed.
+ <br>
+ (David Wells, 2016/05/07)
+ </li>
+
  <li> New: Added function GridOut::write_mesh_per_processor_as_vtu. This allows 
  the visualization of a parallel finite element mesh that can be separated into each 
  processor's owned and ghost cells. It also allows for the visualization of each level
diff --git a/source/grid/grid_generator.cc b/source/grid/grid_generator.cc
index d1aafe9c33..341b82f339 100644
--- a/source/grid/grid_generator.cc
+++ b/source/grid/grid_generator.cc
@@ -856,26 +856,43 @@ namespace GridGenerator
                              const std::vector<unsigned int> &subdivisions,
                              const bool colorize)
   {
-    if (subdivisions.size()==0)
+    // 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)
       {
-        std::vector<unsigned int> new_subdivisions(dim, 1);
-        subdivided_parallelepiped<dim,spacedim>(tria, origin, edges, new_subdivisions, colorize);
-        return;
+        compute_subdivisions.resize(dim, 1);
       }
 
-    Assert(subdivisions.size()==dim, ExcMessage(""));
-
+    Assert(compute_subdivisions.size()==dim,
+           ExcMessage("One subdivision must be provided for each dimension."));
     // check subdivisions
     for (unsigned int i=0; i<dim; ++i)
       {
-        Assert (subdivisions[i]>0, ExcInvalidRepetitions(subdivisions[i]));
-        Assert (edges[i].norm()>0, ExcMessage("Edges in subdivided_parallelepiped() must not be degenerate."));
+        Assert (compute_subdivisions[i]>0, ExcInvalidRepetitions(subdivisions[i]));
+        Assert (edges[i].norm()>0,
+                ExcMessage("Edges in subdivided_parallelepiped() must not be degenerate."));
+      }
+
+    /*
+     * Verify that the vectors are oriented in a counter clockwise direction in 2D.
+     */
+    if (dim == 2)
+      {
+        const double plane_normal = edges[0][0]*edges[1][1] - edges[0][1]*edges[1][0];
+        if (plane_normal < 0.0)
+          {
+            std::swap(compute_edges[0], compute_edges[1]);
+            std::swap(compute_subdivisions[0], compute_subdivisions[1]);
+          }
       }
 
+
     // Check corners do not overlap (unique)
     for (unsigned int i=0; i<dim; ++i)
       for (unsigned int j=i+1; j<dim; ++j)
-        Assert ((edges[i]!=edges[j]),
+        Assert ((compute_edges[i]!=compute_edges[j]),
                 ExcMessage ("Degenerate edges of subdivided_parallelepiped encountered."));
 
     // Create a list of points
@@ -884,27 +901,27 @@ namespace GridGenerator
     switch (dim)
       {
       case 1:
-        for (unsigned int x=0; x<=subdivisions[0]; ++x)
-          points.push_back (origin + edges[0]/subdivisions[0]*x);
+        for (unsigned int x=0; x<=compute_subdivisions[0]; ++x)
+          points.push_back (origin + compute_edges[0]/compute_subdivisions[0]*x);
         break;
 
       case 2:
-        for (unsigned int y=0; y<=subdivisions[1]; ++y)
-          for (unsigned int x=0; x<=subdivisions[0]; ++x)
+        for (unsigned int y=0; y<=compute_subdivisions[1]; ++y)
+          for (unsigned int x=0; x<=compute_subdivisions[0]; ++x)
             points.push_back (origin
-                              + edges[0]/subdivisions[0]*x
-                              + edges[1]/subdivisions[1]*y);
+                              + compute_edges[0]/compute_subdivisions[0]*x
+                              + compute_edges[1]/compute_subdivisions[1]*y);
         break;
 
       case 3:
-        for (unsigned int z=0; z<=subdivisions[2]; ++z)
-          for (unsigned int y=0; y<=subdivisions[1]; ++y)
-            for (unsigned int x=0; x<=subdivisions[0]; ++x)
+        for (unsigned int z=0; z<=compute_subdivisions[2]; ++z)
+          for (unsigned int y=0; y<=compute_subdivisions[1]; ++y)
+            for (unsigned int x=0; x<=compute_subdivisions[0]; ++x)
               points.push_back (
                 origin
-                + edges[0]/subdivisions[0]*x
-                + edges[1]/subdivisions[1]*y
-                + edges[2]/subdivisions[2]*z);
+                + compute_edges[0]/compute_subdivisions[0]*x
+                + compute_edges[1]/compute_subdivisions[1]*y
+                + compute_edges[2]/compute_subdivisions[2]*z);
         break;
 
       default:
@@ -914,14 +931,14 @@ namespace GridGenerator
     // Prepare cell data
     unsigned int n_cells = 1;
     for (unsigned int i=0; i<dim; ++i)
-      n_cells *= subdivisions[i];
+      n_cells *= compute_subdivisions[i];
     std::vector<CellData<dim> > cells (n_cells);
 
     // Create fixed ordering of
     switch (dim)
       {
       case 1:
-        for (unsigned int x=0; x<subdivisions[0]; ++x)
+        for (unsigned int x=0; x<compute_subdivisions[0]; ++x)
           {
             cells[x].vertices[0] = x;
             cells[x].vertices[1] = x+1;
@@ -934,8 +951,8 @@ namespace GridGenerator
       case 2:
       {
         // Shorthand
-        const unsigned int n_dy = subdivisions[1];
-        const unsigned int n_dx = subdivisions[0];
+        const unsigned int n_dy = compute_subdivisions[1];
+        const unsigned int n_dx = compute_subdivisions[0];
 
         for (unsigned int y=0; y<n_dy; ++y)
           for (unsigned int x=0; x<n_dx; ++x)
@@ -955,9 +972,9 @@ namespace GridGenerator
       case 3:
       {
         // Shorthand
-        const unsigned int n_dz = subdivisions[2];
-        const unsigned int n_dy = subdivisions[1];
-        const unsigned int n_dx = subdivisions[0];
+        const unsigned int n_dz = compute_subdivisions[2];
+        const unsigned int n_dy = compute_subdivisions[1];
+        const unsigned int n_dx = compute_subdivisions[0];
 
         for (unsigned int z=0; z<n_dz; ++z)
           for (unsigned int y=0; y<n_dy; ++y)
diff --git a/tests/grid/grid_parallelepiped.output b/tests/grid/grid_parallelepiped.output
index f0422e944d..0d29404628 100644
--- a/tests/grid/grid_parallelepiped.output
+++ b/tests/grid/grid_parallelepiped.output
@@ -3,9 +3,9 @@
 
 
 0 0 0 0
-0 0.5 0 0
-0.5 0.5 0 0
 0.5 0 0 0
+0.5 0.5 0 0
+0 0.5 0 0
 0 0 0 0
 
 
diff --git a/tests/grid/grid_parallelepiped_03.output b/tests/grid/grid_parallelepiped_03.output
index c83641ec48..7d23a5f235 100644
--- a/tests/grid/grid_parallelepiped_03.output
+++ b/tests/grid/grid_parallelepiped_03.output
@@ -11,177 +11,177 @@
 
 
 0 0 0 0
-0.05 0.1 0 0
-0.15 0.15 0 0
 0.1 0.05 0 0
+0.15 0.15 0 0
+0.05 0.1 0 0
 0 0 0 0
 
 
-0.05 0.1 0 0
-0.1 0.2 0 0
-0.2 0.25 0 0
+0.1 0.05 0 0
+0.2 0.1 0 0
+0.25 0.2 0 0
 0.15 0.15 0 0
-0.05 0.1 0 0
+0.1 0.05 0 0
 
 
-0.1 0.2 0 0
-0.15 0.3 0 0
-0.25 0.35 0 0
-0.2 0.25 0 0
-0.1 0.2 0 0
+0.2 0.1 0 0
+0.3 0.15 0 0
+0.35 0.25 0 0
+0.25 0.2 0 0
+0.2 0.1 0 0
 
 
-0.15 0.3 0 0
-0.2 0.4 0 0
-0.3 0.45 0 0
-0.25 0.35 0 0
-0.15 0.3 0 0
+0.3 0.15 0 0
+0.4 0.2 0 0
+0.45 0.3 0 0
+0.35 0.25 0 0
+0.3 0.15 0 0
 
 
-0.2 0.4 0 0
-0.25 0.5 0 0
-0.35 0.55 0 0
-0.3 0.45 0 0
-0.2 0.4 0 0
+0.4 0.2 0 0
+0.5 0.25 0 0
+0.55 0.35 0 0
+0.45 0.3 0 0
+0.4 0.2 0 0
 
 
-0.1 0.05 0 0
+0.05 0.1 0 0
 0.15 0.15 0 0
-0.25 0.2 0 0
-0.2 0.1 0 0
-0.1 0.05 0 0
+0.2 0.25 0 0
+0.1 0.2 0 0
+0.05 0.1 0 0
 
 
 0.15 0.15 0 0
-0.2 0.25 0 0
-0.3 0.3 0 0
 0.25 0.2 0 0
+0.3 0.3 0 0
+0.2 0.25 0 0
 0.15 0.15 0 0
 
 
-0.2 0.25 0 0
-0.25 0.35 0 0
-0.35 0.4 0 0
+0.25 0.2 0 0
+0.35 0.25 0 0
+0.4 0.35 0 0
 0.3 0.3 0 0
-0.2 0.25 0 0
+0.25 0.2 0 0
 
 
-0.25 0.35 0 0
-0.3 0.45 0 0
-0.4 0.5 0 0
-0.35 0.4 0 0
-0.25 0.35 0 0
+0.35 0.25 0 0
+0.45 0.3 0 0
+0.5 0.4 0 0
+0.4 0.35 0 0
+0.35 0.25 0 0
 
 
-0.3 0.45 0 0
-0.35 0.55 0 0
-0.45 0.6 0 0
-0.4 0.5 0 0
-0.3 0.45 0 0
+0.45 0.3 0 0
+0.55 0.35 0 0
+0.6 0.45 0 0
+0.5 0.4 0 0
+0.45 0.3 0 0
 
 
-0.2 0.1 0 0
-0.25 0.2 0 0
-0.35 0.25 0 0
-0.3 0.15 0 0
-0.2 0.1 0 0
+0.1 0.2 0 0
+0.2 0.25 0 0
+0.25 0.35 0 0
+0.15 0.3 0 0
+0.1 0.2 0 0
 
 
-0.25 0.2 0 0
+0.2 0.25 0 0
 0.3 0.3 0 0
-0.4 0.35 0 0
-0.35 0.25 0 0
-0.25 0.2 0 0
+0.35 0.4 0 0
+0.25 0.35 0 0
+0.2 0.25 0 0
 
 
 0.3 0.3 0 0
-0.35 0.4 0 0
-0.45 0.45 0 0
 0.4 0.35 0 0
+0.45 0.45 0 0
+0.35 0.4 0 0
 0.3 0.3 0 0
 
 
-0.35 0.4 0 0
-0.4 0.5 0 0
-0.5 0.55 0 0
+0.4 0.35 0 0
+0.5 0.4 0 0
+0.55 0.5 0 0
 0.45 0.45 0 0
-0.35 0.4 0 0
+0.4 0.35 0 0
 
 
-0.4 0.5 0 0
-0.45 0.6 0 0
-0.55 0.65 0 0
-0.5 0.55 0 0
-0.4 0.5 0 0
+0.5 0.4 0 0
+0.6 0.45 0 0
+0.65 0.55 0 0
+0.55 0.5 0 0
+0.5 0.4 0 0
 
 
-0.3 0.15 0 0
-0.35 0.25 0 0
-0.45 0.3 0 0
-0.4 0.2 0 0
-0.3 0.15 0 0
+0.15 0.3 0 0
+0.25 0.35 0 0
+0.3 0.45 0 0
+0.2 0.4 0 0
+0.15 0.3 0 0
 
 
-0.35 0.25 0 0
-0.4 0.35 0 0
-0.5 0.4 0 0
-0.45 0.3 0 0
-0.35 0.25 0 0
+0.25 0.35 0 0
+0.35 0.4 0 0
+0.4 0.5 0 0
+0.3 0.45 0 0
+0.25 0.35 0 0
 
 
-0.4 0.35 0 0
+0.35 0.4 0 0
 0.45 0.45 0 0
-0.55 0.5 0 0
-0.5 0.4 0 0
-0.4 0.35 0 0
+0.5 0.55 0 0
+0.4 0.5 0 0
+0.35 0.4 0 0
 
 
 0.45 0.45 0 0
-0.5 0.55 0 0
-0.6 0.6 0 0
 0.55 0.5 0 0
+0.6 0.6 0 0
+0.5 0.55 0 0
 0.45 0.45 0 0
 
 
-0.5 0.55 0 0
-0.55 0.65 0 0
-0.65 0.7 0 0
+0.55 0.5 0 0
+0.65 0.55 0 0
+0.7 0.65 0 0
 0.6 0.6 0 0
-0.5 0.55 0 0
+0.55 0.5 0 0
 
 
-0.4 0.2 0 0
-0.45 0.3 0 0
-0.55 0.35 0 0
-0.5 0.25 0 0
-0.4 0.2 0 0
+0.2 0.4 0 0
+0.3 0.45 0 0
+0.35 0.55 0 0
+0.25 0.5 0 0
+0.2 0.4 0 0
 
 
-0.45 0.3 0 0
-0.5 0.4 0 0
-0.6 0.45 0 0
-0.55 0.35 0 0
-0.45 0.3 0 0
+0.3 0.45 0 0
+0.4 0.5 0 0
+0.45 0.6 0 0
+0.35 0.55 0 0
+0.3 0.45 0 0
 
 
-0.5 0.4 0 0
-0.55 0.5 0 0
-0.65 0.55 0 0
-0.6 0.45 0 0
-0.5 0.4 0 0
+0.4 0.5 0 0
+0.5 0.55 0 0
+0.55 0.65 0 0
+0.45 0.6 0 0
+0.4 0.5 0 0
 
 
-0.55 0.5 0 0
+0.5 0.55 0 0
 0.6 0.6 0 0
-0.7 0.65 0 0
-0.65 0.55 0 0
-0.55 0.5 0 0
+0.65 0.7 0 0
+0.55 0.65 0 0
+0.5 0.55 0 0
 
 
 0.6 0.6 0 0
-0.65 0.7 0 0
-0.75 0.75 0 0
 0.7 0.65 0 0
+0.75 0.75 0 0
+0.65 0.7 0 0
 0.6 0.6 0 0
 
 
diff --git a/tests/grid/grid_parallelepiped_05.output b/tests/grid/grid_parallelepiped_05.output
index 753252c3d3..0b04927230 100644
--- a/tests/grid/grid_parallelepiped_05.output
+++ b/tests/grid/grid_parallelepiped_05.output
@@ -16,31 +16,31 @@ DEAL::dim 1 spacedim 3
 
 DEAL::dim 2 spacedim 2
 0.100000 1.10000 0 0
-0.250000 1.60000 0 0
-0.600000 1.72500 0 0
 0.450000 1.22500 0 0
+0.600000 1.72500 0 0
+0.250000 1.60000 0 0
 0.100000 1.10000 0 0
 
 
 0.450000 1.22500 0 0
-0.600000 1.72500 0 0
-0.950000 1.85000 0 0
 0.800000 1.35000 0 0
+0.950000 1.85000 0 0
+0.600000 1.72500 0 0
 0.450000 1.22500 0 0
 
 
 DEAL::dim 2 spacedim 3
 0.100000 1.10000 2.10000 0 0
-0.250000 1.60000 2.13000 0 0
-0.600000 1.72500 2.12500 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.100000 1.10000 2.10000 0 0
 
 
 0.450000 1.22500 2.09500 0 0
-0.600000 1.72500 2.12500 0 0
-0.950000 1.85000 2.12000 0 0
 0.800000 1.35000 2.09000 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
 
 
-- 
2.39.5