Fix bug in normal to mesh projection due to Boundary -> Manifold

diff --git a/source/opencascade/boundary_lib.cc b/source/opencascade/boundary_lib.cc
index 1a3df42f5c26b5dc245d4c065179756e9f35612d..773a30b7c0dd03712bcd01e75c16bd23e40be23d 100644 (file)
--- a/source/opencascade/boundary_lib.cc
+++ b/source/opencascade/boundary_lib.cc
@@ -189,6 +189,27 @@ namespace OpenCASCADE
         average_normal /= average_normal.norm();
         break;
       }
+      case 4:
+      {
+        Tensor<1,3> u = surrounding_points[1]-surrounding_points[0];
+        Tensor<1,3> v = surrounding_points[2]-surrounding_points[0];
+        const double n1_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(n1_coords);
+        n1 = n1/n1.norm();
+        u = surrounding_points[2]-surrounding_points[3];
+        v = surrounding_points[1]-surrounding_points[3];
+        const double n2_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> n2(n2_coords);
+        n2 = n2/n2.norm();
+
+        average_normal = (n1+n2)/2.0;
+
+        Assert(average_normal.norm() > tolerance,
+               ExcMessage("Failed to refine cell: the normal estimated via the surrounding points turns out to be a null vector, making the projection direction undetermined."));
+
+        average_normal /= average_normal.norm();
+        break;
+      }
       case 8:
       {
         Tensor<1,3> u = surrounding_points[1]-surrounding_points[0];
@@ -211,10 +232,6 @@ namespace OpenCASCADE
         const double n4_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> n4(n4_coords);
         n4 = n4/n4.norm();
-        //for (unsigned int i=0; i<surrounding_points.size(); ++i)
-        //    cout<<surrounding_points[i]<<endl;
-        //cout<<"-"<<endl;
-        //cout<<n1<<endl;cout<<n2<<endl;cout<<n3<<endl;cout<<n4<<endl;
 
         average_normal = (n1+n2+n3+n4)/4.0;
 

diff --git a/tests/opencascade/normal_to_mesh_projection_03.cc b/tests/opencascade/normal_to_mesh_projection_03.cc
new file mode 100644 (file)
index 0000000..410e78e
--- /dev/null
+++ b/tests/opencascade/normal_to_mesh_projection_03.cc
@@ -0,0 +1,107 @@
+//-----------------------------------------------------------
+//
+//    Copyright (C) 2014 - 2017 by the deal.II authors
+//
+//    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.
+//
+//-----------------------------------------------------------
+
+// Create a BSpline surface, and test axis projection for Q2 elements.
+
+#include "../tests.h"
+
+#include <deal.II/opencascade/utilities.h>
+#include <deal.II/opencascade/boundary_lib.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_generator.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/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <TopTools.hxx>
+#include <TopoDS_Shape.hxx>
+#include <TopoDS_Edge.hxx>
+#include <TopoDS_Face.hxx>
+#include <BRepFill.hxx>
+#include <Standard_Stream.hxx>
+
+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>(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);
+
+
+  NormalToMeshProjectionBoundary<2,3> manifold(face);
+
+  Triangulation<2,3> tria;
+  GridGenerator::hyper_cube(tria);
+
+  tria.begin_active()->set_all_manifold_ids(1);
+  tria.set_manifold(1, manifold);
+
+  tria.refine_global(1);
+
+  Tensor<1,3> direction({.1, .01, .001});
+  double tolerance = 1e-10;
+
+  {
+    FE_Q<2,3> fe(2);
+    DoFHandler<2,3> dh(tria);
+    dh.distribute_dofs(fe);
+    MappingQGeneric<2,3> mapping2(2);
+    std::vector<Point<3>> spoints(dh.n_dofs());
+    DoFTools::map_dofs_to_support_points(mapping2, dh, spoints);
+
+    std::sort(spoints.begin(), spoints.end(),
+              [&](const Point<3> &p1, const Point<3> &p2)
+    {
+      return OpenCASCADE::point_compare(p1, p2, direction, tolerance);
+    });
+
+    for (auto p: spoints)
+      deallog << p << std::endl;
+    deallog << "============" << std::endl;
+  }
+
+  tria.refine_global(1);
+  {
+    FE_Q<2,3> fe(1);
+    DoFHandler<2,3> dh(tria);
+    dh.distribute_dofs(fe);
+    std::vector<Point<3>> spoints(dh.n_dofs());
+    DoFTools::map_dofs_to_support_points(StaticMappingQ1<2,3>::mapping,
+                                         dh, spoints);
+
+    std::sort(spoints.begin(), spoints.end(),
+              [&](const Point<3> &p1, const Point<3> &p2)
+    {
+      return OpenCASCADE::point_compare(p1, p2, direction, tolerance);
+    });
+
+    for (auto p: spoints)
+      deallog << p << std::endl;
+  }
+}

