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.
-
+ 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)
+
+ -
New: The new function PersistentTriangulation<dim>
::restore(unsigned int)
allows to restore
--
2.39.5