]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Move implementation of normal_vector and get_normal_at_vertices from StraightBoundary...
authorDaniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Tue, 2 Jan 2018 18:56:32 +0000 (19:56 +0100)
committerDaniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Fri, 5 Jan 2018 10:16:35 +0000 (11:16 +0100)
include/deal.II/grid/manifold.h
include/deal.II/grid/tria_boundary.h
source/grid/manifold.cc
source/grid/tria_boundary.cc

index e34149a72dea8089f388a37e54a3b207b0ebc76e..6f57df46e49213cacfef7059e63f9ed0fb6a59b7 100644 (file)
@@ -773,6 +773,31 @@ public:
   get_tangent_vector (const Point<spacedim> &x1,
                       const Point<spacedim> &x2) const override;
 
+  /**
+   * Return the normal vector to the given face at point p taking into account
+   * that quadrilateral faces of hexahedral cells in 3d may not be planar.
+   * In those cases, the face is assumed to have a geometry described by a
+   * bilinear function, and the normal vector is computed by embedding this
+   * bilinear form into a Cartesian space with a flat metric.
+   */
+  virtual
+  Tensor<1,spacedim>
+  normal_vector (const typename Triangulation<dim,spacedim>::face_iterator &face,
+                 const Point<spacedim> &p) const override;
+
+  /**
+   * Compute the normal vectors to the boundary at each vertex of the
+   * given face taking into account that quadrilateral faces of hexahedral
+   * cells in 3d may not be planar. In those cases, the face is assumed to
+   * have a geometry described by a bilinear function, and the normal vector
+   * is computed by embedding this bilinear form into a Cartesian space with
+   * a flat metric.
+   */
+  virtual
+  void
+  get_normals_at_vertices (const typename Triangulation<dim,spacedim>::face_iterator &face,
+                           typename Manifold<dim, spacedim>::FaceVertexNormals &face_vertex_normals) const override;
+
   /**
    * Return the periodicity of this Manifold.
    */
index d13e70246bc7225b65a6980ac89a34ec4238f26a..fe7053c73c09112b0fb432ed619d517084ec849f 100644 (file)
@@ -392,21 +392,6 @@ void
 Boundary<1,3>::
 get_intermediate_points_on_face (const Triangulation<1,3>::face_iterator &,
                                  std::vector<Point<3> > &) const;
-template <>
-void
-StraightBoundary<1,1>::
-get_normals_at_vertices (const Triangulation<1,1>::face_iterator &,
-                         Boundary<1,1>::FaceVertexNormals &) const;
-template <>
-void
-StraightBoundary<2,2>::
-get_normals_at_vertices (const Triangulation<2,2>::face_iterator &face,
-                         Boundary<2,2>::FaceVertexNormals &face_vertex_normals) const;
-template <>
-void
-StraightBoundary<3,3>::
-get_normals_at_vertices (const Triangulation<3,3>::face_iterator &face,
-                         Boundary<3,3>::FaceVertexNormals &face_vertex_normals) const;
 
 template <>
 Point<3>
