]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Moved normal_vector from Boundary to Manifold.
authorLuca Heltai <luca.heltai@sissa.it>
Sun, 20 Mar 2016 16:27:12 +0000 (17:27 +0100)
committerLuca Heltai <luca.heltai@sissa.it>
Sun, 20 Mar 2016 16:33:21 +0000 (17:33 +0100)
include/deal.II/base/exceptions.h
include/deal.II/grid/manifold.h
include/deal.II/grid/tria_boundary.h
source/grid/manifold.cc
source/grid/manifold_lib.inst.in
source/grid/tria_boundary.cc

index 722ab92ee0c31f5a3232a62057449a7f68c5e0d1..c70ba18ccc15516677209effa4277f0c805addda 100644 (file)
@@ -835,6 +835,20 @@ namespace StandardExceptions
                   << "impossible in " << arg1
                   << "d or simply does not make any sense.");
 
+  /**
+   * This exception is raised if a functionality is not possible in the given
+   * combination of dimension and space-dimension.
+   *
+   * The constructor takes two <tt>int</tt>, denoting the dimension and the
+   * space dimension.
+   */
+  DeclException2 (ExcImpossibleInDimSpacedim,
+                  int, int,
+                  << "You are trying to execute functionality that is "
+                  << "impossible in dimensions <" << arg1 << "," << arg2
+                  << "> or simply does not make any sense.");
+
+
   /**
    * A number is zero, but it should not be here.
    */
index 30c650b9f4b5950b5ea434e63ae6906570e2fd5c..2409a01e8531c3fd3ab0fa5e39c69d4e57098106 100644 (file)
@@ -212,6 +212,19 @@ class Manifold : public Subscriptor
 {
 public:
 
+  /**
+   * Type keeping information about the normals at the vertices of a face of a
+   * cell. Thus, there are <tt>GeometryInfo<dim>::vertices_per_face</tt>
+   * normal vectors, that define the tangent spaces of the boundary at the
+   * vertices. Note that the vectors stored in this object are not required to
+   * be normalized, nor to actually point outward, as one often will only want
+   * to check for orthogonality to define the tangent plane; if a function
+   * requires the normals to be normalized, then it must do so itself.
+   *
+   * For obvious reasons, this type is not useful in 1d.
+   */
+  typedef Tensor<1,spacedim> FaceVertexNormals[GeometryInfo<dim>::vertices_per_face];
+
 
   /**
    * Destructor. Does nothing here, but needs to be declared virtual to make
@@ -382,6 +395,60 @@ public:
                       const Point<spacedim> &x2) const;
 
   /// @}
+
+  /**
+   * @name Computing normal vectors
+   */
+  /// @{
+
+  /**
+   * Return the normal vector to a face embedded in this manifold, at
+   * the point p. If p is not in fact on the surface, but only
+   * close-by, try to return something reasonable, for example the
+   * normal vector at the surface point closest to p.  (The point p
+   * will in fact not normally lie on the actual surface, but rather
+   * be a quadrature point mapped by some polynomial mapping; the
+   * mapped surface, however, will not usually coincide with the
+   * actual surface.)
+   *
+   * The face iterator gives an indication which face this function is
+   * supposed to compute the normal vector for.  This is useful if the
+   * boundary of the domain is composed of different nondifferential
+   * pieces (for example when using the StraightBoundary class to
+   * approximate a geometry that is completely described by the coarse
+   * mesh, with piecewise (bi-)linear components between the vertices,
+   * but where the boundary may have a kink at the vertices itself).
+   *
+   * @note The default implementation of this function computes the
+   * normal vector by taking the cross product between the tangent
+   * vectors from p to the most orthogonal and further non consecutive
+   * vertices of the face.
+   */
+  virtual
+  Tensor<1,spacedim>
+  normal_vector (const typename Triangulation<dim,spacedim>::face_iterator &face,
+                 const Point<spacedim> &p) const;
+
+  /**
+   * Compute the normal vectors to the boundary at each vertex of the
+   * given face embedded in the Manifold. It is not required that the
+   * normal vectors be normed somehow.  Neither is it required that
+   * the normals actually point outward.
+   *
+   * This function is needed to compute data for C1 mappings. The
+   * default implementation calls normal_vector() on each vertex.
+   *
+   * Note that when computing normal vectors at a vertex where the
+   * boundary is not differentiable, you have to make sure that you
+   * compute the one-sided limits, i.e. limit with respect to points
+   * inside the given face.
+   */
+  virtual
+  void
+  get_normals_at_vertices (const typename Triangulation<dim,spacedim>::face_iterator &face,
+                           FaceVertexNormals &face_vertex_normals) const;
+
+  /// @}
 };
 
 
