* lines therefore is also on the
* boundary).
*/
- virtual Point<spacedim>
+ virtual
+ Point<spacedim>
get_new_point_on_line (const typename Triangulation<dim,spacedim>::line_iterator &line) const = 0;
/**
* default implementation throws
* an error in any case, however.
*/
- virtual Point<spacedim>
+ virtual
+ Point<spacedim>
get_new_point_on_quad (const typename Triangulation<dim,spacedim>::quad_iterator &quad) const;
/**
* default implementation throws
* an error in any case, however.
*/
- virtual void
+ virtual
+ void
get_intermediate_points_on_line (const typename Triangulation<dim,spacedim>::line_iterator &line,
std::vector<Point<spacedim> > &points) const;
* default implementation throws
* an error in any case, however.
*/
- virtual void
+ virtual
+ void
get_intermediate_points_on_quad (const typename Triangulation<dim,spacedim>::quad_iterator &quad,
std::vector<Point<spacedim> > &points) const;
get_intermediate_points_on_face (const typename Triangulation<dim,spacedim>::face_iterator &face,
std::vector<Point<spacedim> > &points) const;
+ /**
+ * Return the normal vector to the surface
+ * at the point p. If p is not in fact
+ * on the surface, but only 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<dim,spacedim>::face_iterator &face,
+ const Point<spacedim> &p) const;
+
/**
* Compute the normal vectors to
* the boundary at each vertex of
* with respect to points inside
* the given face.
*/
- virtual void
+ virtual
+ void
get_normals_at_vertices (const typename Triangulation<dim,spacedim>::face_iterator &face,
FaceVertexNormals &face_vertex_normals) const;
* base class for more
* information.
*/
- virtual Point<spacedim>
+ virtual
+ Point<spacedim>
get_new_point_on_quad (const typename Triangulation<dim,spacedim>::quad_iterator &quad) const;
/**
* and the documentation of the
* base class.
*/
- virtual void
+ virtual
+ void
get_intermediate_points_on_line (const typename Triangulation<dim,spacedim>::line_iterator &line,
std::vector<Point<spacedim> > &points) const;
* and the documentation of the
* base class.
*/
- virtual void
+ virtual
+ void
get_intermediate_points_on_quad (const typename Triangulation<dim,spacedim>::quad_iterator &quad,
std::vector<Point<spacedim> > &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<dim,spacedim>::face_iterator &face,
+ const Point<spacedim> &p) const;
+
+ /**
* Compute the normals to the
* boundary at the vertices of
* the given face.
* and the documentation of the
* base class.
*/
- virtual void
+ virtual
+ void
get_normals_at_vertices (const typename Triangulation<dim,spacedim>::face_iterator &face,
typename Boundary<dim,spacedim>::FaceVertexNormals &face_vertex_normals) const;
* this class and the documentation of the
* base class.
*/
- virtual Point<spacedim>
+ virtual
+ Point<spacedim>
get_new_point_on_line (const typename Triangulation<dim,spacedim>::line_iterator &line) const;
/**
* this class and the documentation of the
* base class.
*/
- virtual Point<spacedim>
+ virtual
+ Point<spacedim>
get_new_point_on_quad (const typename Triangulation<dim,spacedim>::quad_iterator &quad) const;
/**
* Calls
* @p get_intermediate_points_between_points.
*/
- virtual void
+ virtual
+ void
get_intermediate_points_on_line (const typename Triangulation<dim,spacedim>::line_iterator &line,
std::vector<Point<spacedim> > &points) const;
* Only implemented for <tt>dim=3</tt>
* and for <tt>points.size()==1</tt>.
*/
- virtual void
+ virtual
+ void
get_intermediate_points_on_quad (const typename Triangulation<dim,spacedim>::quad_iterator &quad,
std::vector<Point<spacedim> > &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<dim,spacedim>::face_iterator &face,
+ const Point<spacedim> &p) const;
+
+ /**
* Compute the normals to the
* boundary at the vertices of
* the given face.
* and the documentation of the
* base class.
*/
- virtual void
+ virtual
+ void
get_normals_at_vertices (const typename Triangulation<dim,spacedim>::face_iterator &face,
typename Boundary<dim,spacedim>::FaceVertexNormals &face_vertex_normals) const;
/**
* Return the center of the ball.
*/
- Point<spacedim> get_center () const;
+ Point<spacedim>
+ get_center () const;
/**
* Return the radius of the ball.
*/
- double get_radius () const;
+ double
+ get_radius () const;
/**
* Exception. Thrown by the
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/fe/fe_q.h>
#include <cmath>
DEAL_II_NAMESPACE_OPEN
+template <int dim, int spacedim>
+Tensor<1,spacedim>
+Boundary<dim, spacedim>::
+normal_vector (const typename Triangulation<dim, spacedim>::face_iterator &,
+ const Point<spacedim> &) const
+{
+ Assert (false, ExcPureFunctionCalled());
+ return Tensor<1,spacedim>();
+}
+
+
+
template <int dim, int spacedim>
void
Boundary<dim, spacedim>::
+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 <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;
+
+ FE_Q<facedim> linear_fe(1);
+
+ const double eps = 1e-12;
+ Tensor<1,spacedim> grad_F[facedim];
+ while (true)
+ {
+ Point<spacedim> F;
+ for (unsigned int v=0; v<GeometryInfo<facedim>::vertices_per_cell; ++v)
+ F += face->vertex(v) * linear_fe.shape_value(v, xi);
+
+ 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) * linear_fe.shape_grad(v, xi)[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 Point<facedim> 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>::
return Point<1>();
}
+
template <>
Point<2>
HyperBallBoundary<1,2>::
}
+
template <int dim, int spacedim>
void
HyperBallBoundary<dim,spacedim>::get_intermediate_points_between_points (
+template <int dim, int spacedim>
+Tensor<1,spacedim>
+HyperBallBoundary<dim,spacedim>::
+normal_vector (const typename Triangulation<dim,spacedim>::face_iterator &,
+ const Point<spacedim> &p) const
+{
+ const Tensor<1,spacedim> unnormalized_normal = p-center;
+ return unnormalized_normal/unnormalized_normal.norm();
+}
+
+
+
template <>
void
HyperBallBoundary<1>::
--- /dev/null
+//---------------------------- 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 <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_boundary.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+
+#include <fstream>
+#include <iomanip>
+
+
+
+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<GeometryInfo<3>::faces_per_cell; ++face_no)
+ {
+ face=cell->face(face_no);
+ boundary.get_normals_at_vertices(face, normals);
+ for (unsigned int v=0; v<GeometryInfo<3>::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;
+}
--- /dev/null
+
+DEAL::Case0
+DEAL::Case1
+DEAL::OK
--- /dev/null
+//---------------------------- 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 <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_boundary.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+
+#include <fstream>
+#include <iomanip>
+
+
+
+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<GeometryInfo<2>::faces_per_cell; ++face_no)
+ {
+ face=cell->face(face_no);
+ boundary.get_normals_at_vertices(face, normals);
+ for (unsigned int v=0; v<GeometryInfo<2>::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;
+}
--- /dev/null
+
+DEAL::Case0
+DEAL::Case1
+DEAL::OK
--- /dev/null
+//---------------------------- 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 <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_boundary.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+
+#include <fstream>
+#include <iomanip>
+
+
+
+
+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<GeometryInfo<3>::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<GeometryInfo<3>::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;
+}
--- /dev/null
+
+DEAL::OK
--- /dev/null
+//---------------------------- 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 <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_boundary.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+
+#include <fstream>
+#include <iomanip>
+
+
+
+
+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<GeometryInfo<2>::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<GeometryInfo<2>::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;
+}
--- /dev/null
+
+DEAL::OK
--- /dev/null
+//---------------------------- 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 <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_boundary.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/fe/fe_q.h>
+
+#include <fstream>
+#include <iomanip>
+
+
+
+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<GeometryInfo<3>::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<GeometryInfo<3>::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;
+}
--- /dev/null
+
+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
--- /dev/null
+//---------------------------- 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 <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_boundary.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/fe/fe_q.h>
+
+#include <fstream>
+#include <iomanip>
+
+
+
+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<GeometryInfo<2>::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<GeometryInfo<2>::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;
+}
--- /dev/null
+
+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