From faeece06687430cb8240ea62064a932b5b4f8047 Mon Sep 17 00:00:00 2001 From: Rene Gassmoeller Date: Mon, 6 Nov 2017 14:23:03 -0700 Subject: [PATCH] Test third algorithm --- source/grid/manifold.cc | 77 +++++++++++--------- tests/numerics/no_flux_14.cc | 72 +++++++++++++++++++ tests/numerics/no_flux_14.output | 118 +++++++++++++++++++++++++++++++ 3 files changed, 234 insertions(+), 33 deletions(-) create mode 100644 tests/numerics/no_flux_14.cc create mode 100644 tests/numerics/no_flux_14.output diff --git a/source/grid/manifold.cc b/source/grid/manifold.cc index 8653fc26c1..67c8c4ce81 100644 --- a/source/grid/manifold.cc +++ b/source/grid/manifold.cc @@ -156,56 +156,67 @@ normal_vector (const Triangulation<3, 3>::face_iterator &face, 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, GeometryInfo<3>::vertices_per_face> distances; + for (unsigned int i=0; i::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::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::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::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::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::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::epsilon() * + Assert(normal.norm_square() > 1e4 * std::numeric_limits::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(); } diff --git a/tests/numerics/no_flux_14.cc b/tests/numerics/no_flux_14.cc new file mode 100644 index 0000000000..bfe007a2dc --- /dev/null +++ b/tests/numerics/no_flux_14.cc @@ -0,0 +1,72 @@ +// --------------------------------------------------------------------- +// +// 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 +#include +#include +#include +#include +#include +#include + + +template +void +check () +{ + Triangulation tr; + GridGenerator::quarter_hyper_shell (tr, + Point(), + 0.5, 1.0, + 3, true); + + ConstraintMatrix cm; + MappingQ mapping(1); + + FESystem fe(FE_Q(3),dim); + DoFHandler dofh(tr); + + dofh.distribute_dofs (fe); + + std::set 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> (); +} diff --git a/tests/numerics/no_flux_14.output b/tests/numerics/no_flux_14.output new file mode 100644 index 0000000000..660e9be9d6 --- /dev/null +++ b/tests/numerics/no_flux_14.output @@ -0,0 +1,118 @@ + + 2 = 0 + 5 = 0 + 7 6: -0.4142 + 8 = 0 + 10 9: -1.0000 + 11 = 0 + 19 18: -0.4889 + 19 20: -0.9881 + 22 21: -1.0000 + 22 23: -0.8648 + 28 = 0 + 29 = 0 + 34 = 0 + 35 = 0 + 40 = 0 + 41 = 0 + 44 42: -0.4142 + 46 = 0 + 45 43: -0.4142 + 47 = 0 + 68 66: -0.4889 + 68 70: -0.9881 + 69 67: -0.4394 + 69 71: -0.8417 + 86 84: -0.2948 + 86 88: -0.4142 + 87 85: -0.2948 + 87 89: -0.4142 + 92 90: -1.0000 + 92 94: -0.3789 + 93 91: -1.0000 + 93 95: -0.3789 + 136 132: -0.4142 + 136 140: -0.3831 + 137 133: -0.2825 + 137 141: -0.3778 + 138 134: -0.3716 + 138 142: -0.2599 + 139 135: -0.2916 + 139 143: -0.2447 + 152 = 0 + 153 = 0 + 154 = 0 + 155 = 0 + 193 = 0 + 194 = 0 + 195 = 0 + 196 = 0 + 197 = 0 + 199 = 0 + 202 = 0 + 203 201: -1.0000 + 206 = 0 + 208 = 0 + 207 = 0 + 209 = 0 + 214 = 0 + 215 = 0 + 216 218: -0.4142 + 220 = 0 + 217 219: -0.4142 + 221 = 0 + 224 = 0 + 225 = 0 + 234 236: -0.4394 + 234 238: -0.8417 + 235 237: -0.4889 + 235 239: -0.9881 + 242 = 0 + 243 = 0 + 246 250: -0.4142 + 248 = 0 + 247 251: -0.4142 + 249 = 0 + 256 = 0 + 257 = 0 + 258 = 0 + 259 = 0 + 276 280: -0.3716 + 276 284: -0.2599 + 277 281: -0.2916 + 277 285: -0.2447 + 278 282: -0.4142 + 278 286: -0.3831 + 279 283: -0.2825 + 279 287: -0.3778 + 296 = 0 + 297 = 0 + 298 = 0 + 299 = 0 + 337 = 0 + 340 = 0 + 341 339: -0.4142 + 344 = 0 + 345 = 0 + 350 = 0 + 351 = 0 + 356 = 0 + 358 354: -0.4142 + 357 = 0 + 359 355: -0.4142 + 370 366: -0.4023 + 370 368: -0.4142 + 371 367: -0.4142 + 371 369: -0.4142 + 404 396: -0.4018 + 404 400: -0.4108 + 405 397: -0.4112 + 405 401: -0.4142 + 406 398: -0.4008 + 406 402: -0.4045 + 407 399: -0.4108 + 407 403: -0.4018 + 412 = 0 + 413 = 0 + 414 = 0 + 415 = 0 -- 2.39.5