From 00d7da81adb052f0a85fb6e981bc531ba974852c Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Tue, 2 Jan 2018 19:56:32 +0100 Subject: [PATCH] Move implementation of normal_vector and get_normal_at_vertices from StraightBoundary to FlatManifold --- include/deal.II/grid/manifold.h | 25 +++ include/deal.II/grid/tria_boundary.h | 15 -- source/grid/manifold.cc | 252 +++++++++++++++++++++++++++ source/grid/tria_boundary.cc | 248 +------------------------- 4 files changed, 284 insertions(+), 256 deletions(-) diff --git a/include/deal.II/grid/manifold.h b/include/deal.II/grid/manifold.h index e34149a72d..6f57df46e4 100644 --- a/include/deal.II/grid/manifold.h +++ b/include/deal.II/grid/manifold.h @@ -773,6 +773,31 @@ public: get_tangent_vector (const Point &x1, const Point &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::face_iterator &face, + const Point &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::face_iterator &face, + typename Manifold::FaceVertexNormals &face_vertex_normals) const override; + /** * Return the periodicity of this Manifold. */ diff --git a/include/deal.II/grid/tria_boundary.h b/include/deal.II/grid/tria_boundary.h index d13e70246b..fe7053c73c 100644 --- a/include/deal.II/grid/tria_boundary.h +++ b/include/deal.II/grid/tria_boundary.h @@ -392,21 +392,6 @@ void Boundary<1,3>:: get_intermediate_points_on_face (const Triangulation<1,3>::face_iterator &, std::vector > &) 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> diff --git a/source/grid/manifold.cc b/source/grid/manifold.cc index 13eb7ee81a..b8bbbd2b3c 100644 --- a/source/grid/manifold.cc +++ b/source/grid/manifold.cc @@ -476,6 +476,40 @@ Manifold::get_tangent_vector(const Point &x1, /* -------------------------- 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 FlatManifold::FlatManifold (const Tensor<1,spacedim> &periodicity, @@ -617,6 +651,7 @@ get_new_points (const ArrayView> &surrounding_points, } + template Point FlatManifold::project_to_manifold @@ -662,6 +697,223 @@ FlatManifold::get_tangent_vector (const Point &x1, +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::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 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 +Tensor<1, spacedim> +FlatManifold< dim, spacedim >:: +normal_vector (const typename Triangulation::face_iterator &face, + const Point &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::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(1) + + // we start with the point xi=1/2, xi=(1/2,1/2), ... + const unsigned int facedim = dim-1; + + Point xi; + for (unsigned int i=0; i grad_F[facedim]; + unsigned int iteration = 0; + while (true) + { + Point F; + for (unsigned int v=0; v::vertices_per_cell; ++v) + F += face->vertex(v) * GeometryInfo::d_linear_shape_function(xi, v); + + for (unsigned int i=0; i::vertices_per_cell; ++v) + grad_F[i] += face->vertex(v) * + GeometryInfo::d_linear_shape_function_gradient(xi, v)[i]; + } + + Tensor<1,facedim> J; + for (unsigned int i=0; i H; + for (unsigned int i=0; i 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 ChartManifold::ChartManifold (const Tensor<1,chartdim> &periodicity) diff --git a/source/grid/tria_boundary.cc b/source/grid/tria_boundary.cc index 250dd5dbae..7fb26fe4ab 100644 --- a/source/grid/tria_boundary.cc +++ b/source/grid/tria_boundary.cc @@ -418,258 +418,24 @@ get_intermediate_points_on_quad (const Triangulation<2,3>::quad_iterator &quad, -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 Tensor<1,spacedim> StraightBoundary:: normal_vector (const typename Triangulation::face_iterator &face, - const Point &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::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(1) - - // we start with the point xi=1/2, xi=(1/2,1/2), ... - const unsigned int facedim = dim-1; - - Point xi; - for (unsigned int i=0; i grad_F[facedim]; - unsigned int iteration = 0; - while (true) - { - Point F; - for (unsigned int v=0; v::vertices_per_cell; ++v) - F += face->vertex(v) * GeometryInfo::d_linear_shape_function(xi, v); - - for (unsigned int i=0; i::vertices_per_cell; ++v) - grad_F[i] += face->vertex(v) * - GeometryInfo::d_linear_shape_function_gradient(xi, v)[i]; - } - - Tensor<1,facedim> J; - for (unsigned int i=0; i H; - for (unsigned int i=0; i 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 &p) const { - Assert (false, ExcNotImplemented()); + return FlatManifold::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::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::vertices_per_face; ++vertex) - // compute normals from tangent - face_vertex_normals[vertex] = Point<3>(tangent[1], - -tangent[0],0); - Assert(false, ExcNotImplemented()); -} - - - - -template <> +template void -StraightBoundary<3>:: -get_normals_at_vertices (const Triangulation<3>::face_iterator &face, - Boundary<3,3>::FaceVertexNormals &face_vertex_normals) const +StraightBoundary:: +get_normals_at_vertices (const typename Triangulation::face_iterator &face, + typename Boundary::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 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::get_normals_at_vertices(face, face_vertex_normals); } -- 2.39.5