]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Restored approximation of the normal at arbitrary points for StraightBoundary
authorheltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 18 Jan 2014 17:08:50 +0000 (17:08 +0000)
committerheltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 18 Jan 2014 17:08:50 +0000 (17:08 +0000)
git-svn-id: https://svn.dealii.org/branches/branch_manifold_id@32248 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/include/deal.II/grid/manifold.h
deal.II/include/deal.II/grid/tria_boundary.h
deal.II/include/deal.II/numerics/vector_tools.templates.h
deal.II/source/grid/manifold.cc
deal.II/source/grid/tria_boundary.cc

index e26dd7e4e5a7c198b9310373722fecb746baca2d..ed85ca36786f4c879e96326b94bd7823c4901259 100644 (file)
@@ -103,14 +103,28 @@ public:
     * Given a vector of points, return the normals to the Manifold in
     * those points. Notice that the Manifold may or may not be a
     * codimension one manifold. If it is not, this function will throw
-    * an Error. Input arguments must be of the same size. The default
-    * implementation calls for each point the function
+    * an Error. @points and @normals must be of the same size. The
+    * default implementation calls for each point the function
     * normal_vector. Derived classes can overload that function.
     */
     virtual
     void get_normals_at_points(const std::vector<Point<spacedim> > &points,
                               std::vector<Point<spacedim> > &normals) const;
 
+    /**
+     * Compute the normal to the manifold at the given point. The
+     * default implementation of this function just calls
+     * normal_vector. The difference from this function and that one,
+     * is that one can use this function when only an approximation of
+     * the normal is required, for example when computing normals at
+     * quadrature points. Note that this function only makes sense
+     * when the manifold is a codimension one manifold.
+     */
+    virtual
+    Point<spacedim>
+    normal_vector(const std::vector<Point<spacedim> > &vertices,
+                 const Point<spacedim> &point) const;
+
     /**
      * Compute the normal to the manifold at the given point. The
      * default implementation of this function throws an
index 6c1fc8ee08a6437e19b418ffbe9604e3b4703a71..72b31d4db457cdcfd47747bb8384a221e4093aa1 100644 (file)
@@ -155,11 +155,39 @@ public:
     virtual
     void get_normals_at_points(const std::vector<Point<spacedim> > &points,
                               std::vector<Point<spacedim> > &normals) const;
-    
-    
+  /**
+    * Given a point and a vector of vertices, compute an approximation
+    * of the normal to the manifold.
+    */
+    virtual
+    Point<spacedim> 
+    normal_vector(const std::vector<Point<spacedim> > &vertices,
+                 const Point<spacedim> &point) const;
 };
 
 
+/* -------------- declaration of explicit specializations ------------- */
+
+#ifndef DOXYGEN
+
+template <>
+Point<1>
+StraightBoundary<1,1>::normal_vector(const std::vector<Point<1> > &vertices,
+                                    const Point<1> &point) const;
+
+template <>
+Point<2>
+StraightBoundary<1,2>::normal_vector(const std::vector<Point<2> > &vertices,
+                                    const Point<2> &point) const;
+
+template <>
+Point<3>
+StraightBoundary<1,3>::normal_vector(const std::vector<Point<3> > &vertices,
+                                    const Point<3> &point) const;
+
+
+#endif
+
 DEAL_II_NAMESPACE_CLOSE
 
 #endif
index 8d2ed0861c762771d79f5f4ee4bb9272e8d7f6ca..cc1b81c42e0c7dfd5d3241e4095c5524dcd60438 100644 (file)
@@ -4419,9 +4419,12 @@ namespace VectorTools
                     // sign of the normal vector provided by the boundary
                     // if they should point in different directions. this is the
                     // case in tests/deal.II/no_flux_11.
