From aef23f358d6245828ce5f3d6ede028e932d50516 Mon Sep 17 00:00:00 2001
From: heltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Sat, 4 Jan 2014 00:07:17 +0000
Subject: [PATCH] Made intermediate points of manifold lines and quads equal to
 new implementation of Boundary (Gauss Lobatto).

git-svn-id: https://svn.dealii.org/branches/branch_manifold_id@32163 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/include/deal.II/grid/manifold.h | 29 ++++++++++++++
 deal.II/source/grid/manifold.cc         | 52 ++++++++++++++++++++-----
 2 files changed, 71 insertions(+), 10 deletions(-)

diff --git a/deal.II/include/deal.II/grid/manifold.h b/deal.II/include/deal.II/grid/manifold.h
index 553d9bb0eb..d088f0e813 100644
--- a/deal.II/include/deal.II/grid/manifold.h
+++ b/deal.II/include/deal.II/grid/manifold.h
@@ -22,6 +22,9 @@
 
 #include <deal.II/base/config.h>
 #include <deal.II/base/subscriptor.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/thread_management.h>
+#include <deal.II/base/subscriptor.h>
 #include <deal.II/base/point.h>
 
 DEAL_II_NAMESPACE_OPEN
@@ -145,7 +148,33 @@ public:
     virtual
     Point<spacedim> normal_vector(const Point<spacedim> point) const;
     
+  protected:
+
+  /**
+   * Returns the support points of the Gauss-Lobatto quadrature
+   * formula used for intermediate points.
+   *
+   * @note Since the manifold description is closely tied to the unit
+   * cell support points of MappingQ, new boundary descriptions need
+   * to explicitly use these Gauss-Lobatto points and not equidistant
+   * points.
+   */
+    const std::vector<Point<1> > &
+    get_line_support_points (const unsigned int n_intermediate_points) const;
+    
     
+  private:
+
+				     /**
+				      * Point generator for the intermediate points on a boundary.
+				      */
+    mutable std::vector<std_cxx1x::shared_ptr<QGaussLobatto<1> > > points;
+
+
+				     /**
+				      * Mutex for protecting the points array.
+				      */
+    mutable Threads::Mutex mutex;
 };
 
 
diff --git a/deal.II/source/grid/manifold.cc b/deal.II/source/grid/manifold.cc
index 7a6de0639f..fb02d1bb20 100644
--- a/deal.II/source/grid/manifold.cc
+++ b/deal.II/source/grid/manifold.cc
@@ -72,10 +72,13 @@ Manifold<spacedim>::get_intermediate_points(std::vector<Point<spacedim> > &point
     {
       case 2:
       {
-	const double dx=1./(n+1);
-	double x = dx;
-	for(unsigned int i=0; i<n; ++i, x+=dx) 
+        // Use interior points of QGaussLobatto quadrature formula
+        // support points for consistency with MappingQ
+	const std::vector<Point<1> > &line_points = this->get_line_support_points(n);
+	
+	for(unsigned int i=0; i<n; ++i) 
 	  {
+	    const double x = line_points[i+1][0];
 	    w[1] = x; w[0] = (1-x);
 	    points[i] = get_new_point(surrounding_points, w);
 	  }
@@ -89,15 +92,20 @@ Manifold<spacedim>::get_intermediate_points(std::vector<Point<spacedim> > &point
 	  static_cast<unsigned int>(std::sqrt(static_cast<double>(n)));
 					 // is n a square number
 	Assert(m*m==n, ExcInternalError());
-	const double ds=1./(m+1);
-	double y=ds;
-	for (unsigned int i=0; i<m; ++i, y+=ds)
+
+	const std::vector<Point<1> > &line_points = this->get_line_support_points(m);
+
+	for (unsigned int i=0; i<m; ++i)
 	  {
-	    double x=ds;
-	    for (unsigned int j=0; j<m; ++j, x+=ds)
+	    const double y=line_points[1+i][0];
+	    for (unsigned int j=0; j<m; ++j)
 	      {
-		w[0] =   (1-x); w[1] = x*(1-y);
-		w[2] = y*(1-x); w[3] = x*y;
+		const double x=line_points[1+j][0];
+
+		w[0] = (1-x)*(1-y);
+		w[1] =     x*(1-y);
+		w[2] = (1-x)*y    ;
+		w[3] =     x*y    ;
 		points[i*m+j]=get_new_point(surrounding_points, w);
 	      }
 	  }
@@ -114,6 +122,30 @@ Manifold<spacedim>::get_intermediate_points(std::vector<Point<spacedim> > &point
 }
 
 
+template <int spacedim>
+const std::vector<Point<1> > &
+Manifold<spacedim>::
+get_line_support_points (const unsigned int n_intermediate_points) const
+{
+  if (points.size() <= n_intermediate_points ||
+      points[n_intermediate_points].get() == 0)
+    {
+      Threads::Mutex::ScopedLock lock(mutex);
+      if (points.size() <= n_intermediate_points)
+        points.resize(n_intermediate_points+1);
+
+      // another thread might have created points in the meantime
+      if (points[n_intermediate_points].get() == 0)
+        {
+          std_cxx1x::shared_ptr<QGaussLobatto<1> >
+          quadrature (new QGaussLobatto<1>(n_intermediate_points+2));
+          points[n_intermediate_points] = quadrature;
+        }
+    }
+  return points[n_intermediate_points]->get_points();
+}
+
+
 /* -------------------------- FlatManifold --------------------- */
 
 
-- 
2.39.5