From: Wolfgang Bangerth <bangerth@math.tamu.edu>
Date: Fri, 23 Mar 2001 14:48:42 +0000 (+0000)
Subject: Boundary<dim>::get_tangents_at_vertices to generate a description of the tangent... 
X-Git-Tag: v8.0.0~19502
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=df29b5725891c72c1084228fd96313757218c7cd;p=dealii.git

Boundary<dim>::get_tangents_at_vertices to generate a description of the tangent space to a boundary object at each vertex of a face.


git-svn-id: https://svn.dealii.org/trunk@4278 0785d39b-7218-0410-832d-ea1e28bc413d
---

diff --git a/deal.II/deal.II/include/grid/tria_boundary.h b/deal.II/deal.II/include/grid/tria_boundary.h
index 650bb9a680..9eb1af8c9e 100644
--- a/deal.II/deal.II/include/grid/tria_boundary.h
+++ b/deal.II/deal.II/include/grid/tria_boundary.h
@@ -76,6 +76,48 @@ template <int dim>
 class Boundary : public Subscriptor
 {
   public:
+
+				     /**
+				      * Structure keeping information
+				      * about the tangents at the
+				      * vertices of a face of a
+				      * cell. Thus, there are
+				      * @p{GeometryInfo<dim>::vertices_per_face}
+				      * data sets of @p{dim-1} tangent
+				      * vectors each, that define the
+				      * tangent space of the boundary
+				      * at a given vertex. Note that
+				      * the vectors stored in this
+				      * object does not require that
+				      * the tangent vectors actually
+				      * form an orthogonal basis (as
+				      * long as they form a basis,
+				      * i.e. are not degenerate), not
+				      * are they required to be normed
+				      * somehow.
+				      *
+				      * For obvious reasons, this
+				      * class is not useful in 1d.
+				      */
+    struct FaceVertexTangents 
+    {
+					 /**
+					  * Number of tangent vectors
+					  * that define the tangential
+					  * space. The value is equal
+					  * to @p{dim-1}, with a
+					  * non-sensical value for 1d.
+					  */
+	static const unsigned int n_tangents = (dim>1 ? dim-1 : static_cast<unsigned int>(-1));
+
+					 /**
+					  * Array of tangent vectors,
+					  * as described above.
+					  */
+	Tensor<1,dim> tangents[GeometryInfo<dim>::vertices_per_face][n_tangents];
+    };
+    
+	
 				     /**
 				      * Destructor. Does nothing here, but
 				      * needs to be declared to make it
@@ -182,6 +224,43 @@ class Boundary : public Subscriptor
     get_intermediate_points_on_quad (const typename Triangulation<dim>::quad_iterator &quad,
 				     typename std::vector<Point<dim> > &points) const;
 
+				     /**
+				      * Compute a basis of the tangent
+				      * space at each vertex of the
+				      * given face. It is not
+				      * necessary to compute
+				      * derivatives in any particular
+				      * direction at each vertex, as
+				      * long as the returned basis
+				      * actually spans a space of
+				      * dimension @p{dim-1}. Also, it
+				      * is not required that the
+				      * tangent vectors be normed
+				      * somehow.
+				      *
+				      * This function is
+				      * needed to compute data for C1
+				      * mappings. The default
+				      * implementation is to throw an
+				      * error, so you need not
+				      * overload this function in case
+				      * you don not intend to use C1
+				      * mappings.
+				      *
+				      * Note that when computing
+				      * tangents at a vertex where the
+				      * boundary is not
+				      * differentiable, you have to
+				      * make sure that you compute the
+				      * one-sided derivatives,
+				      * i.e. derivatives with respect
+				      * to the given face.
+				      */
+    virtual void
+    get_tangents_at_vertices (const typename Triangulation<dim>::face_iterator &face,
+			      FaceVertexTangents &face_vertex_tangents) const;
+      
+			     
 				     /**
 				      * Exception.
 				      */
@@ -210,7 +289,8 @@ class Boundary : public Subscriptor
  *   @author Wolfgang Bangerth, 1998, Ralf Hartmann, 2001
  */
 template <int dim>
-class StraightBoundary : public Boundary<dim> {
+class StraightBoundary : public Boundary<dim>
+{
   public:
 				     /**
 				      * Let the new point be the arithmetic
@@ -269,6 +349,19 @@ class StraightBoundary : public Boundary<dim> {
     get_intermediate_points_on_quad (const typename Triangulation<dim>::quad_iterator &quad,
 				     typename std::vector<Point<dim> > &points) const;
 
+				     /**
+				      * Compute a basis of the tangent
+				      * space at each vertex of the
+				      * given face.
+				      *
+				      * Refer to the general
+				      * documentation of this class
+				      * and the documentation of the
+				      * base class.
+				      */
+    virtual void
+    get_tangents_at_vertices (const typename Triangulation<dim>::face_iterator &face,
+			      typename Boundary<dim>::FaceVertexTangents &face_vertex_tangents) const;
 };
 
 
diff --git a/deal.II/deal.II/include/grid/tria_boundary_lib.h b/deal.II/deal.II/include/grid/tria_boundary_lib.h
index 3ccb5d45fa..a972d61bd3 100644
--- a/deal.II/deal.II/include/grid/tria_boundary_lib.h
+++ b/deal.II/deal.II/include/grid/tria_boundary_lib.h
@@ -84,6 +84,20 @@ class HyperBallBoundary : public StraightBoundary<dim>
     get_intermediate_points_on_quad (const typename Triangulation<dim>::quad_iterator &quad,
 				     typename std::vector<Point<dim> > &points) const;
 
+				     /**
+				      * Compute a basis of the tangent
+				      * space at each vertex of the
+				      * given face.
+				      *
+				      * Refer to the general
+				      * documentation of this class
+				      * and the documentation of the
+				      * base class.
+				      */
+    virtual void
+    get_tangents_at_vertices (const typename Triangulation<dim>::face_iterator &face,
+			      typename Boundary<dim>::FaceVertexTangents &face_vertex_tangents) const;
+
 				     /**
 				      * Return the center of the ball.
 				      */
@@ -186,6 +200,21 @@ class HalfHyperBallBoundary : public HyperBallBoundary<dim>
     virtual void
     get_intermediate_points_on_quad (const typename Triangulation<dim>::quad_iterator &quad,
 				     typename std::vector<Point<dim> > &points) const;
+
+
+				     /**
+				      * Compute a basis of the tangent
+				      * space at each vertex of the
+				      * given face.
+				      *
+				      * Refer to the general
+				      * documentation of this class
+				      * and the documentation of the
+				      * base class.
+				      */
+    virtual void
+    get_tangents_at_vertices (const typename Triangulation<dim>::face_iterator &face,
+			      typename Boundary<dim>::FaceVertexTangents &face_vertex_tangents) const;
 };
 
 
diff --git a/deal.II/deal.II/source/grid/tria_boundary.cc b/deal.II/deal.II/source/grid/tria_boundary.cc
index 5ef3b01668..85c8e67fc9 100644
--- a/deal.II/deal.II/source/grid/tria_boundary.cc
+++ b/deal.II/deal.II/source/grid/tria_boundary.cc
@@ -19,11 +19,16 @@
 #include <cmath>
 
 
+
+/* -------------------------- StraightBoundary --------------------- */
+
+
 template <int dim>
 Boundary<dim>::~Boundary ()
 {};
 
 
+
 template <int dim>
 Point<dim>
 Boundary<dim>::get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &) const 
