const Point<3> &p) const
{
const int spacedim=3;
- Tensor<1,spacedim> t1,t2,normal;
- // the counter-clockwise sequence of face vertices (see
- // the GeometryInfo class for more documentation).
- const double face_vertices[4] = {0,1,3,2};
+ // Compute distances from p to vertices
+ std::array<std::pair<double,unsigned int>, GeometryInfo<3>::vertices_per_face> distances;
+ for (unsigned int i=0; i<GeometryInfo<3>::vertices_per_face; ++i)
+ distances[i] = std::make_pair((p-face->vertex(i)).norm_square(),i);
- // Look for the vertex with the largest distance to the given point
- unsigned int max_index=0;
- double max_distance = (p-face->vertex(face_vertices[0])).norm_square();
+ // Sort the distances and figure out, which vertices to use
+ std::sort(distances.begin(),distances.end());
+ unsigned int first_index, second_index, fallback_index;
- for (unsigned int i=1; i<GeometryInfo<3>::vertices_per_face; ++i)
+ const double center_distance = (p-face->center()).norm_square();
+
+ // If we are closer to the center than any vertex, use closest vertices
+ if (center_distance < distances[0].first)
{
- const double distance = (p-face->vertex(face_vertices[i])).norm_square();
- if (distance > max_distance)
- {
- max_index = i;
- max_distance = distance;
- }
+ first_index = distances[0].second;
+ second_index = distances[1].second;
+ fallback_index = distances[2].second;
}
-
- // Compute tangent for max_index vertex.
- t1 = get_tangent_vector(p, face->vertex(face_vertices[max_index]));
-
- // If p is sufficiently far away from next counter-clockwise vertex, compute normal.
- const unsigned int next_index = (max_index + 1) % 4;
-
- if ((p-face->vertex(face_vertices[next_index])).norm_square() >
- 1e4 * std::numeric_limits<double>::epsilon() * max_distance)
+ // Otherwise use vertices further away
+ else
{
- t2 = get_tangent_vector(p, face->vertex(face_vertices[next_index]));
- normal = cross_product_3d(t1,t2);
+ first_index = distances[1].second;
+ second_index = distances[2].second;
+
+ // If we are very close to a vertex use the farthest vertex as fallback,
+ // otherwise use closest vertex.
+ if (distances[0].first < 1e4 * std::numeric_limits<double>::epsilon() * distances[1].first)
+ fallback_index = distances[3].second;
+ else
+ fallback_index = distances[0].second;
}
- // If p is too close to next_index, or the tangents are linearly dependent
- // fall back to previous counter clockwise vertex.
- // In this case the normal direction is flipped.
- if (normal.norm_square() < 1e4 * std::numeric_limits<double>::epsilon() *
+ // Compute tangents and normal for selected vertices
+ Tensor<1,spacedim> t1 = get_tangent_vector(p, face->vertex(first_index));
+ Tensor<1,spacedim> t2 = get_tangent_vector(p, face->vertex(second_index));
+ Tensor<1,spacedim> normal = cross_product_3d(t1,t2);
+
+ // If the tangents are linearly dependent fall back to another tangent
+ if (normal.norm_square() <= 1e4 * std::numeric_limits<double>::epsilon() *
t1.norm_square() * t2.norm_square())
{
- const unsigned int previous_index = (max_index - 1) % 4;
- t2 = get_tangent_vector(p, face->vertex(face_vertices[previous_index]));
- normal = cross_product_3d(t2,t1);
+ t2 = get_tangent_vector(p, face->vertex(fallback_index));
+ normal = cross_product_3d(t1,t2);
}
- Assert(normal.norm_square() >= 1e4 * std::numeric_limits<double>::epsilon() *
+ Assert(normal.norm_square() > 1e4 * 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."));
+ // Now figure out if we need to flip the direction, we do this by comparing to a reference
+ // normal that would be the correct result if all vertices would lie in a plane
+ const Tensor<1,spacedim> rt1 = face->vertex(3) - face->vertex(0);
+ const Tensor<1,spacedim> rt2 = face->vertex(2) - face->vertex(1);
+ const Tensor<1,spacedim> reference_normal = cross_product_3d(rt1,rt2);
+
+ if (reference_normal * normal < 0.0)
+ normal *= -1.0;
+
return normal / normal.norm();
}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2000 - 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Copy of no_flux_09 for a higher order element. The purpose here is to test the
+// Manifold::normal_vector function inside of VectorTools::compute_no_normal_flux_constraints
+// for points other than the vertices.
+
+
+#include "../tests.h"
+
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/numerics/vector_tools.h>
+
+
+template <int dim>
+void
+check ()
+{
+ Triangulation<dim> tr;
+ GridGenerator::quarter_hyper_shell (tr,
+ Point<dim>(),
+ 0.5, 1.0,
+ 3, true);
+
+ ConstraintMatrix cm;
+ MappingQ<dim> mapping(1);
+
+ FESystem<dim> fe(FE_Q<dim>(3),dim);
+ DoFHandler<dim> dofh(tr);
+
+ dofh.distribute_dofs (fe);
+
+ std::set<types::boundary_id> no_normal_flux_boundaries;
+ no_normal_flux_boundaries.insert (1);
+ // no_normal_flux_boundaries.insert (2); // not required for the crash for now, please test with it later!
+ no_normal_flux_boundaries.insert (3);
+ no_normal_flux_boundaries.insert (4);
+ VectorTools::compute_no_normal_flux_constraints (dofh, 0, no_normal_flux_boundaries, cm, mapping);
+
+ cm.print (deallog.get_file_stream ());
+}
+
+
+
+int main ()
+{
+ std::ofstream logfile ("output");
+ logfile.precision (4);
+ logfile.setf(std::ios::fixed);
+ deallog.attach(logfile);
+
+ check<3> ();
+}