From 3bb45d7bd2470386597ed82c87c3d5447dd6f6cc Mon Sep 17 00:00:00 2001
From: heltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Wed, 29 Jan 2014 11:07:50 +0000
Subject: [PATCH] Restored projection to StraightBoundaries

git-svn-id: https://svn.dealii.org/branches/branch_manifold_id@32342 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/include/deal.II/grid/tria_boundary.h |  9 ++
 deal.II/source/grid/tria_boundary.cc         | 95 ++++++++++++++++++++
 2 files changed, 104 insertions(+)

diff --git a/deal.II/include/deal.II/grid/tria_boundary.h b/deal.II/include/deal.II/grid/tria_boundary.h
index 882304fa53..9367c8f3e7 100644
--- a/deal.II/include/deal.II/grid/tria_boundary.h
+++ b/deal.II/include/deal.II/grid/tria_boundary.h
@@ -164,6 +164,15 @@ public:
   Point<spacedim> 
   normal_vector(const std::vector<Point<spacedim> > &vertices,
 		const Point<spacedim> &point) const;
+  /**
+   * Given a point and a vector of vertices, project the point
+   * back to the flat boundary. 
+   */
+  virtual
+  Point<spacedim> 
+  project_to_manifold(const std::vector<Point<spacedim> > &vertices,
+		      const Point<spacedim> &point) const;
+
 };
 
 
diff --git a/deal.II/source/grid/tria_boundary.cc b/deal.II/source/grid/tria_boundary.cc
index d33e1ac01f..9bfca6e8be 100644
--- a/deal.II/source/grid/tria_boundary.cc
+++ b/deal.II/source/grid/tria_boundary.cc
@@ -260,6 +260,101 @@ StraightBoundary<dim, spacedim>::normal_vector(const std::vector<Point<spacedim>
   return internal::normalized_alternating_product(grad_F);
 }
 
+template <int dim, int spacedim>
+Point<spacedim>
+StraightBoundary<dim, spacedim>::project_to_manifold(const std::vector<Point<spacedim> > &v, 
+						     const Point<spacedim> &trial_point) const {
+  if(spacedim <= 1)
+    return trial_point;
+  else
+    switch(v.size()) {
+    case 2:
+      {
+	// find the point that lies on the line p1--p2. the formulas
+	// pan out to something rather simple because the mapping to
+	// the line is linear
+	const double s = (trial_point-v[0])*(v[1]-v[0]) / ((v[1]-v[0])*(v[1]-v[0]));
+	return v[0] + s*(v[1]-v[0]);
+      }
+      break;
+    case 4:
+      {
+	// let's look at this for simplicity for a quad (dim==2) in a
+	// space with spacedim>2:
+
+	// all points on the surface are given by x(\xi) = sum_i v_i
+	//   phi_x(\xi) where v_i are the vertices of the quad, and
+	//   \xi=(\xi_1,\xi_2) are the reference coordinates of the
+	//   quad. so what we are trying to do is find a point x on
+	//   the surface that is closest to the point y. there are
+	//   different ways to solve this problem, but in the end it's
+	//   a nonlinear problem and we have to find reference
+	//   coordinates \xi so that J(\xi) = 1/2 || x(\xi)-y ||^2 is
+	//   minimal. x(\xi) is a function that is dim-linear in \xi,
+	//   so J(\xi) is a polynomial of degree 2*dim that we'd like
+	//   to minimize. unless dim==1, we'll have to use a Newton
+	//   method to find the answer. This leads to the following
+	//   formulation of Newton steps:
+	//
+	// Given \xi_k, find \delta\xi_k so that H_k \delta\xi_k = -
+	//   F_k where H_k is an approximation to the second
+	//   derivatives of J at \xi_k, and F_k is the first
+	//   derivative of J.  We'll iterate this a number of times
+	//   until the right hand side is small enough. As a stopping
+	//   criterion, we terminate if ||\delta\xi||<eps.
+	//
+	// As for the Hessian, the best choice would be H_k =
+	// J''(\xi_k) but we'll opt for the simpler Gauss-Newton form
+	// H_k = A^T A i.e.  (H_k)_{nm} = \sum_{i,j} v_i*v_j *
+	// \partial_n phi_i * \partial_m phi_j we start at
+	// xi=(0.5,0.5).
+	Point<dim> xi;
+	for (unsigned int d=0; d<dim; ++d)
+	  xi[d] = 0.5;
+
+	Point<spacedim> x_k;
+	for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+	  x_k += v[i] * GeometryInfo<dim>::d_linear_shape_function (xi, i);
+
+	do
+	  {
+	    Tensor<1,dim> F_k;
+	    for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+	      F_k += (x_k-trial_point)*v[i] *
+		GeometryInfo<dim>::d_linear_shape_function_gradient (xi, i);
+
+	    Tensor<2,dim> H_k;
+	    for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+	      for (unsigned int j=0; j<GeometryInfo<dim>::vertices_per_cell; ++j)
+		{
+		  Tensor<2,dim> tmp;
+		  outer_product (tmp,
+				 GeometryInfo<dim>::d_linear_shape_function_gradient (xi, i),
+				 GeometryInfo<dim>::d_linear_shape_function_gradient (xi, j));
+		  H_k += (v[i] * v[j]) * tmp;
+		}
+
+	    const Point<dim> delta_xi = - invert(H_k) * F_k;
+	    xi += delta_xi;
+
+	    x_k = Point<spacedim>();
+	    for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+	      x_k += v[i] * GeometryInfo<dim>::d_linear_shape_function (xi, i);
+
+	    if (delta_xi.norm() < 1e-5)
+	      break;
+	  }
+	while (true);
+
+	return x_k;
+      }
+      break;
+    default:
+      Assert(false, ExcNotImplemented());
+      return trial_point;
+    }
+}
+
 
 // explicit instantiations
 #include "tria_boundary.inst"
-- 
2.39.5