added default case for normal to mesh projection manifold

 and indent

diff --git a/source/opencascade/manifold_lib.cc b/source/opencascade/manifold_lib.cc
index dd873d0d15a213d22119184047c2182eb95e4878..701eaada993c61a952a661e150137864fde6748d 100644 (file)
--- a/source/opencascade/manifold_lib.cc
+++ b/source/opencascade/manifold_lib.cc
@@ -302,7 +302,45 @@ namespace OpenCASCADE
           }
         default:
           {
-            AssertThrow(false, ExcNotImplemented());
+            // Given an arbitrary number of points we compute all the possible
+            // normal vectors
+            for (unsigned int i = 0; i < surrounding_points.size(); ++i)
+              for (unsigned int j = 0; j < surrounding_points.size(); ++j)
+                if (j != i)
+                  for (unsigned int k = 0; k < surrounding_points.size(); ++k)
+                    if (k != j && k != i)
+                      {
+                        Tensor<1, 3> u =
+                          surrounding_points[i] - surrounding_points[j];
+                        Tensor<1, 3> v =
+                          surrounding_points[i] - surrounding_points[k];
+                        const double n_coords[3] = {u[1] * v[2] - u[2] * v[1],
+                                                    u[2] * v[0] - u[0] * v[2],
+                                                    u[0] * v[1] - u[1] * v[0]};
+                        Tensor<1, 3> n1(n_coords);
+                        if (n1.norm() > tolerance)
+                          {
+                            n1 = n1 / n1.norm();
+                            if (average_normal.norm() < tolerance)
+                              average_normal = n1;
+                            else
+                              {
+                                auto dot_prod = n1 * average_normal;
+                                // We check that the direction of the normal
+                                // vector w.r.t the current average, and make
+                                // sure we flip it if it is opposite
+                                if (dot_prod > 0)
+                                  average_normal += n1;
+                                else
+                                  average_normal -= n1;
+                              }
+                          }
+                      }
+            Assert(
+              average_normal.norm() > tolerance,
+              ExcMessage(
+                "Failed to compute a normal: the normal estimated via the surrounding points turns out to be a null vector, making the projection direction undetermined."));
+            average_normal = average_normal / average_normal.norm();
             break;
           }
       }

diff --git a/tests/opencascade/normal_to_mesh_projection_04.cc b/tests/opencascade/normal_to_mesh_projection_04.cc
new file mode 100644 (file)
index 0000000..9f303ab
--- /dev/null
+++ b/tests/opencascade/normal_to_mesh_projection_04.cc
@@ -0,0 +1,101 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+// Create a BSpline surface, and test NormalToMeshProjectionManifold  projection
+// using 5 surrounding points.
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_q_generic.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/opencascade/manifold_lib.h>
+#include <deal.II/opencascade/utilities.h>
+
+#include <BRepFill.hxx>
+#include <Standard_Stream.hxx>
+#include <TopTools.hxx>
+#include <TopoDS_Edge.hxx>
+#include <TopoDS_Face.hxx>
+#include <TopoDS_Shape.hxx>
+
+#include "../tests.h"
+
+using namespace OpenCASCADE;
+
+int
+main()
+{
+  initlog();
+
+  // Create a bspline passing through the points
+  std::vector<Point<3>> pts;
+  pts.push_back(Point<3>(0, 0, 0));
+  pts.push_back(Point<3>(1, 0, 0));
+  // pts.push_back(Point<3>(0.5, 2, -0.5));
+  pts.push_back(Point<3>(1, 0, -1));
+  pts.push_back(Point<3>(0, 0, -1));
+
+  TopoDS_Edge edge1 = interpolation_curve(pts);
+  for (unsigned int i = 0; i < pts.size(); ++i)
+    pts[i] += Point<3>(0, 1, 0);
+  TopoDS_Edge edge2 = interpolation_curve(pts);
+
+  TopoDS_Face face = BRepFill::Face(edge1, edge2);
+
+  NormalToMeshProjectionManifold<2, 3> manifold(face);
+  FlatManifold<2, 3>                   flat_manifold;
+  std::vector<Point<3>>                surrounding_points;
+  surrounding_points.push_back(Point<3>(0., 0., 0.));
+  surrounding_points.push_back(Point<3>(1., 0., 0.));
+  Vector<double> weights(surrounding_points.size());
+  weights = 1. / surrounding_points.size();
+  auto intermediate_point =
+    manifold.get_new_point(ArrayView<Point<3>>(&surrounding_points[0],
+                                               surrounding_points.size()),
+                           ArrayView<double>(&weights[0], weights.size()));
+  surrounding_points.push_back(intermediate_point);
+  surrounding_points.push_back(Point<3>(0., 1., 0.));
+  surrounding_points.push_back(Point<3>(1., 1., 0.));
+  weights.reinit(surrounding_points.size());
+  weights = 1. / surrounding_points.size();
+
+
+  auto new_point =
+    manifold.get_new_point(ArrayView<Point<3>>(&surrounding_points[0],
+                                               surrounding_points.size()),
+                           ArrayView<double>(&weights[0], weights.size()));
+  auto new_point_flat =
+    flat_manifold.get_new_point(ArrayView<Point<3>>(&surrounding_points[0],
+                                                    surrounding_points.size()),
+                                ArrayView<double>(&weights[0], weights.size()));
+
+  for (unsigned int i = 0; i < surrounding_points.size(); ++i)
+    {
+      deallog << "Surrunding point " << i << " " << surrounding_points[i]
+              << " weight " << weights[i] << std::endl;
+    }
+  deallog << "Mean point " << new_point_flat << std::endl;
+  deallog << "Projected point "
+          << " " << new_point << std::endl;
+}

diff --git a/tests/opencascade/normal_to_mesh_projection_04.output b/tests/opencascade/normal_to_mesh_projection_04.output
new file mode 100644 (file)
index 0000000..0a8857d
--- /dev/null
+++ b/tests/opencascade/normal_to_mesh_projection_04.output
@@ -0,0 +1,8 @@
+
+DEAL::Surrunding point 0 0.00000 0.00000 0.00000 weight 0.200000
+DEAL::Surrunding point 1 1.00000 0.00000 0.00000 weight 0.200000
+DEAL::Surrunding point 2 0.500000 0.00000 0.245356 weight 0.200000
+DEAL::Surrunding point 3 0.00000 1.00000 0.00000 weight 0.200000
+DEAL::Surrunding point 4 1.00000 1.00000 0.00000 weight 0.200000
+DEAL::Mean point 0.500000 0.400000 0.0490712
+DEAL::Projected point  0.500000 0.447721 0.245356