]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Test third algorithm 5393/head
authorRene Gassmoeller <rene.gassmoeller@mailbox.org>
Mon, 6 Nov 2017 21:23:03 +0000 (14:23 -0700)
committerRene Gassmoeller <rene.gassmoeller@mailbox.org>
Tue, 7 Nov 2017 15:55:45 +0000 (08:55 -0700)
source/grid/manifold.cc
tests/numerics/no_flux_14.cc [new file with mode: 0644]
tests/numerics/no_flux_14.output [new file with mode: 0644]

index 8653fc26c15a097534a0ec5624bada9fb334990e..67c8c4ce81e0a7a6b9f1744777825f3b4b9a748d 100644 (file)
@@ -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<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();
 }
 
diff --git a/tests/numerics/no_flux_14.cc b/tests/numerics/no_flux_14.cc
new file mode 100644 (file)
index 0000000..bfe007a
--- /dev/null
@@ -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 <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> ();
+}
diff --git a/tests/numerics/no_flux_14.output b/tests/numerics/no_flux_14.output
new file mode 100644 (file)
index 0000000..660e9be
--- /dev/null
@@ -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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.