From: heltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Tue, 28 Jan 2014 14:42:29 +0000 (+0000)
Subject: Fixed two distorted cells tests.
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=180e842e48cdf38602a8fcf1989d85cdb0e049a7;p=dealii-svn.git

Fixed two distorted cells tests.

git-svn-id: https://svn.dealii.org/branches/branch_manifold_id@32333 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 955962b69c..48016a8286 100644
--- a/deal.II/include/deal.II/grid/manifold.h
+++ b/deal.II/include/deal.II/grid/manifold.h
@@ -83,6 +83,25 @@ public:
   Point<spacedim>
   get_new_point(const std::vector<Point<spacedim> > &surrounding_points,
 		const std::vector<double> &weights) const = 0;
+  /**
+   * Given a point which lies close to the given manifold, it modifies
+   * it and projects it to manifold itself. This function is used in
+   * some mesh smoothing algorithms that try to move around points in
+   * order to improve the mesh quality but need to ensure that points
+   * that were on the manfold remain on the manifold.
+   *
+   * Derived classes do not need to implement this function unless
+   * mesh smoothing algorithms are used with a particular boundary
+   * object. The default implementation of this function calls the
+   * method with the same name, but without surrounding_points.
+   *
+   * This version of the projection function can be overwritten when
+   * only an approximation of the manifold is available, and this
+   * approximation requires surrounding points to be defined. 
+   */
+  virtual
+  Point<spacedim> project_to_manifold (const std::vector<Point<spacedim> > &surrounding_points,
+				       const Point<spacedim> &candidate) const;
 
   /**
    * Given a point which lies close to the given manifold, it modifies
diff --git a/deal.II/source/grid/grid_tools.cc b/deal.II/source/grid/grid_tools.cc
index e83cc14f5f..b385bb51b2 100644
--- a/deal.II/source/grid/grid_tools.cc
+++ b/deal.II/source/grid/grid_tools.cc
@@ -1751,6 +1751,9 @@ next_cell:
 
         const unsigned int structdim = Iterator::AccessorType::structure_dimension;
         const unsigned int spacedim  = Iterator::AccessorType::space_dimension;
+	std::vector<Point<spacedim> > vertices(GeometryInfo<structdim>::vertices_per_cell);
+	for(unsigned int i=0; i<vertices.size(); ++i)
+	  vertices[i] = object->vertex(i);
 
         // right now we can only deal
         // with cells that have been
@@ -1816,10 +1819,12 @@ next_cell:
                 else
                    gradient[d]
                      = ((objective_function (object,
-                                             manifold->project_to_manifold(object_mid_point + h))
+                                             manifold->project_to_manifold(vertices, 
+									   object_mid_point + h))
                          -
                          objective_function (object,
-                                             manifold->project_to_manifold(object_mid_point - h)))
+                                             manifold->project_to_manifold(vertices, 
+									   object_mid_point - h)))
                         /
                         eps);
               }
@@ -1853,7 +1858,8 @@ next_cell:
                                 gradient;
 
             if (respect_manifold == true)
-              object_mid_point = manifold->project_to_manifold(object_mid_point);
+              object_mid_point = manifold->project_to_manifold(vertices, 
+							       object_mid_point);
 
             // compute current value of the
             // objective function
diff --git a/deal.II/source/grid/manifold.cc b/deal.II/source/grid/manifold.cc
index 60f6fb8f25..325430448f 100644
--- a/deal.II/source/grid/manifold.cc
+++ b/deal.II/source/grid/manifold.cc
@@ -40,6 +40,15 @@ Manifold<spacedim>::project_to_manifold (const Point<spacedim> &candidate) const
   return candidate;
 }
 
+template <int spacedim>
+Point<spacedim>
+Manifold<spacedim>::project_to_manifold (const std::vector<Point<spacedim> > &,
+					 const Point<spacedim> &candidate) const 
+{
+  return project_to_manifold(candidate);
+}
+
+
 template <int spacedim>
 void
 Manifold<spacedim>::get_normals_at_points(const std::vector<Point<spacedim> > &points,
diff --git a/tests/bits/distorted_cells_03.cc b/tests/bits/distorted_cells_03.cc
index 4d32febe71..4ab1d099e2 100644
--- a/tests/bits/distorted_cells_03.cc
+++ b/tests/bits/distorted_cells_03.cc
@@ -23,6 +23,7 @@
 // distorted cells resulting from refinement
 
 #include "../tests.h"
+#include "my_boundary.h"
 #include <deal.II/base/logstream.h>
 #include <deal.II/base/quadrature_lib.h>
 #include <deal.II/grid/tria.h>
@@ -38,44 +39,6 @@
 
 #include <fstream>
 
-
-template <int dim>
-class MyBoundary : public Manifold<dim>
-{
-  virtual Point<dim>
-  get_new_point (const std::vector<Point<dim> > &v,
-		 const std::vector<double> &ws) const
-  {
-    switch(v.size()) {
-    case 2:
-      deallog << "Finding point between "
-	      << v[0] << " and "
-	      << v[1] << std::endl;
-
-      // in 2d, find a point that
-      // lies on the opposite side
-      // of the quad. in 3d, choose
-      // the midpoint of the edge
-      if (dim == 2)
-	return Point<dim>(0,0.75);
-      else
-	return (v[0] + v[1]) / 2;
-      break;
-    default:
-      deallog << "Finding point between "
-	      << v[0] << " and "
-	      << v[1] << " and "
-	      << v[2] << " and "
-	      << v[3] << std::endl;
-
-      return Point<dim>(0,0,.75);
-      break;
-    }
-  }
-};
-
-
-
 template <int dim>
 void check ()
 {
diff --git a/tests/bits/distorted_cells_04.cc b/tests/bits/distorted_cells_04.cc
index 715cbd5759..1810e3cde2 100644
--- a/tests/bits/distorted_cells_04.cc
+++ b/tests/bits/distorted_cells_04.cc
@@ -34,78 +34,7 @@
 #include <deal.II/fe/fe_values.h>
 
 #include <fstream>
-
-
-template <int dim>
-class MyBoundary : public Boundary<dim>
-{
-  virtual Point<dim>
-  get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const
-  {
-    deallog << "Finding point between "
-            << line->vertex(0) << " and "
-            << line->vertex(1) << std::endl;
-
-    // in 2d, find a point that
-    // lies on the opposite side
-    // of the quad. in 3d, choose
-    // the midpoint of the edge
-    if (dim == 2)
-      return Point<dim>(0,0.75);
-    else
-      return (line->vertex(0) + line->vertex(1)) / 2;
-  }
-
-  virtual Point<dim>
-  get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const
-  {
-    deallog << "Finding point between "
-            << quad->vertex(0) << " and "
-            << quad->vertex(1) << " and "
-            << quad->vertex(2) << " and "
-            << quad->vertex(3) << std::endl;
-
-    return Point<dim>(0,0,.75);
-  }
-
-  virtual
-  Point<dim>
-  project_to_surface (const typename Triangulation<dim>::line_iterator &,
-                      const Point<dim> &p) const
-  {
-    deallog << "Projecting line point " << p << std::endl;
-
-    if (dim == 2)
-      return Point<dim>(p[0], 0.75-1.75*std::fabs(p[0]));
-    else
-      return p;
-  }
-
-  virtual
-  Point<dim>
-  project_to_surface (const typename Triangulation<dim>::quad_iterator &,
-                      const Point<dim> &p) const
-  {
-    deallog << "Projecting quad point " << p << std::endl;
-
-    Assert (dim == 3, ExcInternalError());
-
-    return Point<dim>(p[0], p[1],
-                      0.75-1.75*std::max(std::fabs(p[0]),
-                                         std::fabs(p[1])));
-  }
-
-  virtual
-  Point<dim>
-  project_to_surface (const typename Triangulation<dim>::hex_iterator &,
-                      const Point<dim> &) const
-  {
-    Assert (false, ExcInternalError());
-    return Point<dim>();
-  }
-};
-
-
+#include "my_boundary.h"
 
 template <int dim>
 void check ()
diff --git a/tests/bits/my_boundary.h b/tests/bits/my_boundary.h
new file mode 100644
index 0000000000..a251d30406
--- /dev/null
+++ b/tests/bits/my_boundary.h
@@ -0,0 +1,95 @@
+// ---------------------------------------------------------------------
+// $Id: distorted_cells_04.cc 32209 2014-01-15 00:25:27Z heltai $
+//
+// Copyright (C) 2003 - 2013 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// like _03, but catch the exception and pass it to GridTools::fix_up_distorted_child_cells
+
+#include "../tests.h"
+#include <deal.II/grid/manifold.h>
+
+template <int dim>
+class MyBoundary : public Manifold<dim>
+{
+  virtual Point<dim>
+  get_new_point (const std::vector<Point<dim> > &v,
+		 const std::vector<double> &ws) const
+  {
+    switch(v.size()) {
+    case 2:
+      deallog << "Finding point between "
+	      << v[0] << " and "
+	      << v[1] << std::endl;
+
+      // in 2d, find a point that
+      // lies on the opposite side
+      // of the quad. in 3d, choose
+      // the midpoint of the edge
+      if (dim == 2)
+	return Point<dim>(0,0.75);
+      else
+	return (v[0] + v[1]) / 2;
+      break;
+    default:
+      deallog << "Finding point between "
+	      << v[0] << " and "
+	      << v[1] << " and "
+	      << v[2] << " and "
+	      << v[3] << std::endl;
+
+      return Point<dim>(0,0,.75);
+      break;
+    }
+  }
+
+  virtual
+  Point<dim>
+  project_to_manifold (const std::vector<Point<dim> > &v,
+		       const Point<dim> &p) const
+  {
+    switch(v.size()) {
+    case 2:
+      {
+	deallog << "Projecting line point " << p << std::endl;
+	if (dim == 2)
+	  return Point<dim>(p[0], 0.75-1.75*std::fabs(p[0]));
+	else
+	  return p;
+      }
+      break;
+    case 4:
+      {
+	deallog << "Projecting quad point " << p << std::endl;
+
+	Assert (dim == 3, ExcInternalError());
+
+	return Point<dim>(p[0], p[1],
+			  0.75-1.75*std::max(std::fabs(p[0]),
+					     std::fabs(p[1])));
+      }
+      break;
+    default:
+      {
+	Assert (false, ExcInternalError());
+	return Point<dim>();
+      }
+    }
+  }
+};
+
+
+
+