From ff7cc4c8075b745bc835d0e80b0aa6028eca79e6 Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Sun, 20 Mar 2016 17:27:12 +0100 Subject: [PATCH] Moved normal_vector from Boundary to Manifold. --- include/deal.II/base/exceptions.h | 14 ++++ include/deal.II/grid/manifold.h | 67 +++++++++++++++ include/deal.II/grid/tria_boundary.h | 58 ------------- source/grid/manifold.cc | 119 +++++++++++++++++++++++++++ source/grid/manifold_lib.inst.in | 2 - source/grid/tria_boundary.cc | 25 ------ 6 files changed, 200 insertions(+), 85 deletions(-) diff --git a/include/deal.II/base/exceptions.h b/include/deal.II/base/exceptions.h index 722ab92ee0..c70ba18ccc 100644 --- a/include/deal.II/base/exceptions.h +++ b/include/deal.II/base/exceptions.h @@ -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 int, 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. */ diff --git a/include/deal.II/grid/manifold.h b/include/deal.II/grid/manifold.h index 30c650b9f4..2409a01e85 100644 --- a/include/deal.II/grid/manifold.h +++ b/include/deal.II/grid/manifold.h @@ -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 GeometryInfo::vertices_per_face + * 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::vertices_per_face]; + /** * Destructor. Does nothing here, but needs to be declared virtual to make @@ -382,6 +395,60 @@ public: const Point &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::face_iterator &face, + const Point &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::face_iterator &face, + FaceVertexNormals &face_vertex_normals) const; + + /// @} }; diff --git a/include/deal.II/grid/tria_boundary.h b/include/deal.II/grid/tria_boundary.h index 16cd34a611..da49bb72b8 100644 --- a/include/deal.II/grid/tria_boundary.h +++ b/include/deal.II/grid/tria_boundary.h @@ -82,19 +82,6 @@ class Boundary : public FlatManifold { public: - /** - * Type keeping information about the normals at the vertices of a face of a - * cell. Thus, there are GeometryInfo::vertices_per_face - * 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::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::face_iterator &face, std::vector > &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::face_iterator &face, - const Point &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::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 diff --git a/source/grid/manifold.cc b/source/grid/manifold.cc index 80f4aa9970..52cece6198 100644 --- a/source/grid/manifold.cc +++ b/source/grid/manifold.cc @@ -62,6 +62,125 @@ get_new_point (const Quadrature &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 +Tensor<1,spacedim> +Manifold:: +normal_vector (const typename Triangulation::face_iterator &face, + const Point &p) const +{ + Assert(false, ExcPureFunctionCalled()); + return Tensor<1,spacedim>(); +} + + +template +void +Manifold:: +get_normals_at_vertices (const typename Triangulation::face_iterator &face, + FaceVertexNormals &n) const +{ + for (unsigned int v=0; v::vertices_per_face; ++v) + n[v] = normal_vector(face, face->vertex(v)); +} + + + + template Point Manifold:: diff --git a/source/grid/manifold_lib.inst.in b/source/grid/manifold_lib.inst.in index eaaf53d08d..b8d74b7b67 100644 --- a/source/grid/manifold_lib.inst.in +++ b/source/grid/manifold_lib.inst.in @@ -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; -#endif -#if deal_II_dimension <= deal_II_space_dimension template class CylindricalManifold; template class FunctionManifold; template class FunctionManifold; diff --git a/source/grid/tria_boundary.cc b/source/grid/tria_boundary.cc index b1731fcc03..b1744a331d 100644 --- a/source/grid/tria_boundary.cc +++ b/source/grid/tria_boundary.cc @@ -106,31 +106,6 @@ get_intermediate_points_on_face (const Triangulation<1,3>::face_iterator &, } - - -template -Tensor<1,spacedim> -Boundary:: -normal_vector (const typename Triangulation::face_iterator &, - const Point &) const -{ - Assert (false, ExcPureFunctionCalled()); - return Tensor<1,spacedim>(); -} - - - -template -void -Boundary:: -get_normals_at_vertices (const typename Triangulation::face_iterator &, - FaceVertexNormals &) const -{ - Assert (false, ExcPureFunctionCalled()); -} - - - template Point Boundary:: -- 2.39.5