@@ -33,26 +38,43 @@ Boundary<dim>::get_new_point_on_quad (const typename Triangulation<dim>::quad_it
 };
 
 
+
+template <int dim>
+void
+Boundary<dim>::
+get_intermediate_points_on_line (const typename Triangulation<dim>::line_iterator &,
+				 std::vector<Point<dim> > &) const
+{
+  Assert (false, ExcPureVirtualFunctionCalled());
+};
+
+
+
 template <int dim>
 void
-Boundary<dim>::get_intermediate_points_on_line (
-  const typename Triangulation<dim>::line_iterator &,
-  std::vector<Point<dim> > &) const
+Boundary<dim>::
+get_intermediate_points_on_quad (const typename Triangulation<dim>::quad_iterator &,
+				 std::vector<Point<dim> > &) const
 {
   Assert (false, ExcPureVirtualFunctionCalled());
 };
 
 
+
 template <int dim>
 void
-Boundary<dim>::get_intermediate_points_on_quad (
-  const typename Triangulation<dim>::quad_iterator &,
-  std::vector<Point<dim> > &) const
+Boundary<dim>::
+get_tangents_at_vertices (const typename Triangulation<dim>::face_iterator &,
+			  FaceVertexTangents                               &) const
 {
   Assert (false, ExcPureVirtualFunctionCalled());
 };
 
 
