get_tangent_vector (const Point<spacedim> &x1,
const Point<spacedim> &x2) const override;
+ /**
+ * Return the normal vector to the given face at point p taking into account
+ * that quadrilateral faces of hexahedral cells in 3d may not be planar.
+ * In those cases, the face is assumed to have a geometry described by a
+ * bilinear function, and the normal vector is computed by embedding this
+ * bilinear form into a Cartesian space with a flat metric.
+ */
+ virtual
+ Tensor<1,spacedim>
+ normal_vector (const typename Triangulation<dim,spacedim>::face_iterator &face,
+ const Point<spacedim> &p) const override;
+
+ /**
+ * Compute the normal vectors to the boundary at each vertex of the
+ * given face taking into account that quadrilateral faces of hexahedral
+ * cells in 3d may not be planar. In those cases, the face is assumed to
+ * have a geometry described by a bilinear function, and the normal vector
+ * is computed by embedding this bilinear form into a Cartesian space with
+ * a flat metric.
+ */
+ virtual
+ void
+ get_normals_at_vertices (const typename Triangulation<dim,spacedim>::face_iterator &face,
+ typename Manifold<dim, spacedim>::FaceVertexNormals &face_vertex_normals) const override;
+
/**
* Return the periodicity of this Manifold.
*/
Boundary<1,3>::
get_intermediate_points_on_face (const Triangulation<1,3>::face_iterator &,
std::vector<Point<3> > &) const;
-template <>
-void
-StraightBoundary<1,1>::
-get_normals_at_vertices (const Triangulation<1,1>::face_iterator &,
- Boundary<1,1>::FaceVertexNormals &) const;
-template <>
-void
-StraightBoundary<2,2>::
-get_normals_at_vertices (const Triangulation<2,2>::face_iterator &face,
- Boundary<2,2>::FaceVertexNormals &face_vertex_normals) const;
-template <>
-void
-StraightBoundary<3,3>::
-get_normals_at_vertices (const Triangulation<3,3>::face_iterator &face,
- Boundary<3,3>::FaceVertexNormals &face_vertex_normals) const;
template <>
Point<3>
/* -------------------------- FlatManifold --------------------- */
+namespace internal
+{
+ namespace
+ {
+ Tensor<1,2>
+ normalized_alternating_product (const Tensor<1,2> (&basis_vectors)[1])
+ {
+ Tensor<1,2> tmp = cross_product_2d (basis_vectors[0]);
+ return tmp/tmp.norm();
+ }
+
+
+
+ Tensor<1,3>
+ normalized_alternating_product (const Tensor<1,3> ( &)[1])
+ {
+ // we get here from FlatManifold<2,3>::normal_vector, but
+ // the implementation below is bogus for this case anyway
+ // (see the assert at the beginning of that function).
+ Assert (false, ExcNotImplemented());
+ return Tensor<1,3>();
+ }
+
+
+
+ Tensor<1,3>
+ normalized_alternating_product (const Tensor<1,3> (&basis_vectors)[2])
+ {
+ Tensor<1,3> tmp = cross_product_3d (basis_vectors[0], basis_vectors[1]);
+ return tmp/tmp.norm();
+ }
+
+ }
+}
template <int dim, int spacedim>
FlatManifold<dim,spacedim>::FlatManifold (const Tensor<1,spacedim> &periodicity,
}
+
template <int dim, int spacedim>
Point<spacedim>
FlatManifold<dim, spacedim>::project_to_manifold
+template <>
+void
+FlatManifold<1>::
+get_normals_at_vertices (const Triangulation<1>::face_iterator &,
+ Manifold<1,1>::FaceVertexNormals &) const
+{
+ Assert (false, ExcImpossibleInDim(1));
+}
+
+
+
+template <>
+void
+FlatManifold<1,2>::
+get_normals_at_vertices (const Triangulation<1,2>::face_iterator &,
+ Manifold<1,2>::FaceVertexNormals &) const
+{
+ Assert (false, ExcNotImplemented());
+}
+
+
+
+template <>
+void
+FlatManifold<1,3>::
+get_normals_at_vertices (const Triangulation<1,3>::face_iterator &,
+ Manifold<1,3>::FaceVertexNormals &) const
+{
+ Assert (false, ExcNotImplemented());
+}
+
+
+
+template <>
+void
+FlatManifold<2>::
+get_normals_at_vertices (const Triangulation<2>::face_iterator &face,
+ Manifold<2,2>::FaceVertexNormals &face_vertex_normals) const
+{
+ const Tensor<1,2> tangent = face->vertex(1) - face->vertex(0);
+ for (unsigned int vertex=0; vertex<GeometryInfo<2>::vertices_per_face; ++vertex)
+ // compute normals from tangent
+ face_vertex_normals[vertex] = Point<2>(tangent[1],
+ -tangent[0]);
+}
+
+
+
+template <>
+void
+FlatManifold<2,3>::
+get_normals_at_vertices (const Triangulation<2,3>::face_iterator &face,
+ Manifold<2,3>::FaceVertexNormals &face_vertex_normals) const
+{
+ Assert(false, ExcNotImplemented());
+}
+
+
+
+template <>
+void
+FlatManifold<3>::
+get_normals_at_vertices (const Triangulation<3>::face_iterator &face,
+ Manifold<3,3>::FaceVertexNormals &face_vertex_normals) const
+{
+ const unsigned int vertices_per_face = GeometryInfo<3>::vertices_per_face;
+
+ static const unsigned int neighboring_vertices[4][2]=
+ { {1,2},{3,0},{0,3},{2,1}};
+ 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(neighboring_vertices[vertex][0])
+ - face->vertex(vertex),
+ face->vertex(neighboring_vertices[vertex][1])
+ - face->vertex(vertex)
+ };
+
+ // then compute the normal by taking the cross product. since the
+ // normal is not required to be normalized, no problem here
+ face_vertex_normals[vertex] = cross_product_3d(tangents[0], tangents[1]);
+ }
+}
+
+
+
+template <>
+Tensor<1,1>
+FlatManifold<1,1>::
+normal_vector (const Triangulation<1,1>::face_iterator &,
+ const Point<1> &) const
+{
+ Assert (false, ExcNotImplemented());
+ return Tensor<1,1>();
+}
+
+
+
+template <>
+Tensor<1,2>
+FlatManifold<1,2>::
+normal_vector (const Triangulation<1,2>::face_iterator &,
+ const Point<2> &) const
+{
+ Assert (false, ExcNotImplemented());
+ return Tensor<1,2>();
+}
+
+
+
+template <>
+Tensor<1,3>
+FlatManifold<1,3>::
+normal_vector (const Triangulation<1,3>::face_iterator &,
+ const Point<3> &) const
+{
+ Assert (false, ExcNotImplemented());
+ return Tensor<1,3>();
+}
+
+
+
+template <int dim, int spacedim>
+Tensor<1, spacedim>
+FlatManifold< dim, spacedim >::
+normal_vector (const typename Triangulation<dim, spacedim>::face_iterator &face,
+ const Point<spacedim> &p) const
+{
+ // I don't think the implementation below will work when dim!=spacedim;
+ // in fact, I believe that we don't even have enough information here,
+ // because we would need to know not only about the tangent vectors
+ // of the face, but also of the cell, to compute the normal vector.
+ // Someone will have to think about this some more.
+ Assert (dim == spacedim, ExcNotImplemented());
+
+ // in order to find out what the normal vector is, we first need to
+ // find the reference coordinates of the point p on the given face,
+ // or at least the reference coordinates of the closest point on the
+ // face
+ //
+ // in other words, we need to find a point xi so that f(xi)=||F(xi)-p||^2->min
+ // where F(xi) is the mapping. this algorithm is implemented in
+ // MappingQ1<dim,spacedim>::transform_real_to_unit_cell but only for cells,
+ // while we need it for faces here. it's also implemented in somewhat
+ // more generality there using the machinery of the MappingQ1 class
+ // while we really only need it for a specific case here
+ //
+ // in any case, the iteration we use here is a Gauss-Newton's iteration with
+ // xi^{n+1} = xi^n - H(xi^n)^{-1} J(xi^n)
+ // where
+ // J(xi) = (grad F(xi))^T (F(xi)-p)
+ // and
+ // H(xi) = [grad F(xi)]^T [grad F(xi)]
+ // In all this,
+ // F(xi) = sum_v vertex[v] phi_v(xi)
+ // We get the shape functions phi_v from an object of type FE_Q<dim-1>(1)
+
+ // we start with the point xi=1/2, xi=(1/2,1/2), ...
+ const unsigned int facedim = dim-1;
+
+ Point<facedim> xi;
+ for (unsigned int i=0; i<facedim; ++i)
+ xi[i] = 1./2;
+
+ const double eps = 1e-12;
+ Tensor<1,spacedim> grad_F[facedim];
+ unsigned int iteration = 0;
+ while (true)
+ {
+ Point<spacedim> F;
+ for (unsigned int v=0; v<GeometryInfo<facedim>::vertices_per_cell; ++v)
+ F += face->vertex(v) * GeometryInfo<facedim>::d_linear_shape_function(xi, v);
+
+ for (unsigned int i=0; i<facedim; ++i)
+ {
+ grad_F[i] = 0;
+ for (unsigned int v=0; v<GeometryInfo<facedim>::vertices_per_cell; ++v)
+ grad_F[i] += face->vertex(v) *
+ GeometryInfo<facedim>::d_linear_shape_function_gradient(xi, v)[i];
+ }
+
+ Tensor<1,facedim> J;
+ for (unsigned int i=0; i<facedim; ++i)
+ for (unsigned int j=0; j<spacedim; ++j)
+ J[i] += grad_F[i][j] * (F-p)[j];
+
+ Tensor<2,facedim> H;
+ for (unsigned int i=0; i<facedim; ++i)
+ for (unsigned int j=0; j<facedim; ++j)
+ for (unsigned int k=0; k<spacedim; ++k)
+ H[i][j] += grad_F[i][k] * grad_F[j][k];
+
+ const Tensor<1,facedim> delta_xi = -invert(H) * J;
+ xi += delta_xi;
+ ++iteration;
+
+ Assert (iteration<10,
+ ExcMessage("The Newton iteration to find the reference point "
+ "did not converge in 10 iterations. Do you have a "
+ "deformed cell? (See the glossary for a definition "
+ "of what a deformed cell is. You may want to output "
+ "the vertices of your cell."));
+
+ if (delta_xi.norm() < eps)
+ break;
+ }
+
+ // so now we have the reference coordinates xi of the point p.
+ // we then have to compute the normal vector, which we can do
+ // by taking the (normalize) alternating product of all the tangent
+ // vectors given by grad_F
+ return internal::normalized_alternating_product(grad_F);
+}
+
+
/* -------------------------- ChartManifold --------------------- */
template <int dim, int spacedim, int chartdim>
ChartManifold<dim,spacedim,chartdim>::ChartManifold (const Tensor<1,chartdim> &periodicity)
-template <>
-Tensor<1,1>
-StraightBoundary<1,1>::
-normal_vector (const Triangulation<1,1>::face_iterator &,
- const Point<1> &) const
-{
- Assert (false, ExcNotImplemented());
- return Tensor<1,1>();
-}
-
-
-template <>
-Tensor<1,2>
-StraightBoundary<1,2>::
-normal_vector (const Triangulation<1,2>::face_iterator &,
- const Point<2> &) const
-{
- Assert (false, ExcNotImplemented());
- return Tensor<1,2>();
-}
-
-
-template <>
-Tensor<1,3>
-StraightBoundary<1,3>::
-normal_vector (const Triangulation<1,3>::face_iterator &,
- const Point<3> &) const
-{
- Assert (false, ExcNotImplemented());
- return Tensor<1,3>();
-}
-
-
-namespace internal
-{
- namespace
- {
- /**
- * Compute the normalized cross product of a set of dim-1 basis
- * vectors.
- */
- Tensor<1,2>
- normalized_alternating_product (const Tensor<1,2> (&basis_vectors)[1])
- {
- Tensor<1,2> tmp = cross_product_2d (basis_vectors[0]);
- return tmp/tmp.norm();
- }
-
-
-
- Tensor<1,3>
- normalized_alternating_product (const Tensor<1,3> ( &)[1])
- {
- // we get here from StraightBoundary<2,3>::normal_vector, but
- // the implementation below is bogus for this case anyway
- // (see the assert at the beginning of that function).
- Assert (false, ExcNotImplemented());
- return Tensor<1,3>();
- }
-
-
-
- Tensor<1,3>
- normalized_alternating_product (const Tensor<1,3> (&basis_vectors)[2])
- {
- Tensor<1,3> tmp = cross_product_3d (basis_vectors[0], basis_vectors[1]);
- return tmp/tmp.norm();
- }
-
- }
-}
-
-
template <int dim, int spacedim>
Tensor<1,spacedim>
StraightBoundary<dim,spacedim>::
normal_vector (const typename Triangulation<dim,spacedim>::face_iterator &face,
- const Point<spacedim> &p) const
-{
- // I don't think the implementation below will work when dim!=spacedim;
- // in fact, I believe that we don't even have enough information here,
- // because we would need to know not only about the tangent vectors
- // of the face, but also of the cell, to compute the normal vector.
- // Someone will have to think about this some more.
- Assert (dim == spacedim, ExcNotImplemented());
-
- // in order to find out what the normal vector is, we first need to
- // find the reference coordinates of the point p on the given face,
- // or at least the reference coordinates of the closest point on the
- // face
- //
- // in other words, we need to find a point xi so that f(xi)=||F(xi)-p||^2->min
- // where F(xi) is the mapping. this algorithm is implemented in
- // MappingQ1<dim,spacedim>::transform_real_to_unit_cell but only for cells,
- // while we need it for faces here. it's also implemented in somewhat
- // more generality there using the machinery of the MappingQ1 class
- // while we really only need it for a specific case here
- //
- // in any case, the iteration we use here is a Gauss-Newton's iteration with
- // xi^{n+1} = xi^n - H(xi^n)^{-1} J(xi^n)
- // where
- // J(xi) = (grad F(xi))^T (F(xi)-p)
- // and
- // H(xi) = [grad F(xi)]^T [grad F(xi)]
- // In all this,
- // F(xi) = sum_v vertex[v] phi_v(xi)
- // We get the shape functions phi_v from an object of type FE_Q<dim-1>(1)
-
- // we start with the point xi=1/2, xi=(1/2,1/2), ...
- const unsigned int facedim = dim-1;
-
- Point<facedim> xi;
- for (unsigned int i=0; i<facedim; ++i)
- xi[i] = 1./2;
-
- const double eps = 1e-12;
- Tensor<1,spacedim> grad_F[facedim];
- unsigned int iteration = 0;
- while (true)
- {
- Point<spacedim> F;
- for (unsigned int v=0; v<GeometryInfo<facedim>::vertices_per_cell; ++v)
- F += face->vertex(v) * GeometryInfo<facedim>::d_linear_shape_function(xi, v);
-
- for (unsigned int i=0; i<facedim; ++i)
- {
- grad_F[i] = 0;
- for (unsigned int v=0; v<GeometryInfo<facedim>::vertices_per_cell; ++v)
- grad_F[i] += face->vertex(v) *
- GeometryInfo<facedim>::d_linear_shape_function_gradient(xi, v)[i];
- }
-
- Tensor<1,facedim> J;
- for (unsigned int i=0; i<facedim; ++i)
- for (unsigned int j=0; j<spacedim; ++j)
- J[i] += grad_F[i][j] * (F-p)[j];
-
- Tensor<2,facedim> H;
- for (unsigned int i=0; i<facedim; ++i)
- for (unsigned int j=0; j<facedim; ++j)
- for (unsigned int k=0; k<spacedim; ++k)
- H[i][j] += grad_F[i][k] * grad_F[j][k];
-
- const Tensor<1,facedim> delta_xi = -invert(H) * J;
- xi += delta_xi;
- ++iteration;
-
- Assert (iteration<10,
- ExcMessage("The Newton iteration to find the reference point "
- "did not converge in 10 iterations. Do you have a "
- "deformed cell? (See the glossary for a definition "
- "of what a deformed cell is. You may want to output "
- "the vertices of your cell."));
-
- if (delta_xi.norm() < eps)
- break;
- }
-
- // so now we have the reference coordinates xi of the point p.
- // we then have to compute the normal vector, which we can do
- // by taking the (normalize) alternating product of all the tangent
- // vectors given by grad_F
- return internal::normalized_alternating_product(grad_F);
-}
-
-
-
-template <>
-void
-StraightBoundary<1>::
-get_normals_at_vertices (const Triangulation<1>::face_iterator &,
- Boundary<1,1>::FaceVertexNormals &) const
-{
- Assert (false, ExcImpossibleInDim(1));
-}
-
-template <>
-void
-StraightBoundary<1,2>::
-get_normals_at_vertices (const Triangulation<1,2>::face_iterator &,
- Boundary<1,2>::FaceVertexNormals &) const
-{
- Assert (false, ExcNotImplemented());
-}
-
-
-template <>
-void
-StraightBoundary<1,3>::
-get_normals_at_vertices (const Triangulation<1,3>::face_iterator &,
- Boundary<1,3>::FaceVertexNormals &) const
+ const Point<spacedim> &p) const
{
- Assert (false, ExcNotImplemented());
+ return FlatManifold<dim, spacedim>::normal_vector(face, p);
}
-template <>
-void
-StraightBoundary<2>::
-get_normals_at_vertices (const Triangulation<2>::face_iterator &face,
- Boundary<2,2>::FaceVertexNormals &face_vertex_normals) const
-{
- const Tensor<1,2> tangent = face->vertex(1) - face->vertex(0);
- for (unsigned int vertex=0; vertex<GeometryInfo<2>::vertices_per_face; ++vertex)
- // compute normals from tangent
- face_vertex_normals[vertex] = Point<2>(tangent[1],
- -tangent[0]);
-}
-
-template <>
-void
-StraightBoundary<2,3>::
-get_normals_at_vertices (const Triangulation<2,3>::face_iterator &face,
- Boundary<2,3>::FaceVertexNormals &face_vertex_normals) const
-{
- const Tensor<1,3> tangent = face->vertex(1) - face->vertex(0);
- for (unsigned int vertex=0; vertex<GeometryInfo<2>::vertices_per_face; ++vertex)
- // compute normals from tangent
- face_vertex_normals[vertex] = Point<3>(tangent[1],
- -tangent[0],0);
- Assert(false, ExcNotImplemented());
-}
-
-
-
-
-template <>
+template <int dim, int spacedim>
void
-StraightBoundary<3>::
-get_normals_at_vertices (const Triangulation<3>::face_iterator &face,
- Boundary<3,3>::FaceVertexNormals &face_vertex_normals) const
+StraightBoundary<dim, spacedim>::
+get_normals_at_vertices (const typename Triangulation<dim, spacedim>::face_iterator &face,
+ typename Boundary<dim,spacedim>::FaceVertexNormals &face_vertex_normals) const
{
- const unsigned int vertices_per_face = GeometryInfo<3>::vertices_per_face;
-
- static const unsigned int neighboring_vertices[4][2]=
- { {1,2},{3,0},{0,3},{2,1}};
- 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(neighboring_vertices[vertex][0])
- - face->vertex(vertex),
- face->vertex(neighboring_vertices[vertex][1])
- - face->vertex(vertex)
- };
-
- // then compute the normal by taking the cross product. since the
- // normal is not required to be normalized, no problem here
- face_vertex_normals[vertex] = cross_product_3d(tangents[0], tangents[1]);
- };
+ FlatManifold<dim, spacedim>::get_normals_at_vertices(face, face_vertex_normals);
}