From df29b5725891c72c1084228fd96313757218c7cd Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Fri, 23 Mar 2001 14:48:42 +0000 Subject: [PATCH] Boundary::get_tangents_at_vertices to generate a description of the tangent space to a boundary object at each vertex of a face. git-svn-id: https://svn.dealii.org/trunk@4278 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/grid/tria_boundary.h | 95 ++++++++++- .../deal.II/include/grid/tria_boundary_lib.h | 29 ++++ deal.II/deal.II/source/grid/tria_boundary.cc | 119 +++++++++++--- .../deal.II/source/grid/tria_boundary_lib.cc | 154 +++++++++++++++++- deal.II/doc/news/2001/c-3-1.html | 11 +- 5 files changed, 383 insertions(+), 25 deletions(-) diff --git a/deal.II/deal.II/include/grid/tria_boundary.h b/deal.II/deal.II/include/grid/tria_boundary.h index 650bb9a680..9eb1af8c9e 100644 --- a/deal.II/deal.II/include/grid/tria_boundary.h +++ b/deal.II/deal.II/include/grid/tria_boundary.h @@ -76,6 +76,48 @@ template 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::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(-1)); + + /** + * Array of tangent vectors, + * as described above. + */ + Tensor<1,dim> tangents[GeometryInfo::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::quad_iterator &quad, 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. + * + * 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::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 -class StraightBoundary : public Boundary { +class StraightBoundary : public Boundary +{ public: /** * Let the new point be the arithmetic @@ -269,6 +349,19 @@ class StraightBoundary : public Boundary { 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. + * + * Refer to the general + * documentation of this class + * and the documentation of the + * base class. + */ + virtual void + get_tangents_at_vertices (const typename Triangulation::face_iterator &face, + typename Boundary::FaceVertexTangents &face_vertex_tangents) 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 3ccb5d45fa..a972d61bd3 100644 --- a/deal.II/deal.II/include/grid/tria_boundary_lib.h +++ b/deal.II/deal.II/include/grid/tria_boundary_lib.h @@ -84,6 +84,20 @@ class HyperBallBoundary : public StraightBoundary 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. + * + * Refer to the general + * documentation of this class + * and the documentation of the + * base class. + */ + virtual void + get_tangents_at_vertices (const typename Triangulation::face_iterator &face, + typename Boundary::FaceVertexTangents &face_vertex_tangents) const; + /** * Return the center of the ball. */ @@ -186,6 +200,21 @@ class HalfHyperBallBoundary : public HyperBallBoundary virtual void 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. + * + * Refer to the general + * documentation of this class + * and the documentation of the + * base class. + */ + virtual void + get_tangents_at_vertices (const typename Triangulation::face_iterator &face, + typename Boundary::FaceVertexTangents &face_vertex_tangents) const; }; diff --git a/deal.II/deal.II/source/grid/tria_boundary.cc b/deal.II/deal.II/source/grid/tria_boundary.cc index 5ef3b01668..85c8e67fc9 100644 --- a/deal.II/deal.II/source/grid/tria_boundary.cc +++ b/deal.II/deal.II/source/grid/tria_boundary.cc @@ -19,11 +19,16 @@ #include + +/* -------------------------- StraightBoundary --------------------- */ + + template Boundary::~Boundary () {}; + template Point Boundary::get_new_point_on_quad (const typename Triangulation::quad_iterator &) const @@ -33,26 +38,43 @@ Boundary::get_new_point_on_quad (const typename Triangulation::quad_it }; + +template +void +Boundary:: +get_intermediate_points_on_line (const typename Triangulation::line_iterator &, + std::vector > &) const +{ + Assert (false, ExcPureVirtualFunctionCalled()); +}; + + + template void -Boundary::get_intermediate_points_on_line ( - const typename Triangulation::line_iterator &, - std::vector > &) const +Boundary:: +get_intermediate_points_on_quad (const typename Triangulation::quad_iterator &, + std::vector > &) const { Assert (false, ExcPureVirtualFunctionCalled()); }; + template void -Boundary::get_intermediate_points_on_quad ( - const typename Triangulation::quad_iterator &, - std::vector > &) const +Boundary:: +get_tangents_at_vertices (const typename Triangulation::face_iterator &, + FaceVertexTangents &) const { Assert (false, ExcPureVirtualFunctionCalled()); }; + +/* -------------------------- StraightBoundary --------------------- */ + + template Point StraightBoundary::get_new_point_on_line (const typename Triangulation::line_iterator &line) const @@ -77,7 +99,8 @@ StraightBoundary::get_new_point_on_quad (const typename Triangulation: template Point -StraightBoundary::get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const +StraightBoundary:: +get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const { return (quad->vertex(0) + quad->vertex(1) + quad->vertex(2) + quad->vertex(3) + @@ -94,9 +117,9 @@ StraightBoundary::get_new_point_on_quad (const typename Triangulation: template void -StraightBoundary::get_intermediate_points_on_line ( - const typename Triangulation::line_iterator &, - typename std::vector > &) const +StraightBoundary:: +get_intermediate_points_on_line (const typename Triangulation::line_iterator &, + typename std::vector > &) const { Assert(false, typename Boundary::ExcFunctionNotUseful(dim)); } @@ -104,11 +127,12 @@ StraightBoundary::get_intermediate_points_on_line ( #else + template void -StraightBoundary::get_intermediate_points_on_line ( - const typename Triangulation::line_iterator &line, - typename std::vector > &points) const +StraightBoundary:: +get_intermediate_points_on_line (const typename Triangulation::line_iterator &line, + typename std::vector > &points) const { const unsigned int n=points.size(); Assert(n>0, ExcInternalError()); @@ -131,9 +155,9 @@ StraightBoundary::get_intermediate_points_on_line ( template void -StraightBoundary::get_intermediate_points_on_quad ( - const typename Triangulation::quad_iterator &, - typename std::vector > &) const +StraightBoundary:: +get_intermediate_points_on_quad (const typename Triangulation::quad_iterator &, + typename std::vector > &) const { Assert(false, typename Boundary::ExcFunctionNotUseful(dim)); } @@ -142,9 +166,9 @@ StraightBoundary::get_intermediate_points_on_quad ( template void -StraightBoundary::get_intermediate_points_on_quad ( - const typename Triangulation::quad_iterator &quad, - typename std::vector > &points) const +StraightBoundary:: +get_intermediate_points_on_quad (const typename Triangulation::quad_iterator &quad, + typename std::vector > &points) const { const unsigned int n=points.size(), m=static_cast(sqrt(n)); @@ -173,6 +197,63 @@ StraightBoundary::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::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; template class StraightBoundary; 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 aa98b78406..447d2b4b96 100644 --- a/deal.II/deal.II/source/grid/tria_boundary_lib.cc +++ b/deal.II/deal.II/source/grid/tria_boundary_lib.cc @@ -213,6 +213,133 @@ HyperBallBoundary::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::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 <> +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::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 Point HyperBallBoundary::get_center () const @@ -319,9 +446,9 @@ get_intermediate_points_on_quad (const typename Triangulation::quad_iterato // boundary const Point quad_center = quad->center(); if (quad_center(0) == center(0)) - return StraightBoundary::get_intermediate_points_on_quad (quad, points); + StraightBoundary::get_intermediate_points_on_quad (quad, points); else - return HyperBallBoundary::get_intermediate_points_on_quad (quad, points); + HyperBallBoundary::get_intermediate_points_on_quad (quad, points); }; @@ -331,8 +458,8 @@ get_intermediate_points_on_quad (const typename Triangulation::quad_iterato template <> void HalfHyperBallBoundary<1>:: -get_intermediate_points_on_quad (const typename Triangulation<1>::quad_iterator &, - typename std::vector > &) const +get_intermediate_points_on_quad (const Triangulation<1>::quad_iterator &, + std::vector > &) const { Assert (false, ExcInternalError()); }; @@ -340,6 +467,25 @@ get_intermediate_points_on_quad (const typename Triangulation<1>::quad_iterator #endif + +template +void +HyperBallBoundary:: +get_tangents_at_vertices (const typename Triangulation::face_iterator &face, + typename Boundary::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 quad_center = quad->center(); + if (quad_center(0) == center(0)) + StraightBoundary::get_tangents_at_vertices (face, face_vertex_tangents); + else + HyperBallBoundary::get_intermediate_points_on_quad (face, face_vertex_tangents); +}; + + /* ---------------------------------------------------------------------- */ diff --git a/deal.II/doc/news/2001/c-3-1.html b/deal.II/doc/news/2001/c-3-1.html index 88fc805124..e3675419f1 100644 --- a/deal.II/doc/news/2001/c-3-1.html +++ b/deal.II/doc/news/2001/c-3-1.html @@ -272,7 +272,6 @@ documentation, etc.
(WB 2001/01/30)

- @@ -282,7 +281,17 @@ documentation, etc.
  1. + 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 + vertices of a given face. This is used in the construction of + C1 mappings of the boundary. +
    + (WB 2001/03/23) +

    +
  2. New: The new function PersistentTriangulation<dim>::restore(unsigned int) allows to restore -- 2.39.5