index 16cd34a61180a16cebd91777ec71a1128c796a02..da49bb72b87e400aef3395a57b2b8f3df5eb252d 100644 (file)
@@ -82,19 +82,6 @@ class Boundary : public FlatManifold<dim, spacedim>
 {
 public:
 
-  /**
-   * Type keeping information about the normals at the vertices of a face of a
-   * cell. Thus, there are <tt>GeometryInfo<dim>::vertices_per_face</tt>
-   * normal vectors, that define the tangent spaces of the boundary at the
-   * vertices. Note that the vectors stored in this object are not required to
-   * be normalized, nor to actually point outward, as one often will only want
-   * to check for orthogonality to define the tangent plane; if a function
-   * requires the normals to be normalized, then it must do so itself.
-   *
-   * For obvious reasons, this type is not useful in 1d.
-   */
-  typedef Tensor<1,spacedim> FaceVertexNormals[GeometryInfo<dim>::vertices_per_face];
-
   /**
    * Destructor. Does nothing here, but needs to be declared to make it
    * virtual.
@@ -160,51 +147,6 @@ public:
   get_intermediate_points_on_face (const typename Triangulation<dim,spacedim>::face_iterator &face,
                                    std::vector<Point<spacedim> > &points) const;
 
-  /**
-   * Return the normal vector to the surface at the point p. If p is not in
-   * fact on the surface, but only close-by, try to return something
-   * reasonable, for example the normal vector at the surface point closest to
-   * p.  (The point p will in fact not normally lie on the actual surface, but
-   * rather be a quadrature point mapped by some polynomial mapping; the
-   * mapped surface, however, will not usually coincide with the actual
-   * surface.)
-   *
-   * The face iterator gives an indication which face this function is
-   * supposed to compute the normal vector for.  This is useful if the
-   * boundary of the domain is composed of different nondifferential pieces
-   * (for example when using the StraightBoundary class to approximate a
-   * geometry that is completely described by the coarse mesh, with piecewise
-   * (bi-)linear components between the vertices, but where the boundary may
-   * have a kink at the vertices itself).
-   *
-   * @note Implementations of this function should be able to assume that the
-   * point p lies within or close to the face described by the first argument.
-   * In turn, callers of this function should ensure that this is in fact the
-   * case.
-   */
-  virtual
-  Tensor<1,spacedim>
-  normal_vector (const typename Triangulation<dim,spacedim>::face_iterator &face,
-                 const Point<spacedim> &p) const;
-
-  /**
-   * Compute the normal vectors to the boundary at each vertex of the given
-   * face. It is not required that the normal vectors be normed somehow.
-   * Neither is it required that the normals actually point outward.
-   *
-   * 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 do not intend to use C1 mappings.
-   *
-   * Note that when computing normal vectors at a vertex where the boundary is
-   * not differentiable, you have to make sure that you compute the one-sided
-   * limits, i.e. limit with respect to points inside the given face.
-   */
-  virtual
-  void
-  get_normals_at_vertices (const typename Triangulation<dim,spacedim>::face_iterator &face,
-                           FaceVertexNormals &face_vertex_normals) const;
-
   /**
    * Given a candidate point and a line segment characterized by the iterator,
    * return a point that lies on the surface described by this object. This
index 80f4aa997059e492cf69218de15386e786eb7e13..52cece6198804b601cf354e53e71b5b57db66673 100644 (file)
@@ -62,6 +62,125 @@ get_new_point (const Quadrature<spacedim> &quad) const
 }
 
 
+
+template <>
+Tensor<1,2>
+Manifold<2, 2>::
+normal_vector (const typename Triangulation<2, 2>::face_iterator &face,
+               const Point<2> &p) const
+{
+  const int spacedim=2, dim=2;
+
+  Tensor<1,spacedim> tangent = ((p-face->vertex(0)).norm_square() > (p-face->vertex(1)).norm_square() ?
+                                get_tangent_vector(p, face->vertex(0)) :
+                                get_tangent_vector(p, face->vertex(1)));
+  return cross_product_2d(tangent);
+}
+
+template<>
+Tensor<1,3>
+Manifold<3, 3>::
+normal_vector (const typename Triangulation<3, 3>::face_iterator &face,
+               const Point<3> &p) const
+{
+  const int spacedim=3, dim=3;
+  Tensor<1,spacedim> t1,t2;
+
+  // Take the difference between p and all four vertices
+  Tensor<1,spacedim> dp[4];
+  double dpns[4];
+  int min_index=-1;
+  double min_distance = 0;
+  for (unsigned int i=0; i<4; ++i)
+    {
+      dp[i] = p-face->vertex(i);
+      dpns[i] = dp[i].norm_square();
+      min_index = (min_index == -1 ? (int)i : dpns[i] < min_distance ? i : min_index);
+      min_distance = dpns[min_index];
+    }
+  // Verify we have a valid vertex index
+  AssertIndexRange(min_index, 4);
+
+  // now figure out which vertices are better to compute tangent vectors
+  // we split the cell in 4 quadrants, and use v1/v2 for the first quadrant
+  // (the one at the corner with vertex 0)
+  if ((p-face->center()).norm_square() < min_distance)
+    {
+      // we are close to the face center: pick two consecutive vertices,
+      // but not the closest one. We make sure the direction is always
+      // the same.
+      if (min_index < 2)
+        {
+          t1 = get_tangent_vector(p, face->vertex(3));
+          t2 = get_tangent_vector(p, face->vertex(2));
+        }
+      else
+        {
+          t1 = get_tangent_vector(p, face->vertex(0));
+          t2 = get_tangent_vector(p, face->vertex(1));
+        }
+    }
+  else
+    {
+      switch (min_index)
+        {
+        case 0:
+        {
+          t1 = get_tangent_vector(p, face->vertex(1));
+          t2 = get_tangent_vector(p, face->vertex(2));
+          break;
+        }
+        case 1:
+        {
+          t1 = get_tangent_vector(p, face->vertex(3));
+          t2 = get_tangent_vector(p, face->vertex(0));
+          break;
+        }
+        case 2:
+        {
+          t1 = get_tangent_vector(p, face->vertex(0));
+          t2 = get_tangent_vector(p, face->vertex(3));
+          break;
+        }
+        case 3:
+        {
+          t1 = get_tangent_vector(p, face->vertex(2));
+          t2 = get_tangent_vector(p, face->vertex(1));
+          break;
+        }
+        default:
+          Assert(false, ExcInternalError());
+          break;
+        }
+    }
+  return cross_product_3d(t1,t2);
+}
+
+
+template <int dim, int spacedim>
+Tensor<1,spacedim>
+Manifold<dim, spacedim>::
+normal_vector (const typename Triangulation<dim, spacedim>::face_iterator &face,
+               const Point<spacedim> &p) const
+{
+  Assert(false, ExcPureFunctionCalled());
+  return Tensor<1,spacedim>();
+}
+
+
+template <int dim, int spacedim>
+void
+Manifold<dim, spacedim>::
+get_normals_at_vertices (const typename Triangulation<dim, spacedim>::face_iterator &face,
+                         FaceVertexNormals &n) const
+{
+  for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_face; ++v)
+    n[v] = normal_vector(face, face->vertex(v));
+}
+
+
+
+
 template <int dim, int spacedim>
 Point<spacedim>
 Manifold<dim, spacedim>::
index eaaf53d08dd82f00600a85a2f0faae02de83deae..b8d74b7b6789ac18490891faa9f13cc5e24e7456 100644 (file)
@@ -18,8 +18,6 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS
   {
 #if deal_II_dimension <= deal_II_space_dimension
     template class SphericalManifold<deal_II_dimension, deal_II_space_dimension>;
-#endif
-#if deal_II_dimension <= deal_II_space_dimension
     template class CylindricalManifold<deal_II_dimension, deal_II_space_dimension>;
     template class FunctionManifold<deal_II_dimension, deal_II_space_dimension, 1>;                      
     template class FunctionManifold<deal_II_dimension, deal_II_space_dimension, 2>;
index b1731fcc03b06ba43e6efed52d84e3aaeb644691..b1744a331d2f75ce14caf2b8e845b5a9b73f2135 100644 (file)
@@ -106,31 +106,6 @@ get_intermediate_points_on_face (const Triangulation<1,3>::face_iterator &,
 }
 
 
-
-
-template <int dim, int spacedim>
-Tensor<1,spacedim>
-Boundary<dim, spacedim>::
-normal_vector (const typename Triangulation<dim, spacedim>::face_iterator &,
-               const Point<spacedim> &) const
-{
-  Assert (false, ExcPureFunctionCalled());
-  return Tensor<1,spacedim>();
-}
-
-
-
-template <int dim, int spacedim>
-void
-Boundary<dim, spacedim>::
-get_normals_at_vertices (const typename Triangulation<dim, spacedim>::face_iterator &,
-                         FaceVertexNormals &) const
-{
-  Assert (false, ExcPureFunctionCalled());
-}
-
-
-
 template <int dim, int spacedim>
 Point<spacedim>
 Boundary<dim, spacedim>::

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.