diff --git a/tests/opencascade/normal_to_mesh_projection_03.output b/tests/opencascade/normal_to_mesh_projection_03.output
new file mode 100644 (file)
index 0000000..87f13a9
--- /dev/null
+++ b/tests/opencascade/normal_to_mesh_projection_03.output
@@ -0,0 +1,52 @@
+
+DEAL::0.00000 0.00000 0.00000
+DEAL::0.00000 0.250000 0.00000
+DEAL::0.00000 0.500000 0.00000
+DEAL::0.00000 0.750000 0.00000
+DEAL::0.00000 1.00000 0.00000
+DEAL::0.234806 0.00000 0.153640
+DEAL::0.234806 0.250000 0.153640
+DEAL::0.234806 0.500000 0.153640
+DEAL::0.234806 0.750000 0.153640
+DEAL::0.234806 1.00000 0.153640
+DEAL::0.500000 0.00000 0.245356
+DEAL::0.500000 0.250000 0.245356
+DEAL::0.500000 0.500000 0.245356
+DEAL::0.500000 0.750000 0.245356
+DEAL::0.500000 1.00000 0.245356
+DEAL::0.789901 0.00000 0.203990
+DEAL::0.789901 0.250000 0.203990
+DEAL::0.789901 0.500000 0.203990
+DEAL::0.789901 0.750000 0.203990
+DEAL::0.789901 1.00000 0.203990
+DEAL::1.00000 0.00000 0.00000
+DEAL::1.00000 0.250000 -2.79105e-17
+DEAL::1.00000 0.500000 -4.62817e-17
+DEAL::1.00000 0.750000 -3.17705e-17
+DEAL::1.00000 1.00000 0.00000
+DEAL::============
+DEAL::0.00000 0.00000 0.00000
+DEAL::0.00000 0.250000 0.00000
+DEAL::0.00000 0.500000 0.00000
+DEAL::0.00000 0.750000 0.00000
+DEAL::0.00000 1.00000 0.00000
+DEAL::0.234806 0.00000 0.153640
+DEAL::0.234806 0.250000 0.153640
+DEAL::0.234806 0.500000 0.153640
+DEAL::0.234806 0.750000 0.153640
+DEAL::0.234806 1.00000 0.153640
+DEAL::0.500000 0.00000 0.245356
+DEAL::0.500000 0.250000 0.245356
+DEAL::0.500000 0.500000 0.245356
+DEAL::0.500000 0.750000 0.245356
+DEAL::0.500000 1.00000 0.245356
+DEAL::0.789901 0.00000 0.203990
+DEAL::0.789901 0.250000 0.203990
+DEAL::0.789901 0.500000 0.203990
+DEAL::0.789901 0.750000 0.203990
+DEAL::0.789901 1.00000 0.203990
+DEAL::1.00000 0.00000 0.00000
+DEAL::1.00000 0.250000 -2.79105e-17
+DEAL::1.00000 0.500000 -4.62817e-17
+DEAL::1.00000 0.750000 -3.17705e-17
+DEAL::1.00000 1.00000 0.00000