From 37976bf3ea8b5238c3741eb1e1bedae67821c4a7 Mon Sep 17 00:00:00 2001
From: bangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Fri, 5 Apr 2013 17:47:28 +0000
Subject: [PATCH] Fix bug #4: Triangulation::distort_random accidentally only
 ever moved points in positive coordinate directions.

git-svn-id: https://svn.dealii.org/trunk@29203 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/doc/news/changes.h  | 11 ++++++++++-
 deal.II/source/grid/tria.cc | 22 ++++++++++++----------
 2 files changed, 22 insertions(+), 11 deletions(-)

diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h
index c6c59816ac..ae9c06b683 100644
--- a/deal.II/doc/news/changes.h
+++ b/deal.II/doc/news/changes.h
@@ -88,13 +88,22 @@ this function.
 <h3>Specific improvements</h3>
 
 <ol>
+<li> Fixed: GridTools::distort_random (previously called Triangulation::distort_random)
+had a bug where points were only ever moved in <i>positive</i> coordinate
+directions rather than with uniform probability in either direction. The 1d
+implementation also had the problem that it did not move vertices if the
+<i>cell</i> they were on was at the boundary, even if the <i>vertex</i>
+itself was not. All of these problems are now fixed.
+<br>
+(Wolfgang Bangerth, 2013/4/5)
+</li>
+
 <li> New: There is a class VectorFunctionFromTensorFunction that converts
 between objects of type TensorFunction and Function.
 <br>
 (Spencer Patty, 2013/4/2)
 </li>
 
-
 <li> Fixed: The ParameterHandler class could not deal with parameters named
 <code>"value"</code> (and a few other names). This is now fixed.
 <br>
diff --git a/deal.II/source/grid/tria.cc b/deal.II/source/grid/tria.cc
index c6ba3b09de..4b7607565d 100644
--- a/deal.II/source/grid/tria.cc
+++ b/deal.II/source/grid/tria.cc
@@ -2793,20 +2793,12 @@ namespace internal
 
         std::vector<double> minimal_length (triangulation.vertices.size(),
                                             almost_infinite_length);
-        // also note if a vertex is at
-        // the boundary
         std::vector<bool>   at_boundary (triangulation.vertices.size(), false);
 
         for (typename Triangulation<dim,spacedim>::active_line_iterator
              line=triangulation.begin_active_line();
              line != triangulation.end_line(); ++line)
           {
-            if (keep_boundary && line->at_boundary())
-              {
-                at_boundary[line->vertex_index(0)] = true;
-                at_boundary[line->vertex_index(1)] = true;
-              }
-
             minimal_length[line->vertex_index(0)]
               = std::min(line->diameter(),
                          minimal_length[line->vertex_index(0)]);
@@ -2815,6 +2807,16 @@ namespace internal
                          minimal_length[line->vertex_index(1)]);
           }
 
+        // also note if a vertex is at the boundary if we are asked to
+        // keep boundary vertices untouched
+	if (keep_boundary)
+	  for (typename Triangulation<dim,spacedim>::active_line_iterator
+		 line=triangulation.begin_active_line();
+	       line != triangulation.end_line(); ++line)
+	    for (unsigned int vertex=0; vertex<2; ++vertex)
+	      if (line->at_boundary(vertex) == true)
+		at_boundary[line->vertex_index(vertex)] = true;
+
 
         const unsigned int n_vertices = triangulation.vertices.size();
         Point<spacedim> shift_vector;
@@ -2831,7 +2833,7 @@ namespace internal
             // first compute a random shift
             // vector
             for (unsigned int d=0; d<spacedim; ++d)
-              shift_vector(d) = std::rand()*1.0/RAND_MAX;
+              shift_vector(d) = std::rand()*2.0/RAND_MAX-1;
 
             shift_vector *= factor * minimal_length[vertex] /
                             std::sqrt(shift_vector.square());
@@ -2927,7 +2929,7 @@ namespace internal
             // first compute a random shift
             // vector
             for (unsigned int d=0; d<spacedim; ++d)
-              shift_vector(d) = std::rand()*1.0/RAND_MAX;
+              shift_vector(d) = std::rand()*2.0/RAND_MAX-1;
 
             shift_vector *= factor * minimal_length[vertex] /
                             std::sqrt(shift_vector.square());
-- 
2.39.5