From: wolf Date: Thu, 29 Mar 2001 08:52:23 +0000 (+0000) Subject: Compute C1 mapping from normal vectors, rather than from tangentials, as they are... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=022661684158592270d794f1088d86bd8c773493;p=dealii-svn.git Compute C1 mapping from normal vectors, rather than from tangentials, as they are simpler to compute in the beginning, and their future use in 3d is less uncertain than the use of tangentials, as e.g., we don't need both tangentials in the middle of a line in 3d, as the tangential along the line is already determined, we will only need the one perpendicular to the line. git-svn-id: https://svn.dealii.org/trunk@4319 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 9eb1af8c9e..54a37755dd 100644 --- a/deal.II/deal.II/include/grid/tria_boundary.h +++ b/deal.II/deal.II/include/grid/tria_boundary.h @@ -70,7 +70,7 @@ template class Triangulation; * @ref{HyperBallBoundary} creating a hyperball with given radius * around a given center point. * - * @author Wolfgang Bangerth, 1999, Ralf Hartmann, 2001 + * @author Wolfgang Bangerth, 1999, 2001, Ralf Hartmann, 2001 */ template class Boundary : public Subscriptor @@ -78,45 +78,29 @@ class Boundary : public Subscriptor public: /** - * Structure keeping information - * about the tangents at the - * vertices of a face of a - * cell. Thus, there are + * Type keeping information about + * the normals at the vertices of + * a face of a cell. Thus, there + * are * @p{GeometryInfo::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. + * 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 - * class is not useful in 1d. + * type 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(-1)); - - /** - * Array of tangent vectors, - * as described above. - */ - Tensor<1,dim> tangents[GeometryInfo::vertices_per_face][n_tangents]; - }; - + typedef Tensor<1,dim> FaceVertexNormals[GeometryInfo::vertices_per_face]; /** * Destructor. Does nothing here, but @@ -126,37 +110,43 @@ class Boundary : public Subscriptor virtual ~Boundary (); /** - * Return the point which shall become - * the new middle vertex of the two - * children of a regular line. In 2D, - * this line is a line at the boundary, - * while in 3d, it is bounding a face - * at the boundary (the lines therefore - * is also on the boundary). + * Return the point which shall + * become the new middle vertex + * of the two children of a + * regular line. In 2D, this line + * is a line at the boundary, + * while in 3d, it is bounding a + * face at the boundary (the + * lines therefore is also on the + * boundary). */ virtual Point get_new_point_on_line (const typename Triangulation::line_iterator &line) const = 0; /** - * Return the point which shall become - * the common point of the four children - * of a quad at the boundary in three - * or more spatial dimensions. This - * function therefore is only useful in - * at least three dimensions and should - * not be called for lower dimensions. + * Return the point which shall + * become the common point of the + * four children of a quad at the + * boundary in three or more + * spatial dimensions. This + * function therefore is only + * useful in at least three + * dimensions and should not be + * called for lower dimensions. * - * This function is called after the - * four lines bounding the given @p{quad} - * are refined, so you may want to use - * the information provided by - * @p{quad->line(i)->child(j)}, @p{i=0...3}, - * @p{j=0,1}. + * This function is called after + * the four lines bounding the + * given @p{quad} are refined, so + * you may want to use the + * information provided by + * @p{quad->line(i)->child(j)}, + * @p{i=0...3}, @p{j=0,1}. * - * Because in 2D, this function is not - * needed, it is not made pure virtual, - * to avoid the need to overload it. - * The default implementation throws + * Because in 2D, this function + * is not needed, it is not made + * pure virtual, to avoid the + * need to overload it. The + * default implementation throws * an error in any case, however. */ virtual Point @@ -225,18 +215,14 @@ class Boundary : public Subscriptor typename std::vector > &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. + * 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 @@ -248,17 +234,17 @@ class Boundary : public Subscriptor * mappings. * * Note that when computing - * tangents at a vertex where the - * boundary is not + * normal vectors 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. + * one-sided limits, i.e. limit + * with respect to points inside + * the given face. */ virtual void - get_tangents_at_vertices (const typename Triangulation::face_iterator &face, - FaceVertexTangents &face_vertex_tangents) const; + get_normals_at_vertices (const typename Triangulation::face_iterator &face, + FaceVertexNormals &face_vertex_normals) const; /** @@ -277,42 +263,50 @@ class Boundary : public Subscriptor /** - * Specialisation of @ref{Boundary}, which places the new point right - * into the middle of the given points. The middle is defined as the - * arithmetic mean of the points. + * Specialisation of @ref{Boundary}, which places the new point + * right into the middle of the given points. The middle is defined + * as the arithmetic mean of the points. * - * This class does not really describe a boundary in the usual sense. By - * placing new points in the middle of old ones, it rather assumes that the - * boundary of the domain is given by the polygon/polyhedron defined by the - * boundary of the initial coarse triangulation. + * This class does not really describe a boundary in the usual + * sense. By placing new points in the middle of old ones, it rather + * assumes that the boundary of the domain is given by the + * polygon/polyhedron defined by the boundary of the initial coarse + * triangulation. * - * @author Wolfgang Bangerth, 1998, Ralf Hartmann, 2001 + * @author Wolfgang Bangerth, 1998, 2001, Ralf Hartmann, 2001 */ template class StraightBoundary : public Boundary { public: /** - * Let the new point be the arithmetic - * mean of the two vertices of the line. + * Let the new point be the + * arithmetic mean of the two + * vertices of the line. * - * Refer to the general documentation of - * this class and the documentation of the - * base class for more information. + * Refer to the general + * documentation of this class + * and the documentation of the + * base class for more + * information. */ virtual Point get_new_point_on_line (const typename Triangulation::line_iterator &line) const; /** - * Let the new point be the arithmetic mean - * of the four vertices of this quad and - * the four midpoints of the lines, which - * are already created at the time of calling - * this function. + * Let the new point be the + * arithmetic mean of the four + * vertices of this quad and the + * four midpoints of the lines, + * which are already created at + * the time of calling this + * function. * - * Refer to the general documentation of - * this class and the documentation of the - * base class for more information. + * Refer to the general + * documentation of this class + * and the documentation of the + * base class for more + * information. */ virtual Point get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const; @@ -350,9 +344,9 @@ class StraightBoundary : public Boundary typename std::vector > &points) const; /** - * Compute a basis of the tangent - * space at each vertex of the - * given face. + * Compute the normals to the + * boundary at the vertices of + * the given face. * * Refer to the general * documentation of this class @@ -360,8 +354,8 @@ class StraightBoundary : public Boundary * base class. */ virtual void - get_tangents_at_vertices (const typename Triangulation::face_iterator &face, - typename Boundary::FaceVertexTangents &face_vertex_tangents) const; + get_normals_at_vertices (const typename Triangulation::face_iterator &face, + typename Boundary::FaceVertexNormals &face_vertex_normals) 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 a972d61bd3..92cb636666 100644 --- a/deal.II/deal.II/include/grid/tria_boundary_lib.h +++ b/deal.II/deal.II/include/grid/tria_boundary_lib.h @@ -85,9 +85,9 @@ class HyperBallBoundary : public StraightBoundary typename std::vector > &points) const; /** - * Compute a basis of the tangent - * space at each vertex of the - * given face. + * Compute the normals to the + * boundary at the vertices of + * the given face. * * Refer to the general * documentation of this class @@ -95,8 +95,8 @@ class HyperBallBoundary : public StraightBoundary * base class. */ virtual void - get_tangents_at_vertices (const typename Triangulation::face_iterator &face, - typename Boundary::FaceVertexTangents &face_vertex_tangents) const; + get_normals_at_vertices (const typename Triangulation::face_iterator &face, + typename Boundary::FaceVertexNormals &face_vertex_normals) const; /** * Return the center of the ball. @@ -201,11 +201,10 @@ class HalfHyperBallBoundary : public HyperBallBoundary get_intermediate_points_on_quad (const typename Triangulation::quad_iterator &quad, typename std::vector > &points) const; - /** - * Compute a basis of the tangent - * space at each vertex of the - * given face. + * Compute the normals to the + * boundary at the vertices of + * the given face. * * Refer to the general * documentation of this class @@ -213,8 +212,8 @@ class HalfHyperBallBoundary : public HyperBallBoundary * base class. */ virtual void - get_tangents_at_vertices (const typename Triangulation::face_iterator &face, - typename Boundary::FaceVertexTangents &face_vertex_tangents) const; + get_normals_at_vertices (const typename Triangulation::face_iterator &face, + typename Boundary::FaceVertexNormals &face_vertex_normals) const; }; diff --git a/deal.II/deal.II/source/fe/mapping_c1.cc b/deal.II/deal.II/source/fe/mapping_c1.cc index ca846fc36e..cc20ddb6f3 100644 --- a/deal.II/deal.II/source/fe/mapping_c1.cc +++ b/deal.II/deal.II/source/fe/mapping_c1.cc @@ -68,7 +68,7 @@ MappingC1<2>::add_line_support_points (const Triangulation<2>::cell_iterator &ce if (line->at_boundary()) { - // first get the tangential + // first get the normal // vectors at the two // vertices of this line // from the boundary @@ -76,8 +76,8 @@ MappingC1<2>::add_line_support_points (const Triangulation<2>::cell_iterator &ce const Boundary &boundary = line->get_triangulation().get_boundary(line->boundary_indicator()); - Boundary::FaceVertexTangents face_vertex_tangents; - boundary.get_tangents_at_vertices (line, face_vertex_tangents); + Boundary::FaceVertexNormals face_vertex_normals; + boundary.get_normals_at_vertices (line, face_vertex_normals); // then transform them into // interpolation points for @@ -118,7 +118,7 @@ MappingC1<2>::add_line_support_points (const Triangulation<2>::cell_iterator &ce // coefficients from the // tangentials. for that, // first rotate the - // tangentials of @p{s(t)} + // tangents of @p{s(t)} // into the global // coordinate system. they // are @p{A (1,c)} and @p{A @@ -133,7 +133,7 @@ MappingC1<2>::add_line_support_points (const Triangulation<2>::cell_iterator &ce // have to make sure by // matching @p{b,c} that // these tangentials are - // parallel to those + // orthogonal to the normals // returned by the boundary // object const Tensor<1,2> coordinate_vector = line->vertex(1) - line->vertex(0); @@ -142,14 +142,14 @@ MappingC1<2>::add_line_support_points (const Triangulation<2>::cell_iterator &ce coordinate_axis /= h; const double alpha = std::atan2(coordinate_axis[1], coordinate_axis[0]); - const double b = ((face_vertex_tangents.tangents[0][0][0] * sin(alpha) - -face_vertex_tangents.tangents[0][0][1] * cos(alpha)) / - (face_vertex_tangents.tangents[0][0][0] * cos(alpha) - +face_vertex_tangents.tangents[0][0][1] * sin(alpha))), - c = ((face_vertex_tangents.tangents[1][0][0] * sin(alpha) - -face_vertex_tangents.tangents[1][0][1] * cos(alpha)) / - (face_vertex_tangents.tangents[1][0][0] * cos(alpha) - +face_vertex_tangents.tangents[1][0][1] * sin(alpha))); + const double b = ((face_vertex_normals[0][1] * sin(alpha) + +face_vertex_normals[0][0] * cos(alpha)) / + (face_vertex_normals[0][1] * cos(alpha) + -face_vertex_normals[0][0] * sin(alpha))), + c = ((face_vertex_normals[1][1] * sin(alpha) + +face_vertex_normals[1][0] * cos(alpha)) / + (face_vertex_normals[1][1] * cos(alpha) + -face_vertex_normals[1][0] * sin(alpha))); // next evaluate the so // determined cubic diff --git a/deal.II/deal.II/source/grid/tria_boundary.cc b/deal.II/deal.II/source/grid/tria_boundary.cc index 85c8e67fc9..1689110a53 100644 --- a/deal.II/deal.II/source/grid/tria_boundary.cc +++ b/deal.II/deal.II/source/grid/tria_boundary.cc @@ -64,8 +64,8 @@ get_intermediate_points_on_quad (const typename Triangulation::quad_iterato template void Boundary:: -get_tangents_at_vertices (const typename Triangulation::face_iterator &, - FaceVertexTangents &) const +get_normals_at_vertices (const typename Triangulation::face_iterator &, + FaceVertexNormals &) const { Assert (false, ExcPureVirtualFunctionCalled()); }; @@ -203,8 +203,8 @@ get_intermediate_points_on_quad (const typename Triangulation::quad_iterato template <> void StraightBoundary<1>:: -get_tangents_at_vertices (const Triangulation<1>::face_iterator &, - FaceVertexTangents &) const +get_normals_at_vertices (const Triangulation<1>::face_iterator &, + FaceVertexNormals &) const { Assert (false, Boundary<1>::ExcFunctionNotUseful(1)); }; @@ -217,13 +217,15 @@ get_tangents_at_vertices (const Triangulation<1>::face_iterator &, template <> void StraightBoundary<2>:: -get_tangents_at_vertices (const Triangulation<2>::face_iterator &face, - Boundary<2>::FaceVertexTangents &face_vertex_tangents) const +get_normals_at_vertices (const Triangulation<2>::face_iterator &face, + Boundary<2>::FaceVertexNormals &face_vertex_normals) const { const unsigned int dim=2; const Tensor<1,dim> tangent = face->vertex(1) - face->vertex(0); for (unsigned int vertex=0; vertex::vertices_per_face; ++vertex) - face_vertex_tangents.tangents[vertex][0] = tangent; + // compute normals from tangent + face_vertex_normals[vertex] = Point(tangent[1], + -tangent[0]); }; #endif @@ -238,17 +240,29 @@ 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); + const unsigned int vertices_per_face = GeometryInfo<3>::vertices_per_face; - 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); + for (unsigned int vertex=0; vertex tangents[2] + = { face->vertex((vertex+1) % vertices_per_face) + - face->vertex(vertex), + face->vertex((vertex+vertices_per_face-1) % vertices_per_face) + - face->vertex(vertex) }; + + // then compute the normal by + // taking the cross + // product. since the normal is + // not required to be + // normalized, no problem here + cross_product (face_vertex_normals[vertex], + tangents[0], tangents[1]); + }; }; #endif 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 7605ef0ce6..db57d1d01a 100644 --- a/deal.II/deal.II/source/grid/tria_boundary_lib.cc +++ b/deal.II/deal.II/source/grid/tria_boundary_lib.cc @@ -219,8 +219,8 @@ HyperBallBoundary::get_intermediate_points_on_quad ( template <> void HyperBallBoundary<1>:: -get_tangents_at_vertices (const Triangulation<1>::face_iterator &, - FaceVertexTangents &) const +get_normals_at_vertices (const Triangulation<1>::face_iterator &, + FaceVertexNormals &) const { Assert (false, Boundary<1>::ExcFunctionNotUseful(1)); }; @@ -228,116 +228,16 @@ get_tangents_at_vertices (const Triangulation<1>::face_iterator &, #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::vertices_per_face; ++vertex) - { - const Tensor<1,dim> radius = face->vertex(vertex)-center; - face_vertex_tangents.tangents[vertex][0] = Point(radius[1], -radius[0]); - }; -}; - -#endif - - - -#if deal_II_dimension == 3 - -template <> +template void -HyperBallBoundary<3>:: -get_tangents_at_vertices (const Triangulation<3>::face_iterator &face, - Boundary<3>::FaceVertexTangents &face_vertex_tangents) const +HyperBallBoundary:: +get_normals_at_vertices (const Triangulation::face_iterator &face, + typename Boundary::FaceVertexNormals &face_vertex_normals) const { - const unsigned int dim=3; - // construct two tangential vectors - // at each of the vertices for (unsigned int vertex=0; vertex::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); - }; + face_vertex_normals[vertex] = face->vertex(vertex)-center; }; -#endif template @@ -472,8 +372,8 @@ get_intermediate_points_on_quad (const Triangulation<1>::quad_iterator &, template <> void HalfHyperBallBoundary<1>:: -get_tangents_at_vertices (const Triangulation<1>::face_iterator &, - Boundary<1>::FaceVertexTangents &) const +get_normals_at_vertices (const Triangulation<1>::face_iterator &, + Boundary<1>::FaceVertexNormals &) const { Assert (false, Boundary<1>::ExcFunctionNotUseful(1)); }; @@ -484,8 +384,8 @@ get_tangents_at_vertices (const Triangulation<1>::face_iterator &, template void HalfHyperBallBoundary:: -get_tangents_at_vertices (const typename Triangulation::face_iterator &face, - typename Boundary::FaceVertexTangents &face_vertex_tangents) const +get_normals_at_vertices (const typename Triangulation::face_iterator &face, + typename Boundary::FaceVertexNormals &face_vertex_normals) const { // check whether center of object is // at x==0, since then it belongs @@ -493,9 +393,9 @@ get_tangents_at_vertices (const typename Triangulation::face_iterator &face // boundary const Point quad_center = face->center(); if (quad_center(0) == center(0)) - StraightBoundary::get_tangents_at_vertices (face, face_vertex_tangents); + StraightBoundary::get_normals_at_vertices (face, face_vertex_normals); else - HyperBallBoundary::get_tangents_at_vertices (face, face_vertex_tangents); + HyperBallBoundary::get_normals_at_vertices (face, face_vertex_normals); }; diff --git a/deal.II/doc/news/2001/c-3-1.html b/deal.II/doc/news/2001/c-3-1.html index 10141bb268..756ed96c19 100644 --- a/deal.II/doc/news/2001/c-3-1.html +++ b/deal.II/doc/news/2001/c-3-1.html @@ -337,8 +337,8 @@ documentation, etc.
  • New: Boundary and derived classes now have a function get_tangents_at_vertices that returns a - description of the tangent space to the boundary at the + class="member">get_normals_at_vertices that returns a + multiple of the normal vector to the boundary at the vertices of a given face. This is used in the construction of C1 mappings of the boundary.