* @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 <int dim>
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<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.
+ * 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<unsigned int>(-1));
-
- /**
- * Array of tangent vectors,
- * as described above.
- */
- Tensor<1,dim> tangents[GeometryInfo<dim>::vertices_per_face][n_tangents];
- };
-
+ typedef Tensor<1,dim> FaceVertexNormals[GeometryInfo<dim>::vertices_per_face];
/**
* Destructor. Does nothing here, but
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<dim>
get_new_point_on_line (const typename Triangulation<dim>::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<dim>
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.
+ * 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.
*
* 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<dim>::face_iterator &face,
- FaceVertexTangents &face_vertex_tangents) const;
+ get_normals_at_vertices (const typename Triangulation<dim>::face_iterator &face,
+ FaceVertexNormals &face_vertex_normals) const;
/**
/**
- * Specialisation of @ref{Boundary}<dim>, 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}<dim>, 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 <int dim>
class StraightBoundary : public Boundary<dim>
{
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<dim>
get_new_point_on_line (const typename Triangulation<dim>::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<dim>
get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const;
typename std::vector<Point<dim> > &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
* base class.
*/
virtual void
- get_tangents_at_vertices (const typename Triangulation<dim>::face_iterator &face,
- typename Boundary<dim>::FaceVertexTangents &face_vertex_tangents) const;
+ get_normals_at_vertices (const typename Triangulation<dim>::face_iterator &face,
+ typename Boundary<dim>::FaceVertexNormals &face_vertex_normals) const;
};
typename std::vector<Point<dim> > &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
* base class.
*/
virtual void
- get_tangents_at_vertices (const typename Triangulation<dim>::face_iterator &face,
- typename Boundary<dim>::FaceVertexTangents &face_vertex_tangents) const;
+ get_normals_at_vertices (const typename Triangulation<dim>::face_iterator &face,
+ typename Boundary<dim>::FaceVertexNormals &face_vertex_normals) const;
/**
* Return the center of the ball.
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.
+ * Compute the normals to the
+ * boundary at the vertices of
+ * the given face.
*
* Refer to the general
* documentation of this class
* base class.
*/
virtual void
- get_tangents_at_vertices (const typename Triangulation<dim>::face_iterator &face,
- typename Boundary<dim>::FaceVertexTangents &face_vertex_tangents) const;
+ get_normals_at_vertices (const typename Triangulation<dim>::face_iterator &face,
+ typename Boundary<dim>::FaceVertexNormals &face_vertex_normals) const;
};
if (line->at_boundary())
{
- // first get the tangential
+ // first get the normal
// vectors at the two
// vertices of this line
// from the boundary
const Boundary<dim> &boundary
= line->get_triangulation().get_boundary(line->boundary_indicator());
- Boundary<dim>::FaceVertexTangents face_vertex_tangents;
- boundary.get_tangents_at_vertices (line, face_vertex_tangents);
+ Boundary<dim>::FaceVertexNormals face_vertex_normals;
+ boundary.get_normals_at_vertices (line, face_vertex_normals);
// then transform them into
// interpolation points for
// 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
// 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);
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
template <int dim>
void
Boundary<dim>::
-get_tangents_at_vertices (const typename Triangulation<dim>::face_iterator &,
- FaceVertexTangents &) const
+get_normals_at_vertices (const typename Triangulation<dim>::face_iterator &,
+ FaceVertexNormals &) const
{
Assert (false, ExcPureVirtualFunctionCalled());
};
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));
};
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<GeometryInfo<dim>::vertices_per_face; ++vertex)
- face_vertex_tangents.tangents[vertex][0] = tangent;
+ // compute normals from tangent
+ face_vertex_normals[vertex] = Point<dim>(tangent[1],
+ -tangent[0]);
};
#endif
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<vertices_per_face; ++vertex)
+ {
+ // first define the two tangent
+ // vectors at the vertex by
+ // using the two lines
+ // radiating away from this
+ // vertex
+ const Tensor<1,3> 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
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));
};
#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 <>
+template <int dim>
void
-HyperBallBoundary<3>::
-get_tangents_at_vertices (const Triangulation<3>::face_iterator &face,
- Boundary<3>::FaceVertexTangents &face_vertex_tangents) const
+HyperBallBoundary<dim>::
+get_normals_at_vertices (const Triangulation<dim>::face_iterator &face,
+ typename Boundary<dim>::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<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);
- };
+ face_vertex_normals[vertex] = face->vertex(vertex)-center;
};
-#endif
template <int dim>
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));
};
template <int dim>
void
HalfHyperBallBoundary<dim>::
-get_tangents_at_vertices (const typename Triangulation<dim>::face_iterator &face,
- typename Boundary<dim>::FaceVertexTangents &face_vertex_tangents) const
+get_normals_at_vertices (const typename Triangulation<dim>::face_iterator &face,
+ typename Boundary<dim>::FaceVertexNormals &face_vertex_normals) const
{
// check whether center of object is
// at x==0, since then it belongs
// boundary
const Point<dim> quad_center = face->center();
if (quad_center(0) == center(0))
- StraightBoundary<dim>::get_tangents_at_vertices (face, face_vertex_tangents);
+ StraightBoundary<dim>::get_normals_at_vertices (face, face_vertex_normals);
else
- HyperBallBoundary<dim>::get_tangents_at_vertices (face, face_vertex_tangents);
+ HyperBallBoundary<dim>::get_normals_at_vertices (face, face_vertex_normals);
};
<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
+ class="member">get_normals_at_vertices</code> 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
C<sup>1</sup> mappings of the boundary.
<br>