* rotates it outward (with respect to the coordinate system of the edge)
* by 90 degrees. In 3d, the default implementation is more
* complicated, aiming at avoiding problems with numerical round-off
- * for points close to one of the vertices. If the point p is closer
- * to the center of the face than to any of the vertices, the
- * normal vector is computed by the cross product of the tangent
- * vectors from p to either vertex zero and one of the face (if
- * the closest vertex is either vertex two or three), or of the tangent
- * vectors from p to vertices two and three (if the closest vertex is
- * either vertex zero or one). On the other hand, if the point p
- * is closer to one of the vertices than to the center of the face,
- * then we take the cross product of the tangent vectors from p
- * to the two vertices that are adjacent to the closest one.
+ * for points close to one of the vertices, and avoiding tangent directions
+ * that are linearly dependent.
*/
virtual
Tensor<1,spacedim>
{
const int spacedim=3;
Tensor<1,spacedim> t1,t2;
+ Tensor<1,spacedim> normal;
- // Take the difference between p and all four vertices
- unsigned int min_index=0;
- double min_distance = (p-face->vertex(0)).norm_square();
-
- for (unsigned int i=1; i<4; ++i)
- {
- const Tensor<1,spacedim> dp = p-face->vertex(i);
- double distance = dp.norm_square();
- if (distance < min_distance)
- {
- min_index = i;
- min_distance = distance;
- }
- }
- // Verify we have a valid vertex index
- AssertIndexRange(min_index, 4);
-
- // Now figure out which vertices are best to compute tangent vectors.
- // We split the cell in a central diamond of points closer to the
- // center than to any of the vertices, and the 4 triangles in the
- // corner. The central diamond is split into its upper and lower
- // half. For each of these 6 cases, the following encodes a list
- // of two vertices each to which we compute the tangent vectors,
- // and then take the cross product. See the documentation of this
- // function for exact details.
- if ((p-face->center()).norm_square() < min_distance)
+ // Look for a combination of tangent vectors that
+ // are of approximately equal length and not linearly dependent
+ for (unsigned int i=0,j=1 ; i<4 && j<4; ++j)
{
- // we are close to the face center: pick two consecutive vertices,
- // but not the closest one. We make sure the direction is always
- // the same.
- if (min_index < 2)
- {
- t1 = get_tangent_vector(p, face->vertex(3));
- t2 = get_tangent_vector(p, face->vertex(2));
- }
- else
+ // if p is too close to vertex i try again with different i and j
+ if ((p - face->vertex(i)).norm_square() <
+ std::numeric_limits<double>::epsilon() * (p - face->vertex(j)).norm_square())
{
- t1 = get_tangent_vector(p, face->vertex(0));
- t2 = get_tangent_vector(p, face->vertex(1));
- }
- }
- else
- {
- // we are closer to one of the vertices than to the
- // center of the face
- switch (min_index)
- {
- case 0:
- {
- t1 = get_tangent_vector(p, face->vertex(1));
- t2 = get_tangent_vector(p, face->vertex(2));
- break;
- }
- case 1:
- {
- t1 = get_tangent_vector(p, face->vertex(3));
- t2 = get_tangent_vector(p, face->vertex(0));
- break;
- }
- case 2:
- {
- t1 = get_tangent_vector(p, face->vertex(0));
- t2 = get_tangent_vector(p, face->vertex(3));
- break;
- }
- case 3:
- {
- t1 = get_tangent_vector(p, face->vertex(2));
- t2 = get_tangent_vector(p, face->vertex(1));
- break;
- }
- default:
- Assert(false, ExcInternalError());
- break;
+ ++i;
+ continue;
}
+
+ // if p is too close to vertex j try again with different j
+ if ((p - face->vertex(j)).norm_square() <
+ std::numeric_limits<double>::epsilon() * (p - face->vertex(i)).norm_square())
+ continue;
+
+ t1 = get_tangent_vector(p, face->vertex(i));
+ t2 = get_tangent_vector(p, face->vertex(j));
+
+ normal = cross_product_3d(t1,t2);
+
+ // if t1 and t2 are (nearly) linearly dependent try again with different j / t2
+ if (normal.norm_square() < std::numeric_limits<double>::epsilon() *
+ t1.norm_square() * t2.norm_square())
+ continue;
+
+ break;
}
- const Tensor<1,spacedim> normal = cross_product_3d(t1,t2);
+ Assert(normal.norm_square() >= std::numeric_limits<double>::epsilon() *
+ t1.norm_square() * t2.norm_square(),
+ ExcMessage("Manifold::normal_vector was unable to find a suitable combination "
+ "of vertices to compute a normal on this face. Check for distorted "
+ "faces in your triangulation."));
+
+ // Make sure all found normal vectors on this face point in the same direction
+ // as the 'reference' normal vector created at the center position.
+ const Point<spacedim> center = face->center();
+ const Tensor<1,spacedim> reference_t1 = get_tangent_vector(center, face->vertex(0));
+ const Tensor<1,spacedim> reference_t2 = get_tangent_vector(center, face->vertex(1));
+ const Tensor<1,spacedim> reference_normal = cross_product_3d(reference_t1,reference_t2);
+
+ if (reference_normal * normal < 0.0)
+ normal *= -1;
+
return normal/normal.norm();
}
--- /dev/null
+//---------------------------- manifold_id_01.cc ---------------------------
+// Copyright (C) 2011 - 2015 by the mathLab team.
+//
+// This file is subject to LGPL 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.
+//
+//---------------------------- flat_manifold_01.cc ---------------------------
+
+
+// Test that the flat manifold does what it should
+
+#include "../tests.h"
+
+
+// all include files you need here
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/grid/grid_out.h>
+
+// Helper function
+template <int dim, int spacedim>
+void test(unsigned int ref=1)
+{
+ std::vector<Point<spacedim> > vertices (GeometryInfo<dim>::vertices_per_cell);
+ std::vector<CellData<dim> > cells (1);
+ for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+ cells[0].vertices[i] = i;
+ cells[0].material_id = 0;
+
+ vertices[0] = Point<dim>(0,0,0);
+ vertices[1] = Point<dim>(1,0,0);
+ vertices[2] = Point<dim>(0.5,0.4,0);
+ vertices[3] = Point<dim>(1.5,0.4,0);
+
+ vertices[4] = Point<dim>(0,0,1);
+ vertices[5] = Point<dim>(1,0,1);
+ vertices[6] = Point<dim>(0.5,0.4,1);
+ vertices[7] = Point<dim>(1.5,0.4,1);
+
+ Triangulation<dim,spacedim> tria;
+ tria.create_triangulation (vertices, cells, SubCellData());
+
+ typename Triangulation<dim,spacedim>::active_cell_iterator
+ cell = tria.begin_active();
+
+ Point<dim> p1 (0.5,0,0);
+ deallog << "Normal vector of face 4: " << cell->get_manifold().normal_vector(cell->face(4),p1) << std::endl;
+ deallog << "Center of face 4: " << cell->face(4)->center() << std::endl;
+}
+
+int main ()
+{
+ initlog();
+
+ test<3,3>();
+
+ return 0;
+}
+
projects / dealii.git / commitdiff