From 5814157bceeb8fbac44d6694fccad486c9605aef Mon Sep 17 00:00:00 2001
From: Luca Heltai <luca.heltai@sissa.it>
Date: Mon, 28 Jul 2014 15:16:36 +0200
Subject: [PATCH] Fixed documentation, and made default behavior do what we say
 it does.

Fixed compilation bug in manifold.cc

Added instantiations of ManifoldChart, and fixed a little bug.

Fixed instantiation bug, and tolerance error.
---
 include/deal.II/grid/manifold.h | 19 +++++++++++++++----
 source/grid/manifold.cc         | 27 ++++++++++++++++++++-------
 source/grid/manifold.inst.in    |  4 ++++
 3 files changed, 39 insertions(+), 11 deletions(-)

diff --git a/include/deal.II/grid/manifold.h b/include/deal.II/grid/manifold.h
index 4b6de55b67..2cde5218a6 100644
--- a/include/deal.II/grid/manifold.h
+++ b/include/deal.II/grid/manifold.h
@@ -70,8 +70,10 @@ namespace Manifolds {
  *   collection of points in @p spacedim dimension, and a collection of
  *   weights.
  *
- *   Internally, the get_new_point() function calls the project_to_manifold()
- *   function. This allows end users to only overload project_to_manifold().
+ *   Internally, the get_new_point() function calls the
+ *   project_to_manifold() function after computing the weighted
+ *   average of the quadrature poitns. This allows end users to only
+ *   overload project_to_manifold() for simple situations.
  *
  *   Should a finer control be necessary, then get_new_point() can be
  *   overloaded.  For backward compatibility, this function also
@@ -80,6 +82,13 @@ namespace Manifolds {
  *   FlatManifold<dim,spacedim>, allowing old user codes to keep using
  *   their boundary descriptors as Manifold<dim,spacedim> objects.
  *
+ *   The default behavior of these backward compatible interfaces is
+ *   to construct a Quadrature<spacedim> object containting the
+ *   vertices, midpoints of lines, and midpoints of quads with the
+ *   correct weight, and call get_new_point() with this quadrature. If
+ *   you need finer tuning for lines, quads or hexes, you can overload
+ *   any of the get_new_point_on_* functions. 
+ *
  *   FlatManifold is the specialization from which StraigthBoundary is
  *   derived, where the project_to_manifold() function is the identity.
  *
@@ -134,7 +143,9 @@ public:
    * become the new middle vertex of the two children of a regular
    * line. In 2D, this line is a line at the boundary, while in 3d, it
    * is bounding a face at the boundary (the lines therefore is also
-   * on the boundary). The default implementation of this function
+   * on the boundary).
+   *
+   * The default implementation of this function
    * passes its argument to the Manifolds::get_default_quadrature()
    * function, and then calls the
    * Manifold<dim,spacedim>::get_new_point() function. User derived
@@ -415,7 +426,7 @@ public:
    * The sub_manifold object is used to compute the average of the
    * points in the chart coordinates system. 
    */
-  FlatManifold<dim,spacedim> sub_manifold;
+  FlatManifold<dim,chartdim> sub_manifold;
 };
 
 
diff --git a/source/grid/manifold.cc b/source/grid/manifold.cc
index 0d9205e259..1685bbc1be 100644
--- a/source/grid/manifold.cc
+++ b/source/grid/manifold.cc
@@ -142,12 +142,24 @@ project_to_manifold (const std::vector<Point<spacedim> > &,
 template <int dim, int spacedim>
 Point<spacedim>
 Manifold<dim, spacedim>::
-get_new_point (const Quadrature<spacedim> &) const
+get_new_point (const Quadrature<spacedim> &quad) const
 {
-  Assert (false, ExcPureFunctionCalled());
-  return Point<spacedim>();
-}
+  const std::vector<Point<spacedim> > &surrounding_points = quad.get_points();
+  const std::vector<double> &weights = quad.get_weights();
+  Point<spacedim> p;
 
+#ifdef DEBUG
+  double sum=0;
+  for(unsigned int i=0; i<weights.size(); ++i)
+    sum+= weights[i];
+  Assert(std::abs(sum-1.0) < 1e-10, ExcMessage("Weights should sum to 1!"));
+#endif
+  
+  for(unsigned int i=0; i<surrounding_points.size(); ++i) 
+    p += surrounding_points[i]*weights[i];
+
+  return project_to_manifold(surrounding_points, p);
+}
 
 
 template <int dim, int spacedim>
@@ -323,7 +335,7 @@ get_new_point (const Quadrature<spacedim> &quad) const
       for(unsigned int d=0; d<spacedim; ++d) {
 	minP[d] = std::min(minP[d], surrounding_points[i][d]);
 	if(periodicity[d] > 0)
-	  Assert(surrounding_points[i][d] < periodicity[d],
+	  Assert(surrounding_points[i][d] < periodicity[d]+1e-10,
 		 ExcMessage("One of the points does not lye into the periodic box! Bailing out."));
       }
   
@@ -377,8 +389,9 @@ get_new_point (const Quadrature<spacedim> &quad) const
   
   for(unsigned int i=0; i<surrounding_points.size(); ++i) 
     chart_points[i] = pull_back(surrounding_points[i]);
-      
-  Point<chartdim> p_chart = sub_manifold.get_new_point(chart_points, weights);
+
+  Quadrature<chartdim> chart_quad(chart_points, weights);
+  Point<chartdim> p_chart = sub_manifold.get_new_point(chart_quad);
   
   return push_forward(p_chart);
 }
diff --git a/source/grid/manifold.inst.in b/source/grid/manifold.inst.in
index 892f3f92bf..a840f39247 100644
--- a/source/grid/manifold.inst.in
+++ b/source/grid/manifold.inst.in
@@ -21,6 +21,10 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS
 #if deal_II_dimension <= deal_II_space_dimension 	
     template class Manifold<deal_II_dimension, deal_II_space_dimension>;
     template class FlatManifold<deal_II_dimension, deal_II_space_dimension>;
+    
+    template class ManifoldChart<deal_II_dimension, deal_II_space_dimension, 1>;
+    template class ManifoldChart<deal_II_dimension, deal_II_space_dimension, 2>;
+    template class ManifoldChart<deal_II_dimension, deal_II_space_dimension, 3>;
 #endif
   }
 
-- 
2.39.5