From 361d64f09b51f8fd5b94093306ef1041cbff8b25 Mon Sep 17 00:00:00 2001 From: bangerth Date: Thu, 25 Aug 2011 15:51:59 +0000 Subject: [PATCH] Add Boundary::normal_vector. Implement it. Test it. git-svn-id: https://svn.dealii.org/trunk@24197 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/include/deal.II/grid/tria_boundary.h | 81 ++++- .../include/deal.II/grid/tria_boundary_lib.h | 35 +- deal.II/source/grid/tria_boundary.cc | 164 ++++++++++ deal.II/source/grid/tria_boundary_lib.cc | 14 + tests/deal.II/normal_vector_01.cc | 92 ++++++ tests/deal.II/normal_vector_01/cmp/generic | 4 + tests/deal.II/normal_vector_01_2d.cc | 92 ++++++ tests/deal.II/normal_vector_01_2d/cmp/generic | 4 + tests/deal.II/normal_vector_02.cc | 65 ++++ tests/deal.II/normal_vector_02/cmp/generic | 2 + tests/deal.II/normal_vector_02_2d.cc | 65 ++++ tests/deal.II/normal_vector_02_2d/cmp/generic | 2 + tests/deal.II/normal_vector_03.cc | 114 +++++++ tests/deal.II/normal_vector_03/cmp/generic | 304 ++++++++++++++++++ tests/deal.II/normal_vector_03_2d.cc | 114 +++++++ tests/deal.II/normal_vector_03_2d/cmp/generic | 204 ++++++++++++ 16 files changed, 1340 insertions(+), 16 deletions(-) create mode 100644 tests/deal.II/normal_vector_01.cc create mode 100644 tests/deal.II/normal_vector_01/cmp/generic create mode 100644 tests/deal.II/normal_vector_01_2d.cc create mode 100644 tests/deal.II/normal_vector_01_2d/cmp/generic create mode 100644 tests/deal.II/normal_vector_02.cc create mode 100644 tests/deal.II/normal_vector_02/cmp/generic create mode 100644 tests/deal.II/normal_vector_02_2d.cc create mode 100644 tests/deal.II/normal_vector_02_2d/cmp/generic create mode 100644 tests/deal.II/normal_vector_03.cc create mode 100644 tests/deal.II/normal_vector_03/cmp/generic create mode 100644 tests/deal.II/normal_vector_03_2d.cc create mode 100644 tests/deal.II/normal_vector_03_2d/cmp/generic diff --git a/deal.II/include/deal.II/grid/tria_boundary.h b/deal.II/include/deal.II/grid/tria_boundary.h index 046fdf3066..1b2871ca11 100644 --- a/deal.II/include/deal.II/grid/tria_boundary.h +++ b/deal.II/include/deal.II/grid/tria_boundary.h @@ -122,7 +122,8 @@ class Boundary : public Subscriptor * lines therefore is also on the * boundary). */ - virtual Point + virtual + Point get_new_point_on_line (const typename Triangulation::line_iterator &line) const = 0; /** @@ -151,7 +152,8 @@ class Boundary : public Subscriptor * default implementation throws * an error in any case, however. */ - virtual Point + virtual + Point get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const; /** @@ -193,7 +195,8 @@ class Boundary : public Subscriptor * default implementation throws * an error in any case, however. */ - virtual void + virtual + void get_intermediate_points_on_line (const typename Triangulation::line_iterator &line, std::vector > &points) const; @@ -227,7 +230,8 @@ class Boundary : public Subscriptor * default implementation throws * an error in any case, however. */ - virtual void + virtual + void get_intermediate_points_on_quad (const typename Triangulation::quad_iterator &quad, std::vector > &points) const; @@ -248,6 +252,46 @@ class Boundary : public Subscriptor get_intermediate_points_on_face (const typename Triangulation::face_iterator &face, std::vector > &points) const; + /** + * Return the normal vector to the surface + * at the point p. If p is not in fact + * on the surface, but only closeby, + * 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::face_iterator &face, + const Point &p) const; + /** * Compute the normal vectors to * the boundary at each vertex of @@ -276,7 +320,8 @@ class Boundary : public Subscriptor * with respect to points inside * the given face. */ - virtual void + virtual + void get_normals_at_vertices (const typename Triangulation::face_iterator &face, FaceVertexNormals &face_vertex_normals) const; @@ -414,7 +459,8 @@ class StraightBoundary : public Boundary * base class for more * information. */ - virtual Point + virtual + Point get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const; /** @@ -429,7 +475,8 @@ class StraightBoundary : public Boundary * and the documentation of the * base class. */ - virtual void + virtual + void get_intermediate_points_on_line (const typename Triangulation::line_iterator &line, std::vector > &points) const; @@ -445,11 +492,26 @@ class StraightBoundary : public Boundary * and the documentation of the * base class. */ - virtual void + virtual + void get_intermediate_points_on_quad (const typename Triangulation::quad_iterator &quad, std::vector > &points) const; /** + * Implementation of the function + * declared in the base class. + * + * Refer to the general + * documentation of this class + * and the documentation of the + * base class. + */ + virtual + Tensor<1,spacedim> + normal_vector (const typename Triangulation::face_iterator &face, + const Point &p) const; + + /** * Compute the normals to the * boundary at the vertices of * the given face. @@ -459,7 +521,8 @@ class StraightBoundary : public Boundary * and the documentation of the * base class. */ - virtual void + virtual + void get_normals_at_vertices (const typename Triangulation::face_iterator &face, typename Boundary::FaceVertexNormals &face_vertex_normals) const; diff --git a/deal.II/include/deal.II/grid/tria_boundary_lib.h b/deal.II/include/deal.II/grid/tria_boundary_lib.h index 178f5983c0..eb5f68f804 100644 --- a/deal.II/include/deal.II/grid/tria_boundary_lib.h +++ b/deal.II/include/deal.II/grid/tria_boundary_lib.h @@ -388,7 +388,8 @@ class HyperBallBoundary : public StraightBoundary * this class and the documentation of the * base class. */ - virtual Point + virtual + Point get_new_point_on_line (const typename Triangulation::line_iterator &line) const; /** @@ -396,7 +397,8 @@ class HyperBallBoundary : public StraightBoundary * this class and the documentation of the * base class. */ - virtual Point + virtual + Point get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const; /** @@ -408,7 +410,8 @@ class HyperBallBoundary : public StraightBoundary * Calls * @p get_intermediate_points_between_points. */ - virtual void + virtual + void get_intermediate_points_on_line (const typename Triangulation::line_iterator &line, std::vector > &points) const; @@ -421,11 +424,26 @@ class HyperBallBoundary : public StraightBoundary * Only implemented for dim=3 * and for points.size()==1. */ - virtual void + virtual + void get_intermediate_points_on_quad (const typename Triangulation::quad_iterator &quad, std::vector > &points) const; /** + * Implementation of the function + * declared in the base class. + * + * Refer to the general + * documentation of this class + * and the documentation of the + * base class. + */ + virtual + Tensor<1,spacedim> + normal_vector (const typename Triangulation::face_iterator &face, + const Point &p) const; + + /** * Compute the normals to the * boundary at the vertices of * the given face. @@ -435,19 +453,22 @@ class HyperBallBoundary : public StraightBoundary * and the documentation of the * base class. */ - virtual void + virtual + void get_normals_at_vertices (const typename Triangulation::face_iterator &face, typename Boundary::FaceVertexNormals &face_vertex_normals) const; /** * Return the center of the ball. */ - Point get_center () const; + Point + get_center () const; /** * Return the radius of the ball. */ - double get_radius () const; + double + get_radius () const; /** * Exception. Thrown by the diff --git a/deal.II/source/grid/tria_boundary.cc b/deal.II/source/grid/tria_boundary.cc index 6b4a049a11..91b7f1e213 100644 --- a/deal.II/source/grid/tria_boundary.cc +++ b/deal.II/source/grid/tria_boundary.cc @@ -17,6 +17,7 @@ #include #include #include +#include #include DEAL_II_NAMESPACE_OPEN @@ -150,6 +151,18 @@ get_intermediate_points_on_face (const Triangulation<1,2>::face_iterator &, +template +Tensor<1,spacedim> +Boundary:: +normal_vector (const typename Triangulation::face_iterator &, + const Point &) const +{ + Assert (false, ExcPureFunctionCalled()); + return Tensor<1,spacedim>(); +} + + + template void Boundary:: @@ -472,6 +485,157 @@ get_intermediate_points_on_quad (const Triangulation<2,3>::quad_iterator &quad, +template <> +Tensor<1,1> +StraightBoundary<1,1>:: +normal_vector (const typename 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 typename Triangulation<1,2>::face_iterator &, + const Point<2> &) const +{ + Assert (false, ExcNotImplemented()); + return Tensor<1,2>(); +} + + + +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> (&basis_vectors)[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 +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 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 += face->vertex(v) * linear_fe.shape_value(v, xi); + + for (unsigned int i=0; i::vertices_per_cell; ++v) + grad_F[i] += face->vertex(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); +} + + + template <> void StraightBoundary<1>:: diff --git a/deal.II/source/grid/tria_boundary_lib.cc b/deal.II/source/grid/tria_boundary_lib.cc index 774cc39dc0..2cf45eac64 100644 --- a/deal.II/source/grid/tria_boundary_lib.cc +++ b/deal.II/source/grid/tria_boundary_lib.cc @@ -515,6 +515,7 @@ get_new_point_on_quad (const Triangulation<1,1>::quad_iterator &) const return Point<1>(); } + template <> Point<2> HyperBallBoundary<1,2>:: @@ -576,6 +577,7 @@ HyperBallBoundary::get_intermediate_points_on_line ( } + template void HyperBallBoundary::get_intermediate_points_between_points ( @@ -732,6 +734,18 @@ HyperBallBoundary::get_intermediate_points_on_quad ( +template +Tensor<1,spacedim> +HyperBallBoundary:: +normal_vector (const typename Triangulation::face_iterator &, + const Point &p) const +{ + const Tensor<1,spacedim> unnormalized_normal = p-center; + return unnormalized_normal/unnormalized_normal.norm(); +} + + + template <> void HyperBallBoundary<1>:: diff --git a/tests/deal.II/normal_vector_01.cc b/tests/deal.II/normal_vector_01.cc new file mode 100644 index 0000000000..6d8925e4de --- /dev/null +++ b/tests/deal.II/normal_vector_01.cc @@ -0,0 +1,92 @@ +//---------------------------- normal_vector_01.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2005, 2008, 2011 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- normal_vector_01.cc --------------------------- + + +// test that at the vertices, Boundary::normal_vector returns the same as +// Boundary::get_normals_at_vertices once the latter vectors are normalized + + + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include + +#include +#include + + + +void create_triangulation(const unsigned int case_no, + Triangulation<3> &tria) +{ + switch (case_no) + { + case 0: + GridGenerator::hyper_cube(tria, 1., 3.); + break; + case 1: + { + GridGenerator::hyper_cube(tria, 1., 3.); + Point<3> &v0=tria.begin_active()->vertex(0); + v0 = Point<3> (0,-0.5,-1); + Point<3> &v1=tria.begin_active()->vertex(1); + v1 = Point<3> (1.25, 0.25, 0.25); + break; + } + default: + Assert(false, ExcNotImplemented()); + }; +} + + + +int main () +{ + std::ofstream logfile ("normal_vector_01/output"); + deallog << std::setprecision (3); + deallog << std::fixed; + deallog.attach(logfile); + deallog.depth_console (0); + + Triangulation<3> tria; + StraightBoundary<3> boundary; + Boundary<3>::FaceVertexNormals normals; + for (unsigned int case_no=0; case_no<2; ++case_no) + { + deallog << "Case" << case_no << std::endl; + create_triangulation(case_no, tria); + const Triangulation<3>::active_cell_iterator cell=tria.begin_active(); + Triangulation<3>::face_iterator face; + for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) + { + face=cell->face(face_no); + boundary.get_normals_at_vertices(face, normals); + for (unsigned int v=0; v::vertices_per_face; ++v) + Assert ((boundary.normal_vector (face, + face->vertex(v)) + - + normals[v] / normals[v].norm()).norm() + < + 1e-12, + ExcInternalError()); + } + tria.clear(); + } + + deallog << "OK" << std::endl; +} diff --git a/tests/deal.II/normal_vector_01/cmp/generic b/tests/deal.II/normal_vector_01/cmp/generic new file mode 100644 index 0000000000..dfea1647fe --- /dev/null +++ b/tests/deal.II/normal_vector_01/cmp/generic @@ -0,0 +1,4 @@ + +DEAL::Case0 +DEAL::Case1 +DEAL::OK diff --git a/tests/deal.II/normal_vector_01_2d.cc b/tests/deal.II/normal_vector_01_2d.cc new file mode 100644 index 0000000000..39b0a7f744 --- /dev/null +++ b/tests/deal.II/normal_vector_01_2d.cc @@ -0,0 +1,92 @@ +//---------------------------- normal_vector_01_2d.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2005, 2008, 2011 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- normal_vector_01_2d.cc --------------------------- + + +// test that at the vertices, Boundary::normal_vector returns the same as +// Boundary::get_normals_at_vertices once the latter vectors are normalized + + + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include + +#include +#include + + + +void create_triangulation(const unsigned int case_no, + Triangulation<2> &tria) +{ + switch (case_no) + { + case 0: + GridGenerator::hyper_cube(tria, 1., 3.); + break; + case 1: + { + GridGenerator::hyper_cube(tria, 1., 3.); + Point<2> &v0=tria.begin_active()->vertex(0); + v0 = Point<2> (-0.5,-1); + Point<2> &v1=tria.begin_active()->vertex(1); + v1 = Point<2> (0.25, 0.25); + break; + } + default: + Assert(false, ExcNotImplemented()); + }; +} + + + +int main () +{ + std::ofstream logfile ("normal_vector_01_2d/output"); + deallog << std::setprecision (3); + deallog << std::fixed; + deallog.attach(logfile); + deallog.depth_console (0); + + Triangulation<2> tria; + StraightBoundary<2> boundary; + Boundary<2>::FaceVertexNormals normals; + for (unsigned int case_no=0; case_no<2; ++case_no) + { + deallog << "Case" << case_no << std::endl; + create_triangulation(case_no, tria); + const Triangulation<2>::active_cell_iterator cell=tria.begin_active(); + Triangulation<2>::face_iterator face; + for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) + { + face=cell->face(face_no); + boundary.get_normals_at_vertices(face, normals); + for (unsigned int v=0; v::vertices_per_face; ++v) + Assert ((boundary.normal_vector (face, + face->vertex(v)) + - + normals[v] / normals[v].norm()).norm() + < + 1e-12, + ExcInternalError()); + } + tria.clear(); + } + + deallog << "OK" << std::endl; +} diff --git a/tests/deal.II/normal_vector_01_2d/cmp/generic b/tests/deal.II/normal_vector_01_2d/cmp/generic new file mode 100644 index 0000000000..dfea1647fe --- /dev/null +++ b/tests/deal.II/normal_vector_01_2d/cmp/generic @@ -0,0 +1,4 @@ + +DEAL::Case0 +DEAL::Case1 +DEAL::OK diff --git a/tests/deal.II/normal_vector_02.cc b/tests/deal.II/normal_vector_02.cc new file mode 100644 index 0000000000..15c9f26946 --- /dev/null +++ b/tests/deal.II/normal_vector_02.cc @@ -0,0 +1,65 @@ +//---------------------------- normal_vector_02.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2005, 2008, 2011 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- normal_vector_02.cc --------------------------- + + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + + + + +int main () +{ + std::ofstream logfile ("normal_vector_02/output"); + deallog << std::setprecision (3); + deallog << std::fixed; + deallog.attach(logfile); + deallog.depth_console (0); + + HyperBallBoundary<3> boundary (Point<3>(1,0,0)); + + Triangulation<3> tria; + Boundary<3>::FaceVertexNormals normals; + + GridGenerator::hyper_ball (tria, Point<3>(1,0,0), 3); + + Triangulation<3>::active_cell_iterator cell=tria.begin_active(); + for (; cell!=tria.end(); ++cell) + for (unsigned int face_no=0; + face_no::faces_per_cell; ++face_no) + if (cell->at_boundary(face_no)) + { + Triangulation<3>::face_iterator face = cell->face(face_no); + boundary.get_normals_at_vertices(face, normals); + for (unsigned int v=0; v::vertices_per_face; ++v) + Assert ((boundary.normal_vector (face, + face->vertex(v)) + - + normals[v] / normals[v].norm()).norm() + < + 1e-12, + ExcInternalError()); + } + + deallog << "OK" << std::endl; +} diff --git a/tests/deal.II/normal_vector_02/cmp/generic b/tests/deal.II/normal_vector_02/cmp/generic new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/deal.II/normal_vector_02/cmp/generic @@ -0,0 +1,2 @@ + +DEAL::OK diff --git a/tests/deal.II/normal_vector_02_2d.cc b/tests/deal.II/normal_vector_02_2d.cc new file mode 100644 index 0000000000..6366127155 --- /dev/null +++ b/tests/deal.II/normal_vector_02_2d.cc @@ -0,0 +1,65 @@ +//---------------------------- normal_vector_02_2d.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2005, 2008, 2011 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- normal_vector_02_2d.cc --------------------------- + + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + + + + +int main () +{ + std::ofstream logfile ("normal_vector_02_2d/output"); + deallog << std::setprecision (3); + deallog << std::fixed; + deallog.attach(logfile); + deallog.depth_console (0); + + HyperBallBoundary<2> boundary (Point<2>(1,0)); + + Triangulation<2> tria; + Boundary<2>::FaceVertexNormals normals; + + GridGenerator::hyper_ball (tria, Point<2>(1,0), 3); + + Triangulation<2>::active_cell_iterator cell=tria.begin_active(); + for (; cell!=tria.end(); ++cell) + for (unsigned int face_no=0; + face_no::faces_per_cell; ++face_no) + if (cell->at_boundary(face_no)) + { + Triangulation<2>::face_iterator face = cell->face(face_no); + boundary.get_normals_at_vertices(face, normals); + for (unsigned int v=0; v::vertices_per_face; ++v) + Assert ((boundary.normal_vector (face, + face->vertex(v)) + - + normals[v] / normals[v].norm()).norm() + < + 1e-12, + ExcInternalError()); + } + + deallog << "OK" << std::endl; +} diff --git a/tests/deal.II/normal_vector_02_2d/cmp/generic b/tests/deal.II/normal_vector_02_2d/cmp/generic new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/deal.II/normal_vector_02_2d/cmp/generic @@ -0,0 +1,2 @@ + +DEAL::OK diff --git a/tests/deal.II/normal_vector_03.cc b/tests/deal.II/normal_vector_03.cc new file mode 100644 index 0000000000..00f04aa049 --- /dev/null +++ b/tests/deal.II/normal_vector_03.cc @@ -0,0 +1,114 @@ +//---------------------------- normal_vector_03.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2005, 2008, 2011 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- normal_vector_03.cc --------------------------- + + +// like _01 but instead for random points on the face. note that the normal +// vector on a bilinearly mapped face is also a bilinear function + + + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + + + +void create_triangulation(const unsigned int case_no, + Triangulation<3> &tria) +{ + switch (case_no) + { + case 0: + GridGenerator::hyper_cube(tria, 1., 3.); + break; + case 1: + { + GridGenerator::hyper_cube(tria, 1., 3.); + Point<3> &v0=tria.begin_active()->vertex(0); + v0 = Point<3> (0,-0.5,-1); + Point<3> &v1=tria.begin_active()->vertex(1); + v1 = Point<3> (1.25, 0.25, 0.25); + break; + } + default: + Assert(false, ExcNotImplemented()); + }; +} + + + +int main () +{ + std::ofstream logfile ("normal_vector_03/output"); + deallog << std::setprecision (3); + deallog << std::fixed; + deallog.attach(logfile); + deallog.depth_console (0); + + FE_Q<2> linear_interpolator(1); + + Triangulation<3> tria; + StraightBoundary<3> boundary; + Boundary<3>::FaceVertexNormals normals; + for (unsigned int case_no=0; case_no<2; ++case_no) + { + deallog << "Case" << case_no << std::endl; + create_triangulation(case_no, tria); + const Triangulation<3>::active_cell_iterator cell=tria.begin_active(); + Triangulation<3>::face_iterator face; + for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) + { + face=cell->face(face_no); + boundary.get_normals_at_vertices(face, normals); + + for (double xi=0; xi<=1; xi+=0.234) + for (double eta=0; eta<=1; eta+=0.234) + { + Point<3> p; + Tensor<1,3> normal; + + for (unsigned int v=0; v::vertices_per_face; ++v) + { + p += face->vertex(v) * linear_interpolator.shape_value(v,Point<2>(xi,eta)); + normal += normals[v] * + linear_interpolator.shape_value(v,Point<2>(xi,eta)); + } + normal /= normal.norm(); + + deallog << "p=" << p + << ", n=" << boundary.normal_vector (face, p) + << std::endl; + + Assert ((boundary.normal_vector (face, p) + - + normal).norm() + < + 1e-10, + ExcInternalError()); + } + } + + tria.clear(); + } + + deallog << "OK" << std::endl; +} diff --git a/tests/deal.II/normal_vector_03/cmp/generic b/tests/deal.II/normal_vector_03/cmp/generic new file mode 100644 index 0000000000..37c91d231f --- /dev/null +++ b/tests/deal.II/normal_vector_03/cmp/generic @@ -0,0 +1,304 @@ + +DEAL::Case0 +DEAL::p=1.000 1.000 1.000, n=1.000 0.000 0.000 +DEAL::p=1.000 1.000 1.468, n=1.000 0.000 0.000 +DEAL::p=1.000 1.000 1.936, n=1.000 0.000 0.000 +DEAL::p=1.000 1.000 2.404, n=1.000 0.000 0.000 +DEAL::p=1.000 1.000 2.872, n=1.000 0.000 0.000 +DEAL::p=1.000 1.468 1.000, n=1.000 0.000 0.000 +DEAL::p=1.000 1.468 1.468, n=1.000 0.000 0.000 +DEAL::p=1.000 1.468 1.936, n=1.000 0.000 0.000 +DEAL::p=1.000 1.468 2.404, n=1.000 0.000 0.000 +DEAL::p=1.000 1.468 2.872, n=1.000 0.000 0.000 +DEAL::p=1.000 1.936 1.000, n=1.000 0.000 0.000 +DEAL::p=1.000 1.936 1.468, n=1.000 0.000 0.000 +DEAL::p=1.000 1.936 1.936, n=1.000 0.000 0.000 +DEAL::p=1.000 1.936 2.404, n=1.000 0.000 0.000 +DEAL::p=1.000 1.936 2.872, n=1.000 0.000 0.000 +DEAL::p=1.000 2.404 1.000, n=1.000 0.000 0.000 +DEAL::p=1.000 2.404 1.468, n=1.000 0.000 0.000 +DEAL::p=1.000 2.404 1.936, n=1.000 0.000 0.000 +DEAL::p=1.000 2.404 2.404, n=1.000 0.000 0.000 +DEAL::p=1.000 2.404 2.872, n=1.000 0.000 0.000 +DEAL::p=1.000 2.872 1.000, n=1.000 0.000 0.000 +DEAL::p=1.000 2.872 1.468, n=1.000 0.000 0.000 +DEAL::p=1.000 2.872 1.936, n=1.000 0.000 0.000 +DEAL::p=1.000 2.872 2.404, n=1.000 0.000 0.000 +DEAL::p=1.000 2.872 2.872, n=1.000 0.000 0.000 +DEAL::p=3.000 1.000 1.000, n=1.000 0.000 0.000 +DEAL::p=3.000 1.000 1.468, n=1.000 0.000 0.000 +DEAL::p=3.000 1.000 1.936, n=1.000 0.000 0.000 +DEAL::p=3.000 1.000 2.404, n=1.000 0.000 0.000 +DEAL::p=3.000 1.000 2.872, n=1.000 0.000 0.000 +DEAL::p=3.000 1.468 1.000, n=1.000 0.000 0.000 +DEAL::p=3.000 1.468 1.468, n=1.000 0.000 0.000 +DEAL::p=3.000 1.468 1.936, n=1.000 0.000 0.000 +DEAL::p=3.000 1.468 2.404, n=1.000 0.000 0.000 +DEAL::p=3.000 1.468 2.872, n=1.000 0.000 0.000 +DEAL::p=3.000 1.936 1.000, n=1.000 0.000 0.000 +DEAL::p=3.000 1.936 1.468, n=1.000 0.000 0.000 +DEAL::p=3.000 1.936 1.936, n=1.000 0.000 0.000 +DEAL::p=3.000 1.936 2.404, n=1.000 0.000 0.000 +DEAL::p=3.000 1.936 2.872, n=1.000 0.000 0.000 +DEAL::p=3.000 2.404 1.000, n=1.000 0.000 0.000 +DEAL::p=3.000 2.404 1.468, n=1.000 0.000 0.000 +DEAL::p=3.000 2.404 1.936, n=1.000 0.000 0.000 +DEAL::p=3.000 2.404 2.404, n=1.000 0.000 0.000 +DEAL::p=3.000 2.404 2.872, n=1.000 0.000 0.000 +DEAL::p=3.000 2.872 1.000, n=1.000 0.000 0.000 +DEAL::p=3.000 2.872 1.468, n=1.000 0.000 0.000 +DEAL::p=3.000 2.872 1.936, n=1.000 0.000 0.000 +DEAL::p=3.000 2.872 2.404, n=1.000 0.000 0.000 +DEAL::p=3.000 2.872 2.872, n=1.000 0.000 0.000 +DEAL::p=1.000 1.000 1.000, n=0.000 1.000 0.000 +DEAL::p=1.468 1.000 1.000, n=0.000 1.000 0.000 +DEAL::p=1.936 1.000 1.000, n=0.000 1.000 0.000 +DEAL::p=2.404 1.000 1.000, n=0.000 1.000 0.000 +DEAL::p=2.872 1.000 1.000, n=0.000 1.000 0.000 +DEAL::p=1.000 1.000 1.468, n=0.000 1.000 0.000 +DEAL::p=1.468 1.000 1.468, n=0.000 1.000 0.000 +DEAL::p=1.936 1.000 1.468, n=0.000 1.000 0.000 +DEAL::p=2.404 1.000 1.468, n=0.000 1.000 0.000 +DEAL::p=2.872 1.000 1.468, n=0.000 1.000 0.000 +DEAL::p=1.000 1.000 1.936, n=0.000 1.000 0.000 +DEAL::p=1.468 1.000 1.936, n=0.000 1.000 0.000 +DEAL::p=1.936 1.000 1.936, n=0.000 1.000 0.000 +DEAL::p=2.404 1.000 1.936, n=0.000 1.000 0.000 +DEAL::p=2.872 1.000 1.936, n=0.000 1.000 0.000 +DEAL::p=1.000 1.000 2.404, n=0.000 1.000 0.000 +DEAL::p=1.468 1.000 2.404, n=0.000 1.000 0.000 +DEAL::p=1.936 1.000 2.404, n=0.000 1.000 0.000 +DEAL::p=2.404 1.000 2.404, n=0.000 1.000 0.000 +DEAL::p=2.872 1.000 2.404, n=0.000 1.000 0.000 +DEAL::p=1.000 1.000 2.872, n=0.000 1.000 0.000 +DEAL::p=1.468 1.000 2.872, n=0.000 1.000 0.000 +DEAL::p=1.936 1.000 2.872, n=0.000 1.000 0.000 +DEAL::p=2.404 1.000 2.872, n=0.000 1.000 0.000 +DEAL::p=2.872 1.000 2.872, n=0.000 1.000 0.000 +DEAL::p=1.000 3.000 1.000, n=0.000 1.000 0.000 +DEAL::p=1.468 3.000 1.000, n=0.000 1.000 0.000 +DEAL::p=1.936 3.000 1.000, n=0.000 1.000 0.000 +DEAL::p=2.404 3.000 1.000, n=0.000 1.000 0.000 +DEAL::p=2.872 3.000 1.000, n=0.000 1.000 0.000 +DEAL::p=1.000 3.000 1.468, n=0.000 1.000 0.000 +DEAL::p=1.468 3.000 1.468, n=0.000 1.000 0.000 +DEAL::p=1.936 3.000 1.468, n=0.000 1.000 0.000 +DEAL::p=2.404 3.000 1.468, n=0.000 1.000 0.000 +DEAL::p=2.872 3.000 1.468, n=0.000 1.000 0.000 +DEAL::p=1.000 3.000 1.936, n=0.000 1.000 0.000 +DEAL::p=1.468 3.000 1.936, n=0.000 1.000 0.000 +DEAL::p=1.936 3.000 1.936, n=0.000 1.000 0.000 +DEAL::p=2.404 3.000 1.936, n=0.000 1.000 0.000 +DEAL::p=2.872 3.000 1.936, n=0.000 1.000 0.000 +DEAL::p=1.000 3.000 2.404, n=0.000 1.000 0.000 +DEAL::p=1.468 3.000 2.404, n=0.000 1.000 0.000 +DEAL::p=1.936 3.000 2.404, n=0.000 1.000 0.000 +DEAL::p=2.404 3.000 2.404, n=0.000 1.000 0.000 +DEAL::p=2.872 3.000 2.404, n=0.000 1.000 0.000 +DEAL::p=1.000 3.000 2.872, n=0.000 1.000 0.000 +DEAL::p=1.468 3.000 2.872, n=0.000 1.000 0.000 +DEAL::p=1.936 3.000 2.872, n=0.000 1.000 0.000 +DEAL::p=2.404 3.000 2.872, n=0.000 1.000 0.000 +DEAL::p=2.872 3.000 2.872, n=0.000 1.000 0.000 +DEAL::p=1.000 1.000 1.000, n=0.000 0.000 1.000 +DEAL::p=1.000 1.468 1.000, n=0.000 0.000 1.000 +DEAL::p=1.000 1.936 1.000, n=0.000 0.000 1.000 +DEAL::p=1.000 2.404 1.000, n=0.000 0.000 1.000 +DEAL::p=1.000 2.872 1.000, n=0.000 0.000 1.000 +DEAL::p=1.468 1.000 1.000, n=0.000 0.000 1.000 +DEAL::p=1.468 1.468 1.000, n=0.000 0.000 1.000 +DEAL::p=1.468 1.936 1.000, n=0.000 0.000 1.000 +DEAL::p=1.468 2.404 1.000, n=0.000 0.000 1.000 +DEAL::p=1.468 2.872 1.000, n=0.000 0.000 1.000 +DEAL::p=1.936 1.000 1.000, n=0.000 0.000 1.000 +DEAL::p=1.936 1.468 1.000, n=0.000 0.000 1.000 +DEAL::p=1.936 1.936 1.000, n=0.000 0.000 1.000 +DEAL::p=1.936 2.404 1.000, n=0.000 0.000 1.000 +DEAL::p=1.936 2.872 1.000, n=0.000 0.000 1.000 +DEAL::p=2.404 1.000 1.000, n=0.000 0.000 1.000 +DEAL::p=2.404 1.468 1.000, n=0.000 0.000 1.000 +DEAL::p=2.404 1.936 1.000, n=0.000 0.000 1.000 +DEAL::p=2.404 2.404 1.000, n=0.000 0.000 1.000 +DEAL::p=2.404 2.872 1.000, n=0.000 0.000 1.000 +DEAL::p=2.872 1.000 1.000, n=0.000 0.000 1.000 +DEAL::p=2.872 1.468 1.000, n=0.000 0.000 1.000 +DEAL::p=2.872 1.936 1.000, n=0.000 0.000 1.000 +DEAL::p=2.872 2.404 1.000, n=0.000 0.000 1.000 +DEAL::p=2.872 2.872 1.000, n=0.000 0.000 1.000 +DEAL::p=1.000 1.000 3.000, n=0.000 0.000 1.000 +DEAL::p=1.000 1.468 3.000, n=0.000 0.000 1.000 +DEAL::p=1.000 1.936 3.000, n=0.000 0.000 1.000 +DEAL::p=1.000 2.404 3.000, n=0.000 0.000 1.000 +DEAL::p=1.000 2.872 3.000, n=0.000 0.000 1.000 +DEAL::p=1.468 1.000 3.000, n=0.000 0.000 1.000 +DEAL::p=1.468 1.468 3.000, n=0.000 0.000 1.000 +DEAL::p=1.468 1.936 3.000, n=0.000 0.000 1.000 +DEAL::p=1.468 2.404 3.000, n=0.000 0.000 1.000 +DEAL::p=1.468 2.872 3.000, n=0.000 0.000 1.000 +DEAL::p=1.936 1.000 3.000, n=0.000 0.000 1.000 +DEAL::p=1.936 1.468 3.000, n=0.000 0.000 1.000 +DEAL::p=1.936 1.936 3.000, n=0.000 0.000 1.000 +DEAL::p=1.936 2.404 3.000, n=0.000 0.000 1.000 +DEAL::p=1.936 2.872 3.000, n=0.000 0.000 1.000 +DEAL::p=2.404 1.000 3.000, n=0.000 0.000 1.000 +DEAL::p=2.404 1.468 3.000, n=0.000 0.000 1.000 +DEAL::p=2.404 1.936 3.000, n=0.000 0.000 1.000 +DEAL::p=2.404 2.404 3.000, n=0.000 0.000 1.000 +DEAL::p=2.404 2.872 3.000, n=0.000 0.000 1.000 +DEAL::p=2.872 1.000 3.000, n=0.000 0.000 1.000 +DEAL::p=2.872 1.468 3.000, n=0.000 0.000 1.000 +DEAL::p=2.872 1.936 3.000, n=0.000 0.000 1.000 +DEAL::p=2.872 2.404 3.000, n=0.000 0.000 1.000 +DEAL::p=2.872 2.872 3.000, n=0.000 0.000 1.000 +DEAL::Case1 +DEAL::p=0.000 -0.500 -1.000, n=0.968 -0.176 -0.176 +DEAL::p=0.234 -0.149 -0.064, n=0.971 -0.145 -0.189 +DEAL::p=0.468 0.202 0.872, n=0.973 -0.108 -0.203 +DEAL::p=0.702 0.553 1.808, n=0.974 -0.065 -0.219 +DEAL::p=0.936 0.904 2.744, n=0.971 -0.015 -0.237 +DEAL::p=0.234 0.319 -0.532, n=0.970 -0.193 -0.148 +DEAL::p=0.413 0.588 0.294, n=0.974 -0.159 -0.159 +DEAL::p=0.592 0.857 1.121, n=0.978 -0.120 -0.173 +DEAL::p=0.772 1.126 1.947, n=0.979 -0.073 -0.189 +DEAL::p=0.951 1.394 2.774, n=0.978 -0.017 -0.207 +DEAL::p=0.468 1.138 -0.064, n=0.971 -0.213 -0.113 +DEAL::p=0.592 1.325 0.653, n=0.976 -0.178 -0.123 +DEAL::p=0.717 1.511 1.370, n=0.982 -0.135 -0.135 +DEAL::p=0.841 1.698 2.087, n=0.985 -0.084 -0.149 +DEAL::p=0.966 1.885 2.804, n=0.986 -0.020 -0.166 +DEAL::p=0.702 1.957 0.404, n=0.969 -0.237 -0.071 +DEAL::p=0.772 2.062 1.011, n=0.977 -0.200 -0.078 +DEAL::p=0.841 2.166 1.619, n=0.984 -0.154 -0.086 +DEAL::p=0.911 2.271 2.226, n=0.991 -0.097 -0.097 +DEAL::p=0.981 2.375 2.834, n=0.994 -0.024 -0.110 +DEAL::p=0.936 2.776 0.872, n=0.964 -0.266 -0.017 +DEAL::p=0.951 2.798 1.370, n=0.974 -0.228 -0.019 +DEAL::p=0.966 2.821 1.868, n=0.984 -0.179 -0.022 +DEAL::p=0.981 2.843 2.366, n=0.993 -0.115 -0.025 +DEAL::p=0.996 2.866 2.864, n=0.999 -0.029 -0.029 +DEAL::p=1.250 0.250 0.250, n=0.816 -0.408 -0.408 +DEAL::p=1.659 0.425 0.894, n=0.833 -0.336 -0.439 +DEAL::p=2.069 0.601 1.537, n=0.846 -0.250 -0.470 +DEAL::p=2.479 0.777 2.181, n=0.852 -0.149 -0.502 +DEAL::p=2.888 0.952 2.824, n=0.847 -0.034 -0.530 +DEAL::p=1.659 0.894 0.425, n=0.833 -0.439 -0.336 +DEAL::p=1.973 1.028 1.028, n=0.857 -0.365 -0.365 +DEAL::p=2.287 1.162 1.630, n=0.877 -0.274 -0.395 +DEAL::p=2.601 1.297 2.233, n=0.889 -0.166 -0.426 +DEAL::p=2.914 1.431 2.835, n=0.890 -0.038 -0.455 +DEAL::p=2.069 1.537 0.601, n=0.846 -0.470 -0.250 +DEAL::p=2.287 1.630 1.162, n=0.877 -0.395 -0.274 +DEAL::p=2.505 1.724 1.724, n=0.905 -0.301 -0.301 +DEAL::p=2.723 1.817 2.285, n=0.926 -0.184 -0.329 +DEAL::p=2.940 1.910 2.846, n=0.934 -0.043 -0.355 +DEAL::p=2.479 2.181 0.777, n=0.852 -0.502 -0.149 +DEAL::p=2.601 2.233 1.297, n=0.889 -0.426 -0.166 +DEAL::p=2.723 2.285 1.817, n=0.926 -0.329 -0.184 +DEAL::p=2.845 2.337 2.337, n=0.957 -0.204 -0.204 +DEAL::p=2.967 2.390 2.858, n=0.974 -0.048 -0.224 +DEAL::p=2.888 2.824 0.952, n=0.847 -0.530 -0.034 +DEAL::p=2.914 2.835 1.431, n=0.890 -0.455 -0.038 +DEAL::p=2.940 2.846 1.910, n=0.934 -0.355 -0.043 +DEAL::p=2.967 2.858 2.390, n=0.974 -0.224 -0.048 +DEAL::p=2.993 2.869 2.869, n=0.997 -0.053 -0.053 +DEAL::p=0.000 -0.500 -1.000, n=-0.276 0.921 -0.276 +DEAL::p=0.293 -0.325 -0.708, n=-0.326 0.918 -0.225 +DEAL::p=0.585 -0.149 -0.415, n=-0.395 0.906 -0.149 +DEAL::p=0.878 0.027 -0.122, n=-0.491 0.871 -0.031 +DEAL::p=1.170 0.202 0.170, n=-0.616 0.772 0.153 +DEAL::p=0.234 -0.149 -0.064, n=-0.170 0.936 -0.308 +DEAL::p=0.568 -0.015 0.160, n=-0.195 0.942 -0.275 +DEAL::p=0.901 0.120 0.384, n=-0.228 0.947 -0.228 +DEAL::p=1.235 0.254 0.608, n=-0.273 0.948 -0.162 +DEAL::p=1.568 0.389 0.832, n=-0.337 0.940 -0.062 +DEAL::p=0.468 0.202 0.872, n=-0.098 0.940 -0.328 +DEAL::p=0.843 0.295 1.028, n=-0.110 0.947 -0.303 +DEAL::p=1.217 0.389 1.183, n=-0.125 0.954 -0.272 +DEAL::p=1.592 0.482 1.339, n=-0.145 0.963 -0.229 +DEAL::p=1.967 0.575 1.494, n=-0.171 0.970 -0.171 +DEAL::p=0.702 0.553 1.808, n=-0.047 0.939 -0.340 +DEAL::p=1.118 0.605 1.895, n=-0.052 0.946 -0.321 +DEAL::p=1.533 0.658 1.982, n=-0.057 0.953 -0.298 +DEAL::p=1.949 0.710 2.069, n=-0.065 0.961 -0.268 +DEAL::p=2.365 0.762 2.157, n=-0.074 0.971 -0.229 +DEAL::p=0.936 0.904 2.744, n=0.009 0.937 -0.349 +DEAL::p=1.393 0.915 2.763, n=-0.010 0.943 -0.334 +DEAL::p=1.850 0.926 2.781, n=-0.010 0.949 -0.315 +DEAL::p=2.306 0.938 2.800, n=-0.012 0.956 -0.292 +DEAL::p=2.763 0.949 2.819, n=-0.013 0.964 -0.264 +DEAL::p=1.000 3.000 1.000, n=0.000 1.000 0.000 +DEAL::p=1.468 3.000 1.000, n=0.000 1.000 0.000 +DEAL::p=1.936 3.000 1.000, n=0.000 1.000 0.000 +DEAL::p=2.404 3.000 1.000, n=0.000 1.000 0.000 +DEAL::p=2.872 3.000 1.000, n=0.000 1.000 0.000 +DEAL::p=1.000 3.000 1.468, n=0.000 1.000 0.000 +DEAL::p=1.468 3.000 1.468, n=0.000 1.000 0.000 +DEAL::p=1.936 3.000 1.468, n=0.000 1.000 0.000 +DEAL::p=2.404 3.000 1.468, n=0.000 1.000 0.000 +DEAL::p=2.872 3.000 1.468, n=0.000 1.000 0.000 +DEAL::p=1.000 3.000 1.936, n=0.000 1.000 0.000 +DEAL::p=1.468 3.000 1.936, n=0.000 1.000 0.000 +DEAL::p=1.936 3.000 1.936, n=0.000 1.000 0.000 +DEAL::p=2.404 3.000 1.936, n=0.000 1.000 0.000 +DEAL::p=2.872 3.000 1.936, n=0.000 1.000 0.000 +DEAL::p=1.000 3.000 2.404, n=0.000 1.000 0.000 +DEAL::p=1.468 3.000 2.404, n=0.000 1.000 0.000 +DEAL::p=1.936 3.000 2.404, n=0.000 1.000 0.000 +DEAL::p=2.404 3.000 2.404, n=0.000 1.000 0.000 +DEAL::p=2.872 3.000 2.404, n=0.000 1.000 0.000 +DEAL::p=1.000 3.000 2.872, n=0.000 1.000 0.000 +DEAL::p=1.468 3.000 2.872, n=0.000 1.000 0.000 +DEAL::p=1.936 3.000 2.872, n=0.000 1.000 0.000 +DEAL::p=2.404 3.000 2.872, n=0.000 1.000 0.000 +DEAL::p=2.872 3.000 2.872, n=0.000 1.000 0.000 +DEAL::p=0.000 -0.500 -1.000, n=-0.600 -0.261 0.756 +DEAL::p=0.234 0.319 -0.532, n=-0.417 -0.358 0.835 +DEAL::p=0.468 1.138 -0.064, n=-0.255 -0.424 0.869 +DEAL::p=0.702 1.957 0.404, n=-0.125 -0.465 0.876 +DEAL::p=0.936 2.776 0.872, n=-0.024 -0.491 0.871 +DEAL::p=0.293 -0.325 -0.708, n=-0.652 -0.151 0.743 +DEAL::p=0.568 0.453 -0.308, n=-0.458 -0.272 0.846 +DEAL::p=0.843 1.231 0.092, n=-0.281 -0.358 0.890 +DEAL::p=1.118 2.009 0.491, n=-0.137 -0.414 0.900 +DEAL::p=1.393 2.787 0.891, n=-0.026 -0.450 0.893 +DEAL::p=0.585 -0.149 -0.415, n=-0.701 -0.020 0.713 +DEAL::p=0.901 0.588 -0.084, n=-0.503 -0.165 0.848 +DEAL::p=1.217 1.325 0.247, n=-0.309 -0.276 0.910 +DEAL::p=1.533 2.062 0.578, n=-0.150 -0.351 0.924 +DEAL::p=1.850 2.798 0.909, n=-0.028 -0.400 0.916 +DEAL::p=0.878 0.027 -0.122, n=-0.739 0.130 0.661 +DEAL::p=1.235 0.722 0.140, n=-0.548 -0.034 0.836 +DEAL::p=1.592 1.418 0.403, n=-0.340 -0.174 0.924 +DEAL::p=1.949 2.114 0.665, n=-0.164 -0.273 0.948 +DEAL::p=2.306 2.810 0.928, n=-0.030 -0.339 0.940 +DEAL::p=1.170 0.202 0.170, n=-0.758 0.287 0.586 +DEAL::p=1.568 0.857 0.364, n=-0.586 0.119 0.801 +DEAL::p=1.967 1.511 0.558, n=-0.373 -0.048 0.927 +DEAL::p=2.365 2.166 0.753, n=-0.181 -0.177 0.967 +DEAL::p=2.763 2.821 0.947, n=-0.033 -0.266 0.963 +DEAL::p=1.000 1.000 3.000, n=0.000 0.000 1.000 +DEAL::p=1.000 1.468 3.000, n=0.000 0.000 1.000 +DEAL::p=1.000 1.936 3.000, n=0.000 0.000 1.000 +DEAL::p=1.000 2.404 3.000, n=0.000 0.000 1.000 +DEAL::p=1.000 2.872 3.000, n=0.000 0.000 1.000 +DEAL::p=1.468 1.000 3.000, n=0.000 0.000 1.000 +DEAL::p=1.468 1.468 3.000, n=0.000 0.000 1.000 +DEAL::p=1.468 1.936 3.000, n=0.000 0.000 1.000 +DEAL::p=1.468 2.404 3.000, n=0.000 0.000 1.000 +DEAL::p=1.468 2.872 3.000, n=0.000 0.000 1.000 +DEAL::p=1.936 1.000 3.000, n=0.000 0.000 1.000 +DEAL::p=1.936 1.468 3.000, n=0.000 0.000 1.000 +DEAL::p=1.936 1.936 3.000, n=0.000 0.000 1.000 +DEAL::p=1.936 2.404 3.000, n=0.000 0.000 1.000 +DEAL::p=1.936 2.872 3.000, n=0.000 0.000 1.000 +DEAL::p=2.404 1.000 3.000, n=0.000 0.000 1.000 +DEAL::p=2.404 1.468 3.000, n=0.000 0.000 1.000 +DEAL::p=2.404 1.936 3.000, n=0.000 0.000 1.000 +DEAL::p=2.404 2.404 3.000, n=0.000 0.000 1.000 +DEAL::p=2.404 2.872 3.000, n=0.000 0.000 1.000 +DEAL::p=2.872 1.000 3.000, n=0.000 0.000 1.000 +DEAL::p=2.872 1.468 3.000, n=0.000 0.000 1.000 +DEAL::p=2.872 1.936 3.000, n=0.000 0.000 1.000 +DEAL::p=2.872 2.404 3.000, n=0.000 0.000 1.000 +DEAL::p=2.872 2.872 3.000, n=0.000 0.000 1.000 +DEAL::OK diff --git a/tests/deal.II/normal_vector_03_2d.cc b/tests/deal.II/normal_vector_03_2d.cc new file mode 100644 index 0000000000..e79dd1ce73 --- /dev/null +++ b/tests/deal.II/normal_vector_03_2d.cc @@ -0,0 +1,114 @@ +//---------------------------- normal_vector_03_2d.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2005, 2008, 2011 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- normal_vector_03_2d.cc --------------------------- + + +// like _01 but instead for random points on the face. note that the normal +// vector on a bilinearly mapped face is also a bilinear function + + + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + + + +void create_triangulation(const unsigned int case_no, + Triangulation<2> &tria) +{ + switch (case_no) + { + case 0: + GridGenerator::hyper_cube(tria, 1., 3.); + break; + case 1: + { + GridGenerator::hyper_cube(tria, 1., 3.); + Point<2> &v0=tria.begin_active()->vertex(0); + v0 = Point<2> (-0.5,-1); + Point<2> &v1=tria.begin_active()->vertex(1); + v1 = Point<2> (1.25, 0.25); + break; + } + default: + Assert(false, ExcNotImplemented()); + }; +} + + + +int main () +{ + std::ofstream logfile ("normal_vector_03_2d/output"); + deallog << std::setprecision (3); + deallog << std::fixed; + deallog.attach(logfile); + deallog.depth_console (0); + + FE_Q<2> linear_interpolator(1); + + Triangulation<2> tria; + StraightBoundary<2> boundary; + Boundary<2>::FaceVertexNormals normals; + for (unsigned int case_no=0; case_no<2; ++case_no) + { + deallog << "Case" << case_no << std::endl; + create_triangulation(case_no, tria); + const Triangulation<2>::active_cell_iterator cell=tria.begin_active(); + Triangulation<2>::face_iterator face; + for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) + { + face=cell->face(face_no); + boundary.get_normals_at_vertices(face, normals); + + for (double xi=0; xi<=1; xi+=0.234) + for (double eta=0; eta<=1; eta+=0.234) + { + Point<2> p; + Tensor<1,2> normal; + + for (unsigned int v=0; v::vertices_per_face; ++v) + { + p += face->vertex(v) * linear_interpolator.shape_value(v,Point<2>(xi,eta)); + normal += normals[v] * + linear_interpolator.shape_value(v,Point<2>(xi,eta)); + } + normal /= normal.norm(); + + deallog << "p=" << p + << ", n=" << boundary.normal_vector (face, p) + << std::endl; + + Assert ((boundary.normal_vector (face, p) + - + normal).norm() + < + 1e-10, + ExcInternalError()); + } + } + + tria.clear(); + } + + deallog << "OK" << std::endl; +} diff --git a/tests/deal.II/normal_vector_03_2d/cmp/generic b/tests/deal.II/normal_vector_03_2d/cmp/generic new file mode 100644 index 0000000000..209f98de7b --- /dev/null +++ b/tests/deal.II/normal_vector_03_2d/cmp/generic @@ -0,0 +1,204 @@ + +DEAL::Case0 +DEAL::p=1.000 1.000, n=1.000 0.000 +DEAL::p=0.766 0.766, n=1.000 0.000 +DEAL::p=0.532 0.532, n=1.000 0.000 +DEAL::p=0.298 0.298, n=1.000 0.000 +DEAL::p=0.064 0.064, n=1.000 0.000 +DEAL::p=1.000 1.468, n=1.000 0.000 +DEAL::p=0.766 1.124, n=1.000 0.000 +DEAL::p=0.532 0.781, n=1.000 0.000 +DEAL::p=0.298 0.437, n=1.000 0.000 +DEAL::p=0.064 0.094, n=1.000 0.000 +DEAL::p=1.000 1.936, n=1.000 0.000 +DEAL::p=0.766 1.483, n=1.000 0.000 +DEAL::p=0.532 1.030, n=1.000 0.000 +DEAL::p=0.298 0.577, n=1.000 0.000 +DEAL::p=0.064 0.124, n=1.000 0.000 +DEAL::p=1.000 2.404, n=1.000 0.000 +DEAL::p=0.766 1.841, n=1.000 0.000 +DEAL::p=0.532 1.279, n=1.000 0.000 +DEAL::p=0.298 0.716, n=1.000 0.000 +DEAL::p=0.064 0.154, n=1.000 0.000 +DEAL::p=1.000 2.872, n=1.000 0.000 +DEAL::p=0.766 2.200, n=1.000 0.000 +DEAL::p=0.532 1.528, n=1.000 0.000 +DEAL::p=0.298 0.856, n=1.000 0.000 +DEAL::p=0.064 0.184, n=1.000 0.000 +DEAL::p=3.000 1.000, n=1.000 0.000 +DEAL::p=2.298 0.766, n=1.000 0.000 +DEAL::p=1.596 0.532, n=1.000 0.000 +DEAL::p=0.894 0.298, n=1.000 0.000 +DEAL::p=0.192 0.064, n=1.000 0.000 +DEAL::p=3.000 1.468, n=1.000 0.000 +DEAL::p=2.298 1.124, n=1.000 0.000 +DEAL::p=1.596 0.781, n=1.000 0.000 +DEAL::p=0.894 0.437, n=1.000 0.000 +DEAL::p=0.192 0.094, n=1.000 0.000 +DEAL::p=3.000 1.936, n=1.000 0.000 +DEAL::p=2.298 1.483, n=1.000 0.000 +DEAL::p=1.596 1.030, n=1.000 0.000 +DEAL::p=0.894 0.577, n=1.000 0.000 +DEAL::p=0.192 0.124, n=1.000 0.000 +DEAL::p=3.000 2.404, n=1.000 0.000 +DEAL::p=2.298 1.841, n=1.000 0.000 +DEAL::p=1.596 1.279, n=1.000 0.000 +DEAL::p=0.894 0.716, n=1.000 0.000 +DEAL::p=0.192 0.154, n=1.000 0.000 +DEAL::p=3.000 2.872, n=1.000 0.000 +DEAL::p=2.298 2.200, n=1.000 0.000 +DEAL::p=1.596 1.528, n=1.000 0.000 +DEAL::p=0.894 0.856, n=1.000 0.000 +DEAL::p=0.192 0.184, n=1.000 0.000 +DEAL::p=1.000 1.000, n=0.000 -1.000 +DEAL::p=0.766 0.766, n=0.000 -1.000 +DEAL::p=0.532 0.532, n=0.000 -1.000 +DEAL::p=0.298 0.298, n=0.000 -1.000 +DEAL::p=0.064 0.064, n=0.000 -1.000 +DEAL::p=1.468 1.000, n=0.000 -1.000 +DEAL::p=1.124 0.766, n=0.000 -1.000 +DEAL::p=0.781 0.532, n=0.000 -1.000 +DEAL::p=0.437 0.298, n=0.000 -1.000 +DEAL::p=0.094 0.064, n=0.000 -1.000 +DEAL::p=1.936 1.000, n=0.000 -1.000 +DEAL::p=1.483 0.766, n=0.000 -1.000 +DEAL::p=1.030 0.532, n=0.000 -1.000 +DEAL::p=0.577 0.298, n=0.000 -1.000 +DEAL::p=0.124 0.064, n=0.000 -1.000 +DEAL::p=2.404 1.000, n=0.000 -1.000 +DEAL::p=1.841 0.766, n=0.000 -1.000 +DEAL::p=1.279 0.532, n=0.000 -1.000 +DEAL::p=0.716 0.298, n=0.000 -1.000 +DEAL::p=0.154 0.064, n=0.000 -1.000 +DEAL::p=2.872 1.000, n=0.000 -1.000 +DEAL::p=2.200 0.766, n=0.000 -1.000 +DEAL::p=1.528 0.532, n=0.000 -1.000 +DEAL::p=0.856 0.298, n=0.000 -1.000 +DEAL::p=0.184 0.064, n=0.000 -1.000 +DEAL::p=1.000 3.000, n=0.000 -1.000 +DEAL::p=0.766 2.298, n=0.000 -1.000 +DEAL::p=0.532 1.596, n=0.000 -1.000 +DEAL::p=0.298 0.894, n=0.000 -1.000 +DEAL::p=0.064 0.192, n=0.000 -1.000 +DEAL::p=1.468 3.000, n=0.000 -1.000 +DEAL::p=1.124 2.298, n=0.000 -1.000 +DEAL::p=0.781 1.596, n=0.000 -1.000 +DEAL::p=0.437 0.894, n=0.000 -1.000 +DEAL::p=0.094 0.192, n=0.000 -1.000 +DEAL::p=1.936 3.000, n=0.000 -1.000 +DEAL::p=1.483 2.298, n=0.000 -1.000 +DEAL::p=1.030 1.596, n=0.000 -1.000 +DEAL::p=0.577 0.894, n=0.000 -1.000 +DEAL::p=0.124 0.192, n=0.000 -1.000 +DEAL::p=2.404 3.000, n=0.000 -1.000 +DEAL::p=1.841 2.298, n=0.000 -1.000 +DEAL::p=1.279 1.596, n=0.000 -1.000 +DEAL::p=0.716 0.894, n=0.000 -1.000 +DEAL::p=0.154 0.192, n=0.000 -1.000 +DEAL::p=2.872 3.000, n=0.000 -1.000 +DEAL::p=2.200 2.298, n=0.000 -1.000 +DEAL::p=1.528 1.596, n=0.000 -1.000 +DEAL::p=0.856 0.894, n=0.000 -1.000 +DEAL::p=0.184 0.192, n=0.000 -1.000 +DEAL::Case1 +DEAL::p=-0.500 -1.000, n=0.936 -0.351 +DEAL::p=-0.383 -0.766, n=0.936 -0.351 +DEAL::p=-0.266 -0.532, n=0.936 -0.351 +DEAL::p=-0.149 -0.298, n=0.936 -0.351 +DEAL::p=-0.032 -0.064, n=0.936 -0.351 +DEAL::p=-0.149 -0.064, n=0.936 -0.351 +DEAL::p=-0.114 -0.049, n=0.936 -0.351 +DEAL::p=-0.079 -0.034, n=0.936 -0.351 +DEAL::p=-0.044 -0.019, n=0.936 -0.351 +DEAL::p=-0.010 0.004, n=0.936 -0.351 +DEAL::p=0.202 0.872, n=0.936 -0.351 +DEAL::p=0.155 0.668, n=0.936 -0.351 +DEAL::p=0.107 0.464, n=0.936 -0.351 +DEAL::p=0.060 0.260, n=0.936 -0.351 +DEAL::p=0.013 0.056, n=0.936 -0.351 +DEAL::p=0.553 1.808, n=0.936 -0.351 +DEAL::p=0.424 1.385, n=0.936 -0.351 +DEAL::p=0.294 0.962, n=0.936 -0.351 +DEAL::p=0.165 0.539, n=0.936 -0.351 +DEAL::p=0.035 0.116, n=0.936 -0.351 +DEAL::p=0.904 2.744, n=0.936 -0.351 +DEAL::p=0.692 2.102, n=0.936 -0.351 +DEAL::p=0.481 1.460, n=0.936 -0.351 +DEAL::p=0.269 0.818, n=0.936 -0.351 +DEAL::p=0.058 0.176, n=0.936 -0.351 +DEAL::p=1.250 0.250, n=0.844 -0.537 +DEAL::p=0.958 0.192, n=0.844 -0.537 +DEAL::p=0.665 0.133, n=0.844 -0.537 +DEAL::p=0.372 0.074, n=0.844 -0.537 +DEAL::p=0.080 0.016, n=0.844 -0.537 +DEAL::p=1.659 0.894, n=0.844 -0.537 +DEAL::p=1.271 0.684, n=0.844 -0.537 +DEAL::p=0.883 0.475, n=0.844 -0.537 +DEAL::p=0.495 0.266, n=0.844 -0.537 +DEAL::p=0.106 0.057, n=0.844 -0.537 +DEAL::p=2.069 1.537, n=0.844 -0.537 +DEAL::p=1.585 1.177, n=0.844 -0.537 +DEAL::p=1.101 0.818, n=0.844 -0.537 +DEAL::p=0.617 0.458, n=0.844 -0.537 +DEAL::p=0.132 0.098, n=0.844 -0.537 +DEAL::p=2.479 2.181, n=0.844 -0.537 +DEAL::p=1.899 1.670, n=0.844 -0.537 +DEAL::p=1.319 1.160, n=0.844 -0.537 +DEAL::p=0.739 0.650, n=0.844 -0.537 +DEAL::p=0.159 0.140, n=0.844 -0.537 +DEAL::p=2.888 2.824, n=0.844 -0.537 +DEAL::p=2.212 2.163, n=0.844 -0.537 +DEAL::p=1.536 1.502, n=0.844 -0.537 +DEAL::p=0.861 0.842, n=0.844 -0.537 +DEAL::p=0.185 0.181, n=0.844 -0.537 +DEAL::p=-0.500 -1.000, n=0.581 -0.814 +DEAL::p=-0.383 -0.766, n=0.581 -0.814 +DEAL::p=-0.266 -0.532, n=0.581 -0.814 +DEAL::p=-0.149 -0.298, n=0.581 -0.814 +DEAL::p=-0.032 -0.064, n=0.581 -0.814 +DEAL::p=-0.090 -0.708, n=0.581 -0.814 +DEAL::p=-0.069 -0.542, n=0.581 -0.814 +DEAL::p=-0.048 -0.376, n=0.581 -0.814 +DEAL::p=-0.027 -0.211, n=0.581 -0.814 +DEAL::p=0.006 -0.045, n=0.581 -0.814 +DEAL::p=0.319 -0.415, n=0.581 -0.814 +DEAL::p=0.244 -0.318, n=0.581 -0.814 +DEAL::p=0.170 -0.221, n=0.581 -0.814 +DEAL::p=0.095 -0.124, n=0.581 -0.814 +DEAL::p=0.020 -0.027, n=0.581 -0.814 +DEAL::p=0.729 -0.122, n=0.581 -0.814 +DEAL::p=0.558 -0.094, n=0.581 -0.814 +DEAL::p=0.388 -0.065, n=0.581 -0.814 +DEAL::p=0.217 -0.037, n=0.581 -0.814 +DEAL::p=0.047 0.008, n=0.581 -0.814 +DEAL::p=1.138 0.170, n=0.581 -0.814 +DEAL::p=0.872 0.130, n=0.581 -0.814 +DEAL::p=0.605 0.090, n=0.581 -0.814 +DEAL::p=0.339 0.051, n=0.581 -0.814 +DEAL::p=0.073 0.011, n=0.581 -0.814 +DEAL::p=1.000 3.000, n=0.000 -1.000 +DEAL::p=0.766 2.298, n=0.000 -1.000 +DEAL::p=0.532 1.596, n=0.000 -1.000 +DEAL::p=0.298 0.894, n=0.000 -1.000 +DEAL::p=0.064 0.192, n=0.000 -1.000 +DEAL::p=1.468 3.000, n=0.000 -1.000 +DEAL::p=1.124 2.298, n=0.000 -1.000 +DEAL::p=0.781 1.596, n=0.000 -1.000 +DEAL::p=0.437 0.894, n=0.000 -1.000 +DEAL::p=0.094 0.192, n=0.000 -1.000 +DEAL::p=1.936 3.000, n=0.000 -1.000 +DEAL::p=1.483 2.298, n=0.000 -1.000 +DEAL::p=1.030 1.596, n=0.000 -1.000 +DEAL::p=0.577 0.894, n=0.000 -1.000 +DEAL::p=0.124 0.192, n=0.000 -1.000 +DEAL::p=2.404 3.000, n=0.000 -1.000 +DEAL::p=1.841 2.298, n=0.000 -1.000 +DEAL::p=1.279 1.596, n=0.000 -1.000 +DEAL::p=0.716 0.894, n=0.000 -1.000 +DEAL::p=0.154 0.192, n=0.000 -1.000 +DEAL::p=2.872 3.000, n=0.000 -1.000 +DEAL::p=2.200 2.298, n=0.000 -1.000 +DEAL::p=1.528 1.596, n=0.000 -1.000 +DEAL::p=0.856 0.894, n=0.000 -1.000 +DEAL::p=0.184 0.192, n=0.000 -1.000 +DEAL::OK -- 2.39.5