<< "impossible in " << arg1
<< "d or simply does not make any sense.");
+ /**
+ * This exception is raised if a functionality is not possible in the given
+ * combination of dimension and space-dimension.
+ *
+ * The constructor takes two <tt>int</tt>, denoting the dimension and the
+ * space dimension.
+ */
+ DeclException2 (ExcImpossibleInDimSpacedim,
+ int, int,
+ << "You are trying to execute functionality that is "
+ << "impossible in dimensions <" << arg1 << "," << arg2
+ << "> or simply does not make any sense.");
+
+
/**
* A number is zero, but it should not be here.
*/
{
public:
+ /**
+ * Type keeping information about the normals at the vertices of a face of a
+ * cell. Thus, there are <tt>GeometryInfo<dim>::vertices_per_face</tt>
+ * 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 type is not useful in 1d.
+ */
+ typedef Tensor<1,spacedim> FaceVertexNormals[GeometryInfo<dim>::vertices_per_face];
+
/**
* Destructor. Does nothing here, but needs to be declared virtual to make
const Point<spacedim> &x2) const;
/// @}
+
+ /**
+ * @name Computing normal vectors
+ */
+ /// @{
+
+ /**
+ * Return the normal vector to a face embedded in this manifold, at
+ * the point p. If p is not in fact on the surface, but only
+ * close-by, try to return something reasonable, for example the
+ * normal vector at the surface point closest to p. (The point p
+ * will in fact not normally lie on the actual surface, but rather
+ * be a quadrature point mapped by some polynomial mapping; the
+ * mapped surface, however, will not usually coincide with the
+ * actual surface.)
+ *
+ * The face iterator gives an indication which face this function is
+ * supposed to compute the normal vector for. This is useful if the
+ * boundary of the domain is composed of different nondifferential
+ * pieces (for example when using the StraightBoundary class to
+ * approximate a geometry that is completely described by the coarse
+ * mesh, with piecewise (bi-)linear components between the vertices,
+ * but where the boundary may have a kink at the vertices itself).
+ *
+ * @note The default implementation of this function computes the
+ * normal vector by taking the cross product between the tangent
+ * vectors from p to the most orthogonal and further non consecutive
+ * vertices of the face.
+ */
+ virtual
+ Tensor<1,spacedim>
+ normal_vector (const typename Triangulation<dim,spacedim>::face_iterator &face,
+ const Point<spacedim> &p) const;
+
+ /**
+ * Compute the normal vectors to the boundary at each vertex of the
+ * given face embedded in the Manifold. 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. The
+ * default implementation calls normal_vector() on each vertex.
+ *
+ * Note that when computing normal vectors at a vertex where the
+ * boundary is not differentiable, you have to make sure that you
+ * compute the one-sided limits, i.e. limit with respect to points
+ * inside the given face.
+ */
+ virtual
+ void
+ get_normals_at_vertices (const typename Triangulation<dim,spacedim>::face_iterator &face,
+ FaceVertexNormals &face_vertex_normals) const;
+
+ /// @}
};
{
public:
- /**
- * Type keeping information about the normals at the vertices of a face of a
- * cell. Thus, there are <tt>GeometryInfo<dim>::vertices_per_face</tt>
- * 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 type is not useful in 1d.
- */
- typedef Tensor<1,spacedim> FaceVertexNormals[GeometryInfo<dim>::vertices_per_face];
-
/**
* Destructor. Does nothing here, but needs to be declared to make it
* virtual.
get_intermediate_points_on_face (const typename Triangulation<dim,spacedim>::face_iterator &face,
std::vector<Point<spacedim> > &points) const;
- /**
- * Return the normal vector to the surface at the point p. If p is not in
- * fact on the surface, but only close-by, try to return something
- * reasonable, for example the normal vector at the surface point closest to
- * p. (The point p will in fact not normally lie on the actual surface, but
- * rather be a quadrature point mapped by some polynomial mapping; the
- * mapped surface, however, will not usually coincide with the actual
- * surface.)
- *
- * The face iterator gives an indication which face this function is
- * supposed to compute the normal vector for. This is useful if the
- * boundary of the domain is composed of different nondifferential pieces
- * (for example when using the StraightBoundary class to approximate a
- * geometry that is completely described by the coarse mesh, with piecewise
- * (bi-)linear components between the vertices, but where the boundary may
- * have a kink at the vertices itself).
- *
- * @note Implementations of this function should be able to assume that the
- * point p lies within or close to the face described by the first argument.
- * In turn, callers of this function should ensure that this is in fact the
- * case.
- */
- virtual
- Tensor<1,spacedim>
- normal_vector (const typename Triangulation<dim,spacedim>::face_iterator &face,
- const Point<spacedim> &p) const;
-
- /**
- * 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. The default
- * implementation is to throw an error, so you need not overload this
- * function in case you do not intend to use C1 mappings.
- *
- * Note that when computing normal vectors at a vertex where the boundary is
- * not differentiable, you have to make sure that you compute the one-sided
- * limits, i.e. limit with respect to points inside the given face.
- */
- virtual
- void
- get_normals_at_vertices (const typename Triangulation<dim,spacedim>::face_iterator &face,
- FaceVertexNormals &face_vertex_normals) const;
-
/**
* Given a candidate point and a line segment characterized by the iterator,
* return a point that lies on the surface described by this object. This
}
+
+template <>
+Tensor<1,2>
+Manifold<2, 2>::
+normal_vector (const typename Triangulation<2, 2>::face_iterator &face,
+ const Point<2> &p) const
+{
+ const int spacedim=2, dim=2;
+
+ Tensor<1,spacedim> tangent = ((p-face->vertex(0)).norm_square() > (p-face->vertex(1)).norm_square() ?
+ get_tangent_vector(p, face->vertex(0)) :
+ get_tangent_vector(p, face->vertex(1)));
+ return cross_product_2d(tangent);
+}
+
+template<>
+Tensor<1,3>
+Manifold<3, 3>::
+normal_vector (const typename Triangulation<3, 3>::face_iterator &face,
+ const Point<3> &p) const
+{
+ const int spacedim=3, dim=3;
+ Tensor<1,spacedim> t1,t2;
+
+ // Take the difference between p and all four vertices
+ Tensor<1,spacedim> dp[4];
+ double dpns[4];
+ int min_index=-1;
+ double min_distance = 0;
+ for (unsigned int i=0; i<4; ++i)
+ {
+ dp[i] = p-face->vertex(i);
+ dpns[i] = dp[i].norm_square();
+ min_index = (min_index == -1 ? (int)i : dpns[i] < min_distance ? i : min_index);
+ min_distance = dpns[min_index];
+ }
+ // Verify we have a valid vertex index
+ AssertIndexRange(min_index, 4);
+
+ // now figure out which vertices are better to compute tangent vectors
+ // we split the cell in 4 quadrants, and use v1/v2 for the first quadrant
+ // (the one at the corner with vertex 0)
+ if ((p-face->center()).norm_square() < min_distance)
+ {
+ // we are close to the face center: pick two consecutive vertices,
+ // but not the closest one. We make sure the direction is always
+ // the same.
+ if (min_index < 2)
+ {
+ t1 = get_tangent_vector(p, face->vertex(3));
+ t2 = get_tangent_vector(p, face->vertex(2));
+ }
+ else
+ {
+ t1 = get_tangent_vector(p, face->vertex(0));
+ t2 = get_tangent_vector(p, face->vertex(1));
+ }
+ }
+ else
+ {
+ switch (min_index)
+ {
+ case 0:
+ {
+ t1 = get_tangent_vector(p, face->vertex(1));
+ t2 = get_tangent_vector(p, face->vertex(2));
+ break;
+ }
+ case 1:
+ {
+ t1 = get_tangent_vector(p, face->vertex(3));
+ t2 = get_tangent_vector(p, face->vertex(0));
+ break;
+ }
+ case 2:
+ {
+ t1 = get_tangent_vector(p, face->vertex(0));
+ t2 = get_tangent_vector(p, face->vertex(3));
+ break;
+ }
+ case 3:
+ {
+ t1 = get_tangent_vector(p, face->vertex(2));
+ t2 = get_tangent_vector(p, face->vertex(1));
+ break;
+ }
+ default:
+ Assert(false, ExcInternalError());
+ break;
+ }
+ }
+ return cross_product_3d(t1,t2);
+}
+
+
+template <int dim, int spacedim>
+Tensor<1,spacedim>
+Manifold<dim, spacedim>::
+normal_vector (const typename Triangulation<dim, spacedim>::face_iterator &face,
+ const Point<spacedim> &p) const
+{
+ Assert(false, ExcPureFunctionCalled());
+ return Tensor<1,spacedim>();
+}
+
+
+template <int dim, int spacedim>
+void
+Manifold<dim, spacedim>::
+get_normals_at_vertices (const typename Triangulation<dim, spacedim>::face_iterator &face,
+ FaceVertexNormals &n) const
+{
+ for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_face; ++v)
+ n[v] = normal_vector(face, face->vertex(v));
+}
+
+
+
+
template <int dim, int spacedim>
Point<spacedim>
Manifold<dim, spacedim>::
{
#if deal_II_dimension <= deal_II_space_dimension
template class SphericalManifold<deal_II_dimension, deal_II_space_dimension>;
-#endif
-#if deal_II_dimension <= deal_II_space_dimension
template class CylindricalManifold<deal_II_dimension, deal_II_space_dimension>;
template class FunctionManifold<deal_II_dimension, deal_II_space_dimension, 1>;
template class FunctionManifold<deal_II_dimension, deal_II_space_dimension, 2>;
}
-
-
-template <int dim, int spacedim>
-Tensor<1,spacedim>
-Boundary<dim, spacedim>::
-normal_vector (const typename Triangulation<dim, spacedim>::face_iterator &,
- const Point<spacedim> &) const
-{
- Assert (false, ExcPureFunctionCalled());
- return Tensor<1,spacedim>();
-}
-
-
-
-template <int dim, int spacedim>
-void
-Boundary<dim, spacedim>::
-get_normals_at_vertices (const typename Triangulation<dim, spacedim>::face_iterator &,
- FaceVertexNormals &) const
-{
- Assert (false, ExcPureFunctionCalled());
-}
-
-
-
template <int dim, int spacedim>
Point<spacedim>
Boundary<dim, spacedim>::
projects / dealii.git / commitdiff