From 80c6b9bf4f1dccd9ddf56b30defa4111e0379c07 Mon Sep 17 00:00:00 2001 From: heltai Date: Sat, 18 Jan 2014 17:08:50 +0000 Subject: [PATCH] Restored approximation of the normal at arbitrary points for StraightBoundary git-svn-id: https://svn.dealii.org/branches/branch_manifold_id@32248 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/include/deal.II/grid/manifold.h | 18 +- deal.II/include/deal.II/grid/tria_boundary.h | 32 +++- .../deal.II/numerics/vector_tools.templates.h | 5 +- deal.II/source/grid/manifold.cc | 8 + deal.II/source/grid/tria_boundary.cc | 170 ++++++++++++++++++ 5 files changed, 228 insertions(+), 5 deletions(-) diff --git a/deal.II/include/deal.II/grid/manifold.h b/deal.II/include/deal.II/grid/manifold.h index e26dd7e4e5..ed85ca3678 100644 --- a/deal.II/include/deal.II/grid/manifold.h +++ b/deal.II/include/deal.II/grid/manifold.h @@ -103,14 +103,28 @@ public: * Given a vector of points, return the normals to the Manifold in * those points. Notice that the Manifold may or may not be a * codimension one manifold. If it is not, this function will throw - * an Error. Input arguments must be of the same size. The default - * implementation calls for each point the function + * an Error. @points and @normals must be of the same size. The + * default implementation calls for each point the function * normal_vector. Derived classes can overload that function. */ virtual void get_normals_at_points(const std::vector > &points, std::vector > &normals) const; + /** + * Compute the normal to the manifold at the given point. The + * default implementation of this function just calls + * normal_vector. The difference from this function and that one, + * is that one can use this function when only an approximation of + * the normal is required, for example when computing normals at + * quadrature points. Note that this function only makes sense + * when the manifold is a codimension one manifold. + */ + virtual + Point + normal_vector(const std::vector > &vertices, + const Point &point) const; + /** * Compute the normal to the manifold at the given point. The * default implementation of this function throws an diff --git a/deal.II/include/deal.II/grid/tria_boundary.h b/deal.II/include/deal.II/grid/tria_boundary.h index 6c1fc8ee08..72b31d4db4 100644 --- a/deal.II/include/deal.II/grid/tria_boundary.h +++ b/deal.II/include/deal.II/grid/tria_boundary.h @@ -155,11 +155,39 @@ public: virtual void get_normals_at_points(const std::vector > &points, std::vector > &normals) const; - - + /** + * Given a point and a vector of vertices, compute an approximation + * of the normal to the manifold. + */ + virtual + Point + normal_vector(const std::vector > &vertices, + const Point &point) const; }; +/* -------------- declaration of explicit specializations ------------- */ + +#ifndef DOXYGEN + +template <> +Point<1> +StraightBoundary<1,1>::normal_vector(const std::vector > &vertices, + const Point<1> &point) const; + +template <> +Point<2> +StraightBoundary<1,2>::normal_vector(const std::vector > &vertices, + const Point<2> &point) const; + +template <> +Point<3> +StraightBoundary<1,3>::normal_vector(const std::vector > &vertices, + const Point<3> &point) const; + + +#endif + DEAL_II_NAMESPACE_CLOSE #endif diff --git a/deal.II/include/deal.II/numerics/vector_tools.templates.h b/deal.II/include/deal.II/numerics/vector_tools.templates.h index 8d2ed0861c..cc1b81c42e 100644 --- a/deal.II/include/deal.II/numerics/vector_tools.templates.h +++ b/deal.II/include/deal.II/numerics/vector_tools.templates.h @@ -4419,9 +4419,12 @@ namespace VectorTools // sign of the normal vector provided by the boundary // if they should point in different directions. this is the // case in tests/deal.II/no_flux_11. + std::vector > vertices(GeometryInfo::vertices_per_face); + for(unsigned int v=0; v::vertices_per_face; ++v) + vertices[v] = cell->face(face_no)->vertex(v); Point normal_vector = (cell->face(face_no)->get_boundary() - .normal_vector (fe_values.quadrature_point(i))); + .normal_vector (vertices, fe_values.quadrature_point(i))); if (normal_vector * fe_values.normal_vector(i) < 0) normal_vector *= -1; Assert (std::fabs(normal_vector.norm() - 1) < 1e-14, diff --git a/deal.II/source/grid/manifold.cc b/deal.II/source/grid/manifold.cc index ffd6b6dd05..1d598217a8 100644 --- a/deal.II/source/grid/manifold.cc +++ b/deal.II/source/grid/manifold.cc @@ -50,6 +50,14 @@ Manifold::get_normals_at_points(const std::vector > &p normals[i] = normal_vector(points[i]); } +template +Point +Manifold::normal_vector(const std::vector > &, + const Point &point) const +{ + return normal_vector(point); +} + template Point Manifold::normal_vector(const Point &point) const diff --git a/deal.II/source/grid/tria_boundary.cc b/deal.II/source/grid/tria_boundary.cc index e220300f62..d33e1ac01f 100644 --- a/deal.II/source/grid/tria_boundary.cc +++ b/deal.II/source/grid/tria_boundary.cc @@ -16,6 +16,8 @@ #include #include +#include + DEAL_II_NAMESPACE_OPEN /* -------------------------- Boundary --------------------- */ @@ -28,6 +30,45 @@ Boundary::~Boundary () /* -------------------------- StraightBoundary --------------------- */ +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 (tmp, 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 (tmp, basis_vectors[0], basis_vectors[1]); + return tmp/tmp.norm(); + } + } +} + + + template StraightBoundary::StraightBoundary () {} @@ -91,6 +132,135 @@ StraightBoundary::get_normals_at_points(const std::vector +Point<1> +StraightBoundary<1,1>:: +normal_vector (const std::vector > &vertices, + const Point<1> &) const +{ + Assert (false, ExcNotImplemented()); + return Point<1>(); +} + + + +template <> +Point<2> +StraightBoundary<1,2>:: +normal_vector (const std::vector > &vertices, + const Point<2> &) const +{ + Assert (false, ExcNotImplemented()); + return Point<2>(); +} + + + +template <> +Point<3> +StraightBoundary<1,3>:: +normal_vector (const std::vector > &vertices, + const Point<3> &) const +{ + Assert (false, ExcNotImplemented()); + return Point<3>(); +} + + +template +Point +StraightBoundary::normal_vector(const std::vector > &vertices, + const Point &p) const { + AssertDimension(vertices.size(), GeometryInfo::vertices_per_face); + + // 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) + // + // 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 linear_fe(1); + + const double eps = 1e-12; + Tensor<1,spacedim> grad_F[facedim]; + while (true) + { + Point F; + for (unsigned int v=0; v::vertices_per_cell; ++v) + F += vertices[v] * linear_fe.shape_value(v, xi); + + for (unsigned int i=0; i::vertices_per_cell; ++v) + grad_F[i] += vertices[v] * linear_fe.shape_grad(v, xi)[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; + + 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); +} + + // explicit instantiations #include "tria_boundary.inst" -- 2.39.5