Restored projection to StraightBoundaries

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

diff --git a/deal.II/include/deal.II/grid/tria_boundary.h b/deal.II/include/deal.II/grid/tria_boundary.h
index 882304fa536a66ed31cb7abecb5f9760d9f6336d..9367c8f3e72724d743ab4ab47ea66adeea5993a2 100644 (file)
--- 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 d33e1ac01fda8edad0b10dfe33f885404171ed9c..9bfca6e8be96af7aa0801a0b3951a0f61c782988 100644 (file)
--- 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"