]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Boundary<dim>::get_tangents_at_vertices to generate a description of the tangent...
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 23 Mar 2001 14:48:42 +0000 (14:48 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 23 Mar 2001 14:48:42 +0000 (14:48 +0000)
git-svn-id: https://svn.dealii.org/trunk@4278 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/grid/tria_boundary.h
deal.II/deal.II/include/grid/tria_boundary_lib.h
deal.II/deal.II/source/grid/tria_boundary.cc
deal.II/deal.II/source/grid/tria_boundary_lib.cc
deal.II/doc/news/2001/c-3-1.html

index 650bb9a6802337c0d40b2b2ef612822f70a6e4ed..9eb1af8c9eb64170a129bf396ec0153f205d59ac 100644 (file)
@@ -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;
 };
 
 
index 3ccb5d45fa5197fe4b4e2f2563470b1c31e08f11..a972d61bd33cf46db47f493c87bc4c00dce0ec34 100644 (file)
@@ -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;
 };
 
 
index 5ef3b01668e13d2ea9672071fa392e84ee89da05..85c8e67fc92b95a9b80eab2687cc9a8d139b7093 100644 (file)
 #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>;
index aa98b7840624e40b6915ee389dd7cbf80ec1bb66..447d2b4b96ff600d9f1f91fb387edab2e0f1300a 100644 (file)
@@ -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);
+};
+
+
 /* ---------------------------------------------------------------------- */
 
 
index 88fc80512449cc1cfa293b7c924d596617f37ff8..e3675419f1a468e74d70148b47fe79c567998d93 100644 (file)
@@ -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

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.