From: Luca Heltai <luca.heltai@sissa.it>
Date: Sun, 20 Mar 2016 22:49:22 +0000 (+0100)
Subject: Fixed @bangerth comments.
X-Git-Tag: v8.5.0-rc1~1176^2~2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=4e96c13c2c286a3d9bc14da7d3637c9a3efe48b6;p=dealii.git

Fixed @bangerth comments.
---

diff --git a/source/grid/manifold.cc b/source/grid/manifold.cc
index e3fb14c7f7..129b7d1555 100644
--- a/source/grid/manifold.cc
+++ b/source/grid/manifold.cc
@@ -74,7 +74,8 @@ normal_vector (const typename Triangulation<2, 2>::face_iterator &face,
   Tensor<1,spacedim> tangent = ((p-face->vertex(0)).norm_square() > (p-face->vertex(1)).norm_square() ?
                                 get_tangent_vector(p, face->vertex(0)) :
                                 -get_tangent_vector(p, face->vertex(1)));
-  return cross_product_2d(tangent);
+  Tensor<1,spacedim> normal = cross_product_2d(tangent);
+  return normal/normal.norm();
 }
 
 template<>
@@ -87,16 +88,19 @@ normal_vector (const typename Triangulation<3, 3>::face_iterator &face,
   Tensor<1,spacedim> t1,t2;
 
   // Take the difference between p and all four vertices
-  Tensor<1,spacedim> dp[4];
-  double dpns[4];
-  int min_index=-1;
-  double min_distance = 0;
-  for (unsigned int i=0; i<4; ++i)
+  int min_index=0;
+  Tensor<1,spacedim> dp = p-face->vertex(0);
+  double min_distance = dp.norm_square();
+
+  for (unsigned int i=1; i<4; ++i)
     {
-      dp[i] = p-face->vertex(i);
-      dpns[i] = dp[i].norm_square();
-      min_index = (min_index == -1 ? (int)i : dpns[i] < min_distance ? i : min_index);
-      min_distance = dpns[min_index];
+      dp = p-face->vertex(i);
+      double distance = dp.norm_square();
+      if (distance < min_distance)
+        {
+          min_index = i;
+          min_distance = distance;
+        }
     }
   // Verify we have a valid vertex index
   AssertIndexRange(min_index, 4);
@@ -153,7 +157,8 @@ normal_vector (const typename Triangulation<3, 3>::face_iterator &face,
           break;
         }
     }
-  return cross_product_3d(t1,t2);
+  Tensor<1,spacedim> normal = cross_product_3d(t1,t2);
+  return normal/normal.norm();
 }
 
 
@@ -168,6 +173,48 @@ normal_vector (const typename Triangulation<dim, spacedim>::face_iterator &face,
 }
 
 
+
+template <>
+void
+Manifold<2, 2>::
+get_normals_at_vertices (const typename Triangulation<2, 2>::face_iterator &face,
+                         FaceVertexNormals &n) const
+{
+  n[0] = cross_product_2d(get_tangent_vector(face->vertex(0), face->vertex(1)));
+  n[1] = -cross_product_2d(get_tangent_vector(face->vertex(1), face->vertex(0)));
+
+  n[0] /= n[0].norm();
+  n[1] /= n[1].norm();
+}
+
+
+template <>
+void
+Manifold<3, 3>::
+get_normals_at_vertices (const typename Triangulation<3, 3>::face_iterator &face,
+                         FaceVertexNormals &n) const
+{
+  n[0] = cross_product_3d
+         (get_tangent_vector(face->vertex(0), face->vertex(1)),
+          get_tangent_vector(face->vertex(0), face->vertex(2)));
+
+  n[1] = cross_product_3d
+         (get_tangent_vector(face->vertex(1), face->vertex(3)),
+          get_tangent_vector(face->vertex(1), face->vertex(0)));
+
+  n[2] = cross_product_3d
+         (get_tangent_vector(face->vertex(2), face->vertex(0)),
+          get_tangent_vector(face->vertex(2), face->vertex(3)));
+
+  n[3] = cross_product_3d
+         (get_tangent_vector(face->vertex(3), face->vertex(2)),
+          get_tangent_vector(face->vertex(3), face->vertex(1)));
+
+  for (unsigned int i=0; i<4; ++i)
+    n[i] /=n[i].norm();
+}
+
+
 template <int dim, int spacedim>
 void
 Manifold<dim, spacedim>::