index 13eb7ee81a3b9925a15e538e07e1c15598d811cb..b8bbbd2b3cb554ca0efdcaa70203555c20c49cf4 100644 (file)
@@ -476,6 +476,40 @@ Manifold<dim,spacedim>::get_tangent_vector(const Point<spacedim> &x1,
 
 /* -------------------------- FlatManifold --------------------- */
 
+namespace internal
+{
+  namespace
+  {
+    Tensor<1,2>
+    normalized_alternating_product (const Tensor<1,2> (&basis_vectors)[1])
+    {
+      Tensor<1,2> tmp = cross_product_2d (basis_vectors[0]);
+      return tmp/tmp.norm();
+    }
+
+
+
+    Tensor<1,3>
+    normalized_alternating_product (const Tensor<1,3> ( &)[1])
+    {
+      // we get here from FlatManifold<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_3d (basis_vectors[0], basis_vectors[1]);
+      return tmp/tmp.norm();
+    }
+
+  }
+}
 
 template <int dim, int spacedim>
 FlatManifold<dim,spacedim>::FlatManifold (const Tensor<1,spacedim> &periodicity,
@@ -617,6 +651,7 @@ get_new_points (const ArrayView<const Point<spacedim>> &surrounding_points,
 }
 
 
+
 template <int dim, int spacedim>
 Point<spacedim>
 FlatManifold<dim, spacedim>::project_to_manifold
@@ -662,6 +697,223 @@ FlatManifold<dim, spacedim>::get_tangent_vector (const Point<spacedim> &x1,
 
 
 
+template <>
+void
+FlatManifold<1>::
+get_normals_at_vertices (const Triangulation<1>::face_iterator &,
+                         Manifold<1,1>::FaceVertexNormals &) const
+{
+  Assert (false, ExcImpossibleInDim(1));
+}
+
+
+
+template <>
+void
+FlatManifold<1,2>::
+get_normals_at_vertices (const Triangulation<1,2>::face_iterator &,
+                         Manifold<1,2>::FaceVertexNormals &) const
+{
+  Assert (false, ExcNotImplemented());
+}
+
+
+
+template <>
+void
+FlatManifold<1,3>::
+get_normals_at_vertices (const Triangulation<1,3>::face_iterator &,
+                         Manifold<1,3>::FaceVertexNormals &) const
+{
+  Assert (false, ExcNotImplemented());
+}
+
+
+
+template <>
+void
+FlatManifold<2>::
+get_normals_at_vertices (const Triangulation<2>::face_iterator &face,
+                         Manifold<2,2>::FaceVertexNormals &face_vertex_normals) const
+{
+  const Tensor<1,2> tangent = face->vertex(1) - face->vertex(0);
+  for (unsigned int vertex=0; vertex<GeometryInfo<2>::vertices_per_face; ++vertex)
+    // compute normals from tangent
+    face_vertex_normals[vertex] = Point<2>(tangent[1],
+                                           -tangent[0]);
+}
+
+
+
+template <>
+void
+FlatManifold<2,3>::
+get_normals_at_vertices (const Triangulation<2,3>::face_iterator &face,
+                         Manifold<2,3>::FaceVertexNormals &face_vertex_normals) const
+{
+  Assert(false, ExcNotImplemented());
+}
+
+
+
+template <>
+void
+FlatManifold<3>::
+get_normals_at_vertices (const Triangulation<3>::face_iterator &face,
+                         Manifold<3,3>::FaceVertexNormals &face_vertex_normals) const
+{
+  const unsigned int vertices_per_face = GeometryInfo<3>::vertices_per_face;
+
+  static const unsigned int neighboring_vertices[4][2]=
+  { {1,2},{3,0},{0,3},{2,1}};
+  for (unsigned int vertex=0; vertex<vertices_per_face; ++vertex)
+    {
+      // first define the two tangent vectors at the vertex by using the
+      // two lines radiating away from this vertex
+      const Tensor<1,3> tangents[2]
+        = { face->vertex(neighboring_vertices[vertex][0])
+            - face->vertex(vertex),
+            face->vertex(neighboring_vertices[vertex][1])
+            - face->vertex(vertex)
+          };
+
+      // then compute the normal by taking the cross product. since the
+      // normal is not required to be normalized, no problem here
+      face_vertex_normals[vertex] = cross_product_3d(tangents[0], tangents[1]);
+    }
+}
+
+
+
+template <>
+Tensor<1,1>
+FlatManifold<1,1>::
+normal_vector (const Triangulation<1,1>::face_iterator &,
+               const Point<1> &) const
+{
+  Assert (false, ExcNotImplemented());
+  return Tensor<1,1>();
+}
+
+
+
+template <>
+Tensor<1,2>
+FlatManifold<1,2>::
+normal_vector (const Triangulation<1,2>::face_iterator &,
+               const Point<2> &) const
+{
+  Assert (false, ExcNotImplemented());
+  return Tensor<1,2>();
+}
+
+
+
+template <>
+Tensor<1,3>
+FlatManifold<1,3>::
+normal_vector (const Triangulation<1,3>::face_iterator &,
+               const Point<3> &) const
+{
+  Assert (false, ExcNotImplemented());
+  return Tensor<1,3>();
+}
+
+
+
+template <int dim, int spacedim>
+Tensor<1, spacedim>
+FlatManifold< dim, spacedim >::
+normal_vector (const typename Triangulation<dim, spacedim>::face_iterator &face,
+               const Point<spacedim>                                      &p) const
+{
+  // 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)
+  // 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;
+
+  const double eps = 1e-12;
+  Tensor<1,spacedim> grad_F[facedim];
+  unsigned int iteration = 0;
+  while (true)
+    {
+      Point<spacedim> F;
+      for (unsigned int v=0; v<GeometryInfo<facedim>::vertices_per_cell; ++v)
+        F += face->vertex(v) * GeometryInfo<facedim>::d_linear_shape_function(xi, v);
+
+      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] += face->vertex(v) *
+                         GeometryInfo<facedim>::d_linear_shape_function_gradient(xi, v)[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 Tensor<1,facedim> delta_xi = -invert(H) * J;
+      xi += delta_xi;
+      ++iteration;
+
+      Assert (iteration<10,
+              ExcMessage("The Newton iteration to find the reference point "
+                         "did not converge in 10 iterations. Do you have a "
+                         "deformed cell? (See the glossary for a definition "
+                         "of what a deformed cell is. You may want to output "
+                         "the vertices of your cell."));
+
+      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);
+}
+
+
 /* -------------------------- ChartManifold --------------------- */
 template <int dim, int spacedim, int chartdim>
 ChartManifold<dim,spacedim,chartdim>::ChartManifold (const Tensor<1,chartdim> &periodicity)
index 250dd5dbae7604e6d39932a3c876f7bdd79b34b8..7fb26fe4abb04db3b9ed0fe579a574fcfd2986d6 100644 (file)
@@ -418,258 +418,24 @@ get_intermediate_points_on_quad (const Triangulation<2,3>::quad_iterator &quad,
 
 
 
-template <>
-Tensor<1,1>
-StraightBoundary<1,1>::
-normal_vector (const Triangulation<1,1>::face_iterator &,
-               const Point<1> &) const
-{
-  Assert (false, ExcNotImplemented());
-  return Tensor<1,1>();
-}
-
-
-template <>
-Tensor<1,2>
-StraightBoundary<1,2>::
-normal_vector (const Triangulation<1,2>::face_iterator &,
-               const Point<2> &) const
-{
-  Assert (false, ExcNotImplemented());
-  return Tensor<1,2>();
-}
-
-
-template <>
-Tensor<1,3>
-StraightBoundary<1,3>::
-normal_vector (const Triangulation<1,3>::face_iterator &,
-               const Point<3> &) const
-{
-  Assert (false, ExcNotImplemented());
-  return Tensor<1,3>();
-}
-
-
-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_2d (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_3d (basis_vectors[0], basis_vectors[1]);
-      return tmp/tmp.norm();
-    }
-
-  }
-}
-
-
 template <int dim, int spacedim>
 Tensor<1,spacedim>
 StraightBoundary<dim,spacedim>::
 normal_vector (const typename Triangulation<dim,spacedim>::face_iterator &face,
-               const Point<spacedim> &p) const
-{
-  // 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)
-  // 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;
-
-  const double eps = 1e-12;
-  Tensor<1,spacedim> grad_F[facedim];
-  unsigned int iteration = 0;
-  while (true)
-    {
-      Point<spacedim> F;
-      for (unsigned int v=0; v<GeometryInfo<facedim>::vertices_per_cell; ++v)
-        F += face->vertex(v) * GeometryInfo<facedim>::d_linear_shape_function(xi, v);
-
-      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] += face->vertex(v) *
-                         GeometryInfo<facedim>::d_linear_shape_function_gradient(xi, v)[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 Tensor<1,facedim> delta_xi = -invert(H) * J;
-      xi += delta_xi;
-      ++iteration;
-
-      Assert (iteration<10,
-              ExcMessage("The Newton iteration to find the reference point "
-                         "did not converge in 10 iterations. Do you have a "
-                         "deformed cell? (See the glossary for a definition "
-                         "of what a deformed cell is. You may want to output "
-                         "the vertices of your cell."));
-
-      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);
-}
-
-
-
-template <>
-void
-StraightBoundary<1>::
-get_normals_at_vertices (const Triangulation<1>::face_iterator &,
-                         Boundary<1,1>::FaceVertexNormals &) const
-{
-  Assert (false, ExcImpossibleInDim(1));
-}
-
-template <>
-void
-StraightBoundary<1,2>::
-get_normals_at_vertices (const Triangulation<1,2>::face_iterator &,
-                         Boundary<1,2>::FaceVertexNormals &) const
-{
-  Assert (false, ExcNotImplemented());
-}
-
-
-template <>
-void
-StraightBoundary<1,3>::
-get_normals_at_vertices (const Triangulation<1,3>::face_iterator &,
-                         Boundary<1,3>::FaceVertexNormals &) const
+               const Point<spacedim>                                     &p) const
 {
-  Assert (false, ExcNotImplemented());
+  return FlatManifold<dim, spacedim>::normal_vector(face, p);
 }
 
 
 
-template <>
-void
-StraightBoundary<2>::
-get_normals_at_vertices (const Triangulation<2>::face_iterator &face,
-                         Boundary<2,2>::FaceVertexNormals &face_vertex_normals) const
-{
-  const Tensor<1,2> tangent = face->vertex(1) - face->vertex(0);
-  for (unsigned int vertex=0; vertex<GeometryInfo<2>::vertices_per_face; ++vertex)
-    // compute normals from tangent
-    face_vertex_normals[vertex] = Point<2>(tangent[1],
-                                           -tangent[0]);
-}
-
-template <>
-void
-StraightBoundary<2,3>::
-get_normals_at_vertices (const Triangulation<2,3>::face_iterator &face,
-                         Boundary<2,3>::FaceVertexNormals &face_vertex_normals) const
-{
-  const Tensor<1,3> tangent = face->vertex(1) - face->vertex(0);
-  for (unsigned int vertex=0; vertex<GeometryInfo<2>::vertices_per_face; ++vertex)
-    // compute normals from tangent
-    face_vertex_normals[vertex] = Point<3>(tangent[1],
-                                           -tangent[0],0);
-  Assert(false, ExcNotImplemented());
-}
-
-
-
-
-template <>
+template <int dim, int spacedim>
 void
-StraightBoundary<3>::
-get_normals_at_vertices (const Triangulation<3>::face_iterator &face,
-                         Boundary<3,3>::FaceVertexNormals &face_vertex_normals) const
+StraightBoundary<dim, spacedim>::
+get_normals_at_vertices (const typename Triangulation<dim, spacedim>::face_iterator &face,
+                         typename Boundary<dim,spacedim>::FaceVertexNormals         &face_vertex_normals) const
 {
-  const unsigned int vertices_per_face = GeometryInfo<3>::vertices_per_face;
-
-  static const unsigned int neighboring_vertices[4][2]=
-  { {1,2},{3,0},{0,3},{2,1}};
-  for (unsigned int vertex=0; vertex<vertices_per_face; ++vertex)
-    {
-      // first define the two tangent vectors at the vertex by using the
-      // two lines radiating away from this vertex
-      const Tensor<1,3> tangents[2]
-        = { face->vertex(neighboring_vertices[vertex][0])
-            - face->vertex(vertex),
-            face->vertex(neighboring_vertices[vertex][1])
-            - face->vertex(vertex)
-          };
-
-      // then compute the normal by taking the cross product. since the
-      // normal is not required to be normalized, no problem here
-      face_vertex_normals[vertex] = cross_product_3d(tangents[0], tangents[1]);
-    };
+  FlatManifold<dim, spacedim>::get_normals_at_vertices(face, face_vertex_normals);
 }
 
 

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.