+                   std::vector<Point<dim> > vertices(GeometryInfo<dim>::vertices_per_face);
+                   for(unsigned int v=0; v<GeometryInfo<dim>::vertices_per_face; ++v)
+                     vertices[v] = cell->face(face_no)->vertex(v);
                     Point<dim> normal_vector
                       = (cell->face(face_no)->get_boundary()
-                         .normal_vector (fe_values.quadrature_point(i)));
+                         .normal_vector (vertices, fe_values.quadrature_point(i)));
                     if (normal_vector * fe_values.normal_vector(i) < 0)
                       normal_vector *= -1;
                     Assert (std::fabs(normal_vector.norm() - 1) < 1e-14,
index ffd6b6dd054d2992c4456a250cba9705a4afd25b..1d598217a8bea86b0d87326cc905481b3aa8477c 100644 (file)
@@ -50,6 +50,14 @@ Manifold<spacedim>::get_normals_at_points(const std::vector<Point<spacedim> > &p
     normals[i] = normal_vector(points[i]);
 }
 
+template <int spacedim>
+Point<spacedim>
+Manifold<spacedim>::normal_vector(const std::vector<Point<spacedim> > &,
+                                 const Point<spacedim> &point) const
+{
+  return normal_vector(point);
+}
+
 template <int spacedim>
 Point<spacedim>
 Manifold<spacedim>::normal_vector(const Point<spacedim> &point) const
index e220300f6226e011c5082c59febd6a4ca3ee7645..d33e1ac01fda8edad0b10dfe33f885404171ed9c 100644 (file)
@@ -16,6 +16,8 @@
   
 #include <deal.II/base/tensor.h>
 #include <deal.II/grid/tria_boundary.h>
+#include <deal.II/fe/fe_q.h>
+
 DEAL_II_NAMESPACE_OPEN
 
 /* -------------------------- Boundary --------------------- */
@@ -28,6 +30,45 @@ Boundary<dim, spacedim>::~Boundary ()
 /* -------------------------- StraightBoundary --------------------- */
 
 
+namespace internal
+{
+  namespace
+  {
+    /**                                                                                             
+     * Compute the normalized cross product of a set of dim-1 basis                                 
+     * vectors.                                                                                     
+     */
+    Tensor<1,2>
+    normalized_alternating_product (const Tensor<1,2> (&basis_vectors)[1])
+    {
+      Tensor<1,2> tmp;
+      cross_product (tmp, basis_vectors[0]);
+      return tmp/tmp.norm();
+    }
+
+
+    Tensor<1,3>
+    normalized_alternating_product (const Tensor<1,3> ( &)[1])
+    {
+      // we get here from StraightBoundary<2,3>::normal_vector, but                                 
+      // the implementation below is bogus for this case anyway                                     
+      // (see the assert at the beginning of that function).                                        
+      Assert (false, ExcNotImplemented());
+      return Tensor<1,3>();
+    }
+
+
+    Tensor<1,3>
+    normalized_alternating_product (const Tensor<1,3> (&basis_vectors)[2])    {
+      Tensor<1,3> tmp;
+      cross_product (tmp, basis_vectors[0], basis_vectors[1]);
+      return tmp/tmp.norm();
+    }
+  }
+}
+
+
+
 template <int dim, int spacedim>
 StraightBoundary<dim, spacedim>::StraightBoundary ()
 {}
@@ -91,6 +132,135 @@ StraightBoundary<dim, spacedim>::get_normals_at_points(const std::vector<Point<s
 }
 
 
+template <>
+Point<1>
+StraightBoundary<1,1>::
+normal_vector (const std::vector<Point<1> > &vertices, 
+               const Point<1> &) const
+{
+  Assert (false, ExcNotImplemented());
+  return Point<1>();
+}
+
+
+
+template <>
+Point<2>
+StraightBoundary<1,2>::
+normal_vector (const std::vector<Point<2> > &vertices, 
+               const Point<2> &) const
+{
+  Assert (false, ExcNotImplemented());
+  return Point<2>();
+}
+
+
+
+template <>
+Point<3>
+StraightBoundary<1,3>::
+normal_vector (const std::vector<Point<3> > &vertices, 
+               const Point<3> &) const
+{
+  Assert (false, ExcNotImplemented());
+  return Point<3>();
+}
+
+
+template <int dim, int spacedim>
+Point<spacedim>
+StraightBoundary<dim, spacedim>::normal_vector(const std::vector<Point<spacedim> > &vertices, 
+                                              const Point<spacedim> &p) const {
+  AssertDimension(vertices.size(), GeometryInfo<dim>::vertices_per_face);
+  
+  // I don't think the implementation below will work when dim!=spacedim;                           
+  // in fact, I believe that we don't even have enough information here,                            
+  // because we would need to know not only about the tangent vectors                               
+  // of the face, but also of the cell, to compute the normal vector.                               
+  // Someone will have to think about this some more.                                               
+  Assert (dim == spacedim, ExcNotImplemented());
+
+  // in order to find out what the normal vector is, we first need to                               
+  // find the reference coordinates of the point p on the given face,                               
+  // or at least the reference coordinates of the closest point on the                              
+  // face                                                                                           
+  //                                                                                                
+  // in other words, we need to find a point xi so that f(xi)=||F(xi)-p||^2->min                    
+  // where F(xi) is the mapping. this algorithm is implemented in                                   
+  // MappingQ1<dim,spacedim>::transform_real_to_unit_cell but only for cells,                       
+  // while we need it for faces here. it's also implemented in somewhat                             
+  // more generality there using the machinery of the MappingQ1 class          
+  // while we really only need it for a specific case here                                          
+  //                                                                                                
+  // in any case, the iteration we use here is a Gauss-Newton's iteration with                      
+  //   xi^{n+1} = xi^n - H(xi^n)^{-1} J(xi^n)                                                       
+  // where                                                                                          
+  //   J(xi) = (grad F(xi))^T (F(xi)-p)                                                             
+  // and                                                                                            
+  //   H(xi) = [grad F(xi)]^T [grad F(xi)]                                                          
+  // In all this,                                                                                   
+  //   F(xi) = sum_v vertex[v] phi_v(xi)     
+  //                                                                                                
+  // in any case, the iteration we use here is a Gauss-Newton's iteration with                      
+  //   xi^{n+1} = xi^n - H(xi^n)^{-1} J(xi^n)                                                       
+  // where                                                                                          
+  //   J(xi) = (grad F(xi))^T (F(xi)-p)                                                             
+  // and                                                                                            
+  //   H(xi) = [grad F(xi)]^T [grad F(xi)]                                                          
+  // In all this,                                                                                   
+  //   F(xi) = sum_v vertex[v] phi_v(xi)                                                            
+  // We get the shape functions phi_v from an object of type FE_Q<dim-1>(1)                         
+
+  // we start with the point xi=1/2, xi=(1/2,1/2), ...                                              
+  const unsigned int facedim = dim-1;
+
+  Point<facedim> xi;
+  for (unsigned int i=0; i<facedim; ++i)
+    xi[i] = 1./2;
+  
+  FE_Q<facedim> linear_fe(1);
+
+  const double eps = 1e-12;
+  Tensor<1,spacedim> grad_F[facedim];
+  while (true)
+    {
+      Point<spacedim> F;
+      for (unsigned int v=0; v<GeometryInfo<facedim>::vertices_per_cell; ++v)
+        F += vertices[v] * linear_fe.shape_value(v, xi);
+
+      for (unsigned int i=0; i<facedim; ++i)
+        {
+          grad_F[i] = 0;
+          for (unsigned int v=0; v<GeometryInfo<facedim>::vertices_per_cell; ++v)
+            grad_F[i] += vertices[v] * linear_fe.shape_grad(v, xi)[i];
+        }
+
+      Tensor<1,facedim> J;
+      for (unsigned int i=0; i<facedim; ++i)
+        for (unsigned int j=0; j<spacedim; ++j)
+          J[i] += grad_F[i][j] * (F-p)[j];
+
+      Tensor<2,facedim> H;
+      for (unsigned int i=0; i<facedim; ++i)
+       for (unsigned int j=0; j<facedim; ++j)
+          for (unsigned int k=0; k<spacedim; ++k)
+            H[i][j] += grad_F[i][k] * grad_F[j][k];
+
+      const Point<facedim> delta_xi = -invert(H) * J;
+      xi += delta_xi;
+
+      if (delta_xi.norm() < eps)
+        break;
+    }
+
+  // so now we have the reference coordinates xi of the point p.                                    
+  // we then have to compute the normal vector, which we can do                                     
+  // by taking the (normalize) alternating product of all the tangent                               
+  // vectors given by grad_F                                                                        
+  return internal::normalized_alternating_product(grad_F);
+}
+
+
 // explicit instantiations
 #include "tria_boundary.inst"
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.