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<dim>::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<unsigned int>(-1));
+
+ /**
+ * Array of tangent vectors,
+ * as described above.
+ */
+ Tensor<1,dim> tangents[GeometryInfo<dim>::vertices_per_face][n_tangents];
+ };
+
+
/**
* Destructor. Does nothing here, but
* needs to be declared to make it
get_intermediate_points_on_quad (const typename Triangulation<dim>::quad_iterator &quad,
typename std::vector<Point<dim> > &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<dim>::face_iterator &face,
+ FaceVertexTangents &face_vertex_tangents) const;
+
+
/**
* Exception.
*/
* @author Wolfgang Bangerth, 1998, Ralf Hartmann, 2001
*/
template <int dim>
-class StraightBoundary : public Boundary<dim> {
+class StraightBoundary : public Boundary<dim>
+{
public:
/**
* Let the new point be the arithmetic
get_intermediate_points_on_quad (const typename Triangulation<dim>::quad_iterator &quad,
typename std::vector<Point<dim> > &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<dim>::face_iterator &face,
+ typename Boundary<dim>::FaceVertexTangents &face_vertex_tangents) const;
};
get_intermediate_points_on_quad (const typename Triangulation<dim>::quad_iterator &quad,
typename std::vector<Point<dim> > &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<dim>::face_iterator &face,
+ typename Boundary<dim>::FaceVertexTangents &face_vertex_tangents) const;
+
/**
* Return the center of the ball.
*/
virtual void
get_intermediate_points_on_quad (const typename Triangulation<dim>::quad_iterator &quad,
typename std::vector<Point<dim> > &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<dim>::face_iterator &face,
+ typename Boundary<dim>::FaceVertexTangents &face_vertex_tangents) const;
};
#include <cmath>
+
+/* -------------------------- StraightBoundary --------------------- */
+
+
template <int dim>
Boundary<dim>::~Boundary ()
{};
+
template <int dim>
Point<dim>
Boundary<dim>::get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &) const
};
+
+template <int dim>
+void
+Boundary<dim>::
+get_intermediate_points_on_line (const typename Triangulation<dim>::line_iterator &,
+ std::vector<Point<dim> > &) const
+{
+ Assert (false, ExcPureVirtualFunctionCalled());
+};
+
+
+
template <int dim>
void
-Boundary<dim>::get_intermediate_points_on_line (
- const typename Triangulation<dim>::line_iterator &,
- std::vector<Point<dim> > &) const
+Boundary<dim>::
+get_intermediate_points_on_quad (const typename Triangulation<dim>::quad_iterator &,
+ std::vector<Point<dim> > &) const
{
Assert (false, ExcPureVirtualFunctionCalled());
};
+
template <int dim>
void
-Boundary<dim>::get_intermediate_points_on_quad (
- const typename Triangulation<dim>::quad_iterator &,
- std::vector<Point<dim> > &) const
+Boundary<dim>::
+get_tangents_at_vertices (const typename Triangulation<dim>::face_iterator &,
+ FaceVertexTangents &) const
{
Assert (false, ExcPureVirtualFunctionCalled());
};
+
+/* -------------------------- StraightBoundary --------------------- */
+
+
template <int dim>
Point<dim>
StraightBoundary<dim>::get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const
template <int dim>
Point<dim>
-StraightBoundary<dim>::get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const
+StraightBoundary<dim>::
+get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const
{
return (quad->vertex(0) + quad->vertex(1) +
quad->vertex(2) + quad->vertex(3) +
template <int dim>
void
-StraightBoundary<dim>::get_intermediate_points_on_line (
- const typename Triangulation<dim>::line_iterator &,
- typename std::vector<Point<dim> > &) const
+StraightBoundary<dim>::
+get_intermediate_points_on_line (const typename Triangulation<dim>::line_iterator &,
+ typename std::vector<Point<dim> > &) const
{
Assert(false, typename Boundary<dim>::ExcFunctionNotUseful(dim));
}
#else
+
template <int dim>
void
-StraightBoundary<dim>::get_intermediate_points_on_line (
- const typename Triangulation<dim>::line_iterator &line,
- typename std::vector<Point<dim> > &points) const
+StraightBoundary<dim>::
+get_intermediate_points_on_line (const typename Triangulation<dim>::line_iterator &line,
+ typename std::vector<Point<dim> > &points) const
{
const unsigned int n=points.size();
Assert(n>0, ExcInternalError());
template <int dim>
void
-StraightBoundary<dim>::get_intermediate_points_on_quad (
- const typename Triangulation<dim>::quad_iterator &,
- typename std::vector<Point<dim> > &) const
+StraightBoundary<dim>::
+get_intermediate_points_on_quad (const typename Triangulation<dim>::quad_iterator &,
+ typename std::vector<Point<dim> > &) const
{
Assert(false, typename Boundary<dim>::ExcFunctionNotUseful(dim));
}
template <int dim>
void
-StraightBoundary<dim>::get_intermediate_points_on_quad (
- const typename Triangulation<dim>::quad_iterator &quad,
- typename std::vector<Point<dim> > &points) const
+StraightBoundary<dim>::
+get_intermediate_points_on_quad (const typename Triangulation<dim>::quad_iterator &quad,
+ typename std::vector<Point<dim> > &points) const
{
const unsigned int n=points.size(),
m=static_cast<unsigned int>(sqrt(n));
#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<GeometryInfo<dim>::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<deal_II_dimension>;
template class StraightBoundary<deal_II_dimension>;
}
+
+#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<GeometryInfo<dim>::vertices_per_face; ++vertex)
+ {
+ const Tensor<1,dim> radius = face->vertex(vertex)-center;
+ face_vertex_tangents.tangents[vertex][0] = Point<dim>(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<GeometryInfo<dim>::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 <int dim>
Point<dim>
HyperBallBoundary<dim>::get_center () const
// boundary
const Point<dim> quad_center = quad->center();
if (quad_center(0) == center(0))
- return StraightBoundary<dim>::get_intermediate_points_on_quad (quad, points);
+ StraightBoundary<dim>::get_intermediate_points_on_quad (quad, points);
else
- return HyperBallBoundary<dim>::get_intermediate_points_on_quad (quad, points);
+ HyperBallBoundary<dim>::get_intermediate_points_on_quad (quad, points);
};
template <>
void
HalfHyperBallBoundary<1>::
-get_intermediate_points_on_quad (const typename Triangulation<1>::quad_iterator &,
- typename std::vector<Point<1> > &) const
+get_intermediate_points_on_quad (const Triangulation<1>::quad_iterator &,
+ std::vector<Point<1> > &) const
{
Assert (false, ExcInternalError());
};
#endif
+
+template <int dim>
+void
+HyperBallBoundary<dim>::
+get_tangents_at_vertices (const typename Triangulation<dim>::face_iterator &face,
+ typename Boundary<dim>::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<dim> quad_center = quad->center();
+ if (quad_center(0) == center(0))
+ StraightBoundary<dim>::get_tangents_at_vertices (face, face_vertex_tangents);
+ else
+ HyperBallBoundary<dim>::get_intermediate_points_on_quad (face, face_vertex_tangents);
+};
+
+
/* ---------------------------------------------------------------------- */
<br>
(WB 2001/01/30)
</p>
-
</ol>
<ol>
<li> <p>
+ New: <code class="class">Boundary</code> and derived classes
+ now have a function <code
+ class="member">get_tangents_at_vertices</code> 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
+ C<sup>1</sup> mappings of the boundary.
+ <br>
+ (WB 2001/03/23)
+ </p>
+ <li> <p>
New: The new function <code
class="class">PersistentTriangulation<dim></code>::<code
class="member">restore(unsigned int)</code> allows to restore