diff --git a/tests/manifold/spherical_manifold_07.cc b/tests/manifold/spherical_manifold_07.cc
index f92a567fdf..a058cb6fa3 100644
--- a/tests/manifold/spherical_manifold_07.cc
+++ b/tests/manifold/spherical_manifold_07.cc
@@ -13,7 +13,7 @@
 //
 // ---------------------------------------------------------------------
 
-// test GridTools::torus() and TorusManifold, output visually checked
+// test get_normals_at_vertices for a SphericalManifold.
 
 #include "../tests.h"
 #include <deal.II/base/logstream.h>
@@ -46,7 +46,8 @@ void test ()
       << "set view equal xyz" << std::endl
       << "set size ratio -1" << std::endl
       << "set multiplot" << std::endl
-      << "splot '-' with vectors " << std::endl;
+      << (dim == 3 ? "s" : "")
+      << "plot '-' with vectors " << std::endl;
 
   for (typename Triangulation<dim,spacedim>::active_cell_iterator
        cell = triangulation.begin_active();
diff --git a/tests/manifold/spherical_manifold_07.output b/tests/manifold/spherical_manifold_07.output
index 9bd8aaac63..88cdbd3fb7 100644
--- a/tests/manifold/spherical_manifold_07.output
+++ b/tests/manifold/spherical_manifold_07.output
@@ -2,23 +2,23 @@
 set view equal xyz
 set size ratio -1
 set multiplot
-splot '-' with vectors 
--0.707107 -0.707107 0.0707107 0.0707107
--1.83697e-16 -1.00000 1.83697e-17 0.100000
--1.83697e-16 -1.00000 1.83697e-17 0.100000
-0.707107 -0.707107 -0.0707107 0.0707107
+plot '-' with vectors 
 -0.707107 -0.707107 -0.0707107 -0.0707107
--1.00000 1.22465e-16 -0.100000 1.22465e-17
--1.00000 1.22465e-16 -0.100000 1.22465e-17
--0.707107 0.707107 -0.0707107 0.0707107
-0.707107 -0.707107 -0.0707107 0.0707107
-1.00000 0.00000 -0.100000 0.00000
-1.00000 0.00000 -0.100000 0.00000
-0.707107 0.707107 -0.0707107 -0.0707107
--0.707107 0.707107 -0.0707107 0.0707107
-6.12323e-17 1.00000 6.12323e-18 0.100000
-6.12323e-17 1.00000 6.12323e-18 0.100000
+-1.83697e-16 -1.00000 -1.83697e-17 -0.100000
+-1.83697e-16 -1.00000 -1.83697e-17 -0.100000
+0.707107 -0.707107 0.0707107 -0.0707107
+-0.707107 -0.707107 0.0707107 0.0707107
+-1.00000 1.22465e-16 0.100000 -1.22465e-17
+-1.00000 1.22465e-16 0.100000 -1.22465e-17
+-0.707107 0.707107 0.0707107 -0.0707107
+0.707107 -0.707107 0.0707107 -0.0707107
+1.00000 0.00000 0.100000 0.00000
+1.00000 0.00000 0.100000 0.00000
 0.707107 0.707107 0.0707107 0.0707107
+-0.707107 0.707107 0.0707107 -0.0707107
+6.12323e-17 1.00000 -6.12323e-18 -0.100000
+6.12323e-17 1.00000 -6.12323e-18 -0.100000
+0.707107 0.707107 -0.0707107 -0.0707107
 e
 set view equal xyz
 set size ratio -1