+
+/* -------------------------- StraightBoundary --------------------- */
+
+
 template <int dim>
 Point<dim>
 StraightBoundary<dim>::get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const 
@@ -77,7 +99,8 @@ StraightBoundary<dim>::get_new_point_on_quad (const typename Triangulation<dim>:
 
 template <int dim>
 Point<dim>
-StraightBoundary<dim>::get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const 
+StraightBoundary<dim>::
+get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const 
 {
   return (quad->vertex(0) + quad->vertex(1) +
 	  quad->vertex(2) + quad->vertex(3) +
@@ -94,9 +117,9 @@ StraightBoundary<dim>::get_new_point_on_quad (const typename Triangulation<dim>:
 
 template <int dim>
 void
-StraightBoundary<dim>::get_intermediate_points_on_line (
-  const typename Triangulation<dim>::line_iterator &,
-  typename std::vector<Point<dim> > &) const
+StraightBoundary<dim>::
+get_intermediate_points_on_line (const typename Triangulation<dim>::line_iterator &,
+				 typename std::vector<Point<dim> > &) const
 {
   Assert(false, typename Boundary<dim>::ExcFunctionNotUseful(dim));
 }
@@ -104,11 +127,12 @@ StraightBoundary<dim>::get_intermediate_points_on_line (
 
 #else
 
+
 template <int dim>
 void
-StraightBoundary<dim>::get_intermediate_points_on_line (
-  const typename Triangulation<dim>::line_iterator &line,
-  typename std::vector<Point<dim> > &points) const
+StraightBoundary<dim>::
+get_intermediate_points_on_line (const typename Triangulation<dim>::line_iterator &line,
+				 typename std::vector<Point<dim> > &points) const
 {
   const unsigned int n=points.size();
   Assert(n>0, ExcInternalError());
@@ -131,9 +155,9 @@ StraightBoundary<dim>::get_intermediate_points_on_line (
 
 template <int dim>
 void
-StraightBoundary<dim>::get_intermediate_points_on_quad (
-  const typename Triangulation<dim>::quad_iterator &,
-  typename std::vector<Point<dim> > &) const
+StraightBoundary<dim>::
+get_intermediate_points_on_quad (const typename Triangulation<dim>::quad_iterator &,
+				 typename std::vector<Point<dim> > &) const
 {
   Assert(false, typename Boundary<dim>::ExcFunctionNotUseful(dim));
 }
@@ -142,9 +166,9 @@ StraightBoundary<dim>::get_intermediate_points_on_quad (
 
 template <int dim>
 void
-StraightBoundary<dim>::get_intermediate_points_on_quad (
-  const typename Triangulation<dim>::quad_iterator &quad,
-  typename std::vector<Point<dim> > &points) const
+StraightBoundary<dim>::
+get_intermediate_points_on_quad (const typename Triangulation<dim>::quad_iterator &quad,
+				 typename std::vector<Point<dim> > &points) const
 {
   const unsigned int n=points.size(),
 		     m=static_cast<unsigned int>(sqrt(n));
@@ -173,6 +197,63 @@ StraightBoundary<dim>::get_intermediate_points_on_quad (
 #endif
 
 
+
+#if deal_II_dimension == 1
+
+template <>
+void
+StraightBoundary<1>::
+get_tangents_at_vertices (const Triangulation<1>::face_iterator &,
+			  FaceVertexTangents &) const
+{
+  Assert (false, Boundary<1>::ExcFunctionNotUseful(1));
+};
+
+#endif
+
+
+#if deal_II_dimension == 2
+
+template <>
+void
+StraightBoundary<2>::
+get_tangents_at_vertices (const Triangulation<2>::face_iterator &face,
+			  Boundary<2>::FaceVertexTangents &face_vertex_tangents) const
+{
+  const unsigned int dim=2;
+  const Tensor<1,dim> tangent = face->vertex(1) - face->vertex(0);
+  for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_face; ++vertex)
+    face_vertex_tangents.tangents[vertex][0] = tangent;
+};
+
+#endif
+
+
+
+#if deal_II_dimension == 3
+
+template <>
+void
+StraightBoundary<3>::
+get_tangents_at_vertices (const Triangulation<3>::face_iterator &face,
+			  Boundary<3>::FaceVertexTangents &face_vertex_tangents) const
+{
+  face_vertex_tangents.tangents[0][0] = face->vertex(1)-face->vertex(0);
+  face_vertex_tangents.tangents[0][1] = face->vertex(3)-face->vertex(0);
+  
+  face_vertex_tangents.tangents[1][0] = face->vertex(0)-face->vertex(1);
+  face_vertex_tangents.tangents[1][1] = face->vertex(2)-face->vertex(1);
+  
+  face_vertex_tangents.tangents[2][0] = face->vertex(1)-face->vertex(2);
+  face_vertex_tangents.tangents[2][1] = face->vertex(3)-face->vertex(2);
+  
+  face_vertex_tangents.tangents[3][0] = face->vertex(0)-face->vertex(3);
+  face_vertex_tangents.tangents[3][1] = face->vertex(2)-face->vertex(3);  
+};
+
+#endif
+
+
 // explicit instantiations
 template class Boundary<deal_II_dimension>;
 template class StraightBoundary<deal_II_dimension>;
diff --git a/deal.II/deal.II/source/grid/tria_boundary_lib.cc b/deal.II/deal.II/source/grid/tria_boundary_lib.cc
index aa98b78406..447d2b4b96 100644
--- a/deal.II/deal.II/source/grid/tria_boundary_lib.cc
+++ b/deal.II/deal.II/source/grid/tria_boundary_lib.cc
@@ -213,6 +213,133 @@ HyperBallBoundary<dim>::get_intermediate_points_on_quad (
 }
 
 
+
+#if deal_II_dimension == 1
+
+template <>
+void
+HyperBallBoundary<1>::
+get_tangents_at_vertices (const Triangulation<1>::face_iterator &,
+			  FaceVertexTangents &) const
+{
+  Assert (false, Boundary<1>::ExcFunctionNotUseful(1));
+};
+
+#endif
+
+
+#if deal_II_dimension == 2
+
+template <>
+void
+HyperBallBoundary<2>::
+get_tangents_at_vertices (const Triangulation<2>::face_iterator &face,
+			  Boundary<2>::FaceVertexTangents &face_vertex_tangents) const
+{
+  const unsigned int dim=2;
+				   // construct a tangential vector at
+				   // each of the vertices
+  for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_face; ++vertex)
+    {
+      const Tensor<1,dim> radius = face->vertex(vertex)-center;
+      face_vertex_tangents.tangents[vertex][0] = Point<dim>(radius[1], -radius[0]);
+    };
+};
+
+#endif
+
+
+
+#if deal_II_dimension == 3
+
+template <>
+void
+HyperBallBoundary<3>::
+get_tangents_at_vertices (const Triangulation<3>::face_iterator &face,
+			  Boundary<3>::FaceVertexTangents &face_vertex_tangents) const
+{
+  const unsigned int dim=3;
+				   // construct two tangential vectors
+				   // at each of the vertices
+  for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_face; ++vertex)
+    {
+				       // first take the radial vector
+      const Tensor<1,dim> radius = face->vertex(vertex)-center;
+				       // next we need to construct
+				       // two vectors that are
+				       // perpendicular to the radial
+				       // vector. to do so, first get
+				       // _one_ vector that is
+				       // perpendicular to the radial
+				       // one, and construct the
+				       // second by cross product.
+				       // construct the first
+				       // perpendicular vector by
+				       // forming the cross product of
+				       // the radial vector with
+				       // _some_ vector that is not
+				       // parallel to the radial
+				       // vector
+				       //
+				       // the problem therefore
+				       // reduces to finding an
+				       // auxiliary vector that is not
+				       // linearly dependent to the
+				       // radial vector. this we do as
+				       // follows (note that we know
+				       // that the length of the
+				       // radius vector is non-zero):
+				       //
+				       // 1. find the smallest
+				       // component of r by modulus
+				       //
+				       // 2. exchange the other two
+				       // components and change the
+				       // sign of one.
+				       //
+				       // this guarantees that the
+				       // scalar product of r and aux
+				       // is the square of the
+				       // smallest entry squared,
+				       // which is obviously less than
+				       // the sum of the entries
+				       // squared; the latter would
+				       // indicate that r and aux are
+				       // linearly dependent
+      Tensor<1,dim> aux = radius;
+      if ((std::fabs(radius[0]) < std::fabs(radius[1])) &&
+	  (std::fabs(radius[0]) < std::fabs(radius[2])))
+	{
+	  swap (aux[1], aux[2]);
+	  aux[1] = -aux[1];
+	}
+      else
+	if ((std::fabs(radius[1]) < std::fabs(radius[0])) &&
+	    (std::fabs(radius[1]) < std::fabs(radius[2])))
+	  {
+	    swap (aux[0], aux[2]);
+	    aux[0] = -aux[0];
+	  }
+	else
+	  {
+	    swap (aux[0], aux[1]);
+	    aux[0] = -aux[0];
+	  };
+
+				       // now construct the two
+				       // tangents by forming cross
+				       // products as discussed above
+      cross_product (face_vertex_tangents.tangents[vertex][0],
+		     radius, aux);
+      cross_product (face_vertex_tangents.tangents[vertex][1],
+		     face_vertex_tangents.tangents[vertex][0],
+		     radius);
+    };
+};
+
+#endif
+
+
 template <int dim>
 Point<dim>
 HyperBallBoundary<dim>::get_center () const 
@@ -319,9 +446,9 @@ get_intermediate_points_on_quad (const typename Triangulation<dim>::quad_iterato
 				   // boundary
   const Point<dim> quad_center = quad->center();
   if (quad_center(0) == center(0))
-    return StraightBoundary<dim>::get_intermediate_points_on_quad (quad, points);
+    StraightBoundary<dim>::get_intermediate_points_on_quad (quad, points);
   else
-    return HyperBallBoundary<dim>::get_intermediate_points_on_quad (quad, points);
+    HyperBallBoundary<dim>::get_intermediate_points_on_quad (quad, points);
 };
 
 
@@ -331,8 +458,8 @@ get_intermediate_points_on_quad (const typename Triangulation<dim>::quad_iterato
 template <>
 void
 HalfHyperBallBoundary<1>::
-get_intermediate_points_on_quad (const typename Triangulation<1>::quad_iterator &,
-				 typename std::vector<Point<1> > &) const
+get_intermediate_points_on_quad (const Triangulation<1>::quad_iterator &,
+				 std::vector<Point<1> > &) const
 {
   Assert (false, ExcInternalError());
 };
@@ -340,6 +467,25 @@ get_intermediate_points_on_quad (const typename Triangulation<1>::quad_iterator
 #endif
 
 
+
+template <int dim>
+void
+HyperBallBoundary<dim>::
+get_tangents_at_vertices (const typename Triangulation<dim>::face_iterator &face,
+			  typename Boundary<dim>::FaceVertexTangents &face_vertex_tangents) const
+{
+				   // check whether center of object is
+				   // at x==0, since then it belongs
+				   // to the plane part of the
+				   // boundary
+  const Point<dim> quad_center = quad->center();
+  if (quad_center(0) == center(0))
+    StraightBoundary<dim>::get_tangents_at_vertices (face, face_vertex_tangents);
+  else
+    HyperBallBoundary<dim>::get_intermediate_points_on_quad (face, face_vertex_tangents);
+};
+
+
 /* ---------------------------------------------------------------------- */
 
 
diff --git a/deal.II/doc/news/2001/c-3-1.html b/deal.II/doc/news/2001/c-3-1.html
index 88fc805124..e3675419f1 100644
--- a/deal.II/doc/news/2001/c-3-1.html
+++ b/deal.II/doc/news/2001/c-3-1.html
@@ -272,7 +272,6 @@ documentation, etc</a>.
        <br>
        (WB 2001/01/30)
        </p>
-
 </ol>
 
 
@@ -282,7 +281,17 @@ documentation, etc</a>.
 
 <ol>
   <li> <p>
+       New: <code class="class">Boundary</code> and derived classes
+       now have a function <code
+       class="member">get_tangents_at_vertices</code> that returns a
+       description of the tangent space to the boundary at the
+       vertices of a given face. This is used in the construction of
+       C<sup>1</sup> mappings of the boundary.
+       <br>
+       (WB 2001/03/23)
+       </p>
 
+  <li> <p>
        New: The new function <code
        class="class">PersistentTriangulation&lt;dim&gt;</code>::<code
        class="member">restore(unsigned int)</code> allows to restore