Made intermediate points of manifold lines and quads equal to new implementation...

git-svn-id: https://svn.dealii.org/branches/branch_manifold_id@32163 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/include/deal.II/grid/manifold.h b/deal.II/include/deal.II/grid/manifold.h
index 553d9bb0ebe072f54926bf34068867c0ba88612f..d088f0e813a50de0ca25192d3d0cbbe8000cf626 100644 (file)
--- 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 7a6de0639f7225c25484e1ef8d4307a299ae98ff..fb02d1bb2098e19100658bf4eb656a63c354ada5 100644 (file)
--- 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 --------------------- */