fix all 2d examples including naca generated geometries

library/NekMeshUtils/2DGenerator/2DGenerator.cpp

@@ -101,7 +101,6 @@ void Generator2D::Process()
         }
 
         m_curvemeshes[i]->Mesh();
-
     }
 
     ////////////////////////////////////////
@@ -139,7 +138,7 @@ void Generator2D::Process()
 
     ////////////////////////////////////
 
-    /*EdgeSet::iterator it;
+    EdgeSet::iterator it;
     for (it = m_mesh->m_edgeSet.begin(); it != m_mesh->m_edgeSet.end(); it++)
     {
         vector<NodeSharedPtr> ns;
@@ -159,7 +158,7 @@ void Generator2D::Process()
             LibUtilities::eSegment, conf, ns, tags);
 
         m_mesh->m_element[1].push_back(E2);
-    }*/
+    }
 
     ProcessVertices();
     ProcessEdges();
@@ -200,32 +199,45 @@ void Generator2D::MakeBL(int faceid)
     for (vector<unsigned>::iterator it = m_blCurves.begin();
          it != m_blCurves.end(); ++it)
     {
-        CADSystem::Orientation edgeo = m_mesh->m_cad->GetCurve(*it)->GetOrienationWRT(faceid);
+        CADSystem::Orientation edgeo =
+            m_mesh->m_cad->GetCurve(*it)->GetOrienationWRT(faceid);
 
         vector<EdgeSharedPtr> es = m_curvemeshes[*it]->GetMeshEdges();
 
         // on each !!!EDGE!!! calculate a normal
         // always to the left unless edgeo is 1
+        // normal must be done in the parametric space (and then projected back)
+        // because of face orientation
         for (int j = 0; j < es.size(); j++)
         {
             es[j]->m_id = eid++;
             Array<OneD, NekDouble> p1, p2;
-            p1 = es[j]->m_n1->GetLoc();
-            p2 = es[j]->m_n2->GetLoc();
-            if(edgeo == CADSystem::eBackwards)
+            p1 = es[j]->m_n1->GetCADSurfInfo(faceid);
+            p2 = es[j]->m_n2->GetCADSurfInfo(faceid);
+            if (edgeo == CADSystem::eBackwards)
             {
                 swap(p1, p2);
             }
             Array<OneD, NekDouble> n(2);
-            n[0] = p1[1] - p2[1];
-            n[1] = p2[0] - p1[0];
-            if(m_mesh->m_cad->GetSurf(faceid)->IsReversedNormal())
-            {
-                n[0] *= -1.0; n[1] *= -1.0;
-            }
+            n[0]          = p1[1] - p2[1];
+            n[1]          = p2[0] - p1[0];
             NekDouble mag = sqrt(n[0] * n[0] + n[1] * n[1]);
             n[0] /= mag;
             n[1] /= mag;
+
+            Array<OneD, NekDouble> np = es[j]->m_n1->GetCADSurfInfo(faceid);
+            np[0] += n[0];
+            np[1] += n[1];
+
+            Array<OneD, NekDouble> loc  = es[j]->m_n1->GetLoc();
+            Array<OneD, NekDouble> locp = m_mesh->m_cad->GetSurf(faceid)->P(np);
+
+            n[0] = locp[0] - loc[0];
+            n[1] = locp[1] - loc[1];
+            mag  = sqrt(n[0] * n[0] + n[1] * n[1]);
+            n[0] /= mag;
+            n[1] /= mag;
+
             edgeNormals[es[j]->m_id] = n;
         }
     }
@@ -264,7 +276,8 @@ void Generator2D::MakeBL(int faceid)
     for (vector<unsigned>::iterator it = m_blCurves.begin();
          it != m_blCurves.end(); ++it)
     {
-        CADSystem::Orientation edgeo = m_mesh->m_cad->GetCurve(*it)->GetOrienationWRT(faceid);
+        CADSystem::Orientation edgeo =
+            m_mesh->m_cad->GetCurve(*it)->GetOrienationWRT(faceid);
 
         vector<NodeSharedPtr> ns = m_curvemeshes[*it]->GetMeshPoints();
         vector<NodeSharedPtr> newNs;

library/NekMeshUtils/CADSystem/CADSurf.cpp

@@ -88,8 +88,6 @@ void CADSurf::OrientateEdges(CADSurfSharedPtr surf,
         loopt.push_back(loop);
     }
 
-    cout << surf->GetId() << endl;
-
     for (int i = 0; i < loopt.size(); i++)
     {
         NekDouble area = 0.0;
@@ -108,8 +106,6 @@ void CADSurf::OrientateEdges(CADSurfSharedPtr surf,
         area = bg::area(polygon);
 
         ein[i]->area = area;
-
-        cout << area << endl;
     }
 
     // order by absoulte area
@@ -130,49 +126,6 @@ void CADSurf::OrientateEdges(CADSurfSharedPtr surf,
 
     } while (ct > 0);
 
-    /*if (ein[0]->area < 0) // reverse the first uvLoop
-    {
-        ein[0]->area *= -1.0;
-        reverse(ein[0]->edgeo.begin(), ein[0]->edgeo.end());
-        reverse(ein[0]->edges.begin(), ein[0]->edges.end());
-        reverse(loopt[0].begin(), loopt[0].end());
-        // need to flip edgeo
-        for (int i = 0; i < ein[0]->edgeo.size(); i++)
-        {
-            if (ein[0]->edgeo[i] == CADSystem::eForwards)
-            {
-                ein[0]->edgeo[i] = CADSystem::eBackwards;
-            }
-            else
-            {
-                ein[0]->edgeo[i] = CADSystem::eForwards;
-            }
-        }
-    }
-
-    for (int i = 1; i < ein.size(); i++)
-    {
-        if (ein[i]->area > 0) // reverse the loop
-        {
-            ein[i]->area *= -1.0;
-            reverse(ein[i]->edgeo.begin(), ein[i]->edgeo.end());
-            reverse(ein[i]->edges.begin(), ein[i]->edges.end());
-            reverse(loopt[i].begin(), loopt[i].end());
-            // need to flip edgeo
-            for (int j = 0; j < ein[i]->edgeo.size(); j++)
-            {
-                if (ein[i]->edgeo[j] == CADSystem::eForwards)
-                {
-                    ein[i]->edgeo[j] = CADSystem::eBackwards;
-                }
-                else
-                {
-                    ein[i]->edgeo[j] = CADSystem::eForwards;
-                }
-            }
-        }
-    }*/
-
     //only need center points for inner loops
     for(int i = 1; i < ein.size(); i++)
     {

library/NekMeshUtils/CADSystem/CADSurf.h

@@ -60,7 +60,7 @@ public:
     /**
      * @brief Default constructor.
      */
-    CADSurf() : m_correctNormal(true)
+    CADSurf() : m_orient(CADSystem::eForwards)
     {
         m_type = CADType::eSurf;
     }
@@ -155,24 +155,13 @@ public:
     /**
      * @brief query reversed normal
      */
-    bool IsReversedNormal()
+    CADSystem::Orientation Orientation()
     {
-        return m_correctNormal;
+        return m_orient;
     }
 
 protected:
-    /**
-     * @brief sets the flag to reverse the normal for this suface,
-     * this is determined in cadsystem and ensures all surface normals,
-     * point internaly
-     */
-    void SetReverseNomral()
-    {
-        m_correctNormal = false;
-    }
-
-    /// normal
-    bool m_correctNormal;
+    CADSystem::Orientation m_orient;
     /// List of bounding edges in loops with orientation.
     std::vector<CADSystem::EdgeLoopSharedPtr> m_edges;
 

library/NekMeshUtils/CADSystem/OCE/CADSurfOCE.cpp

@@ -51,11 +51,10 @@ void CADSurfOCE::Initialise(int i, TopoDS_Shape in)
     // defualt to m
 
     m_s = BRep_Tool::Surface(TopoDS::Face(in));
-    //reverse the face
 
     if(in.Orientation() == 1)
     {
-        SetReverseNomral();
+        m_orient = CADSystem::eBackwards;
     }
 
     gp_Trsf transform;
@@ -65,12 +64,7 @@ void CADSurfOCE::Initialise(int i, TopoDS_Shape in)
 
     in.Move(mv);
     m_occSurface    = BRepAdaptor_Surface(TopoDS::Face(in));
-    m_correctNormal = true;
     m_id            = i;
-
-    Array<OneD, NekDouble> uv(2);
-    uv[0] = GetBounds()[0];
-    uv[1] = GetBounds()[2];
 }
 
 Array<OneD, NekDouble> CADSurfOCE::GetBounds()
@@ -275,7 +269,7 @@ Array<OneD, NekDouble> CADSurfOCE::N(Array<OneD, NekDouble> uv)
 
     Array<OneD, NekDouble> normal(3);
 
-    if (!m_correctNormal)
+    if (m_orient == CADSystem::eBackwards)
     {
         d.Reverse();
     }

library/NekMeshUtils/CADSystem/OCE/CADSystemOCE.cpp

@@ -78,7 +78,7 @@ bool CADSystemOCE::LoadCAD()
     for (explr.Init(shape, TopAbs_VERTEX); explr.More(); explr.Next())
     {
         TopoDS_Shape v = explr.Current();
-        if(mapOfVerts.Contains(v))
+        if (mapOfVerts.Contains(v))
         {
             continue;
         }
@@ -92,7 +92,7 @@ bool CADSystemOCE::LoadCAD()
     for (explr.Init(shape, TopAbs_EDGE); explr.More(); explr.Next())
     {
         TopoDS_Shape e = explr.Current().Oriented(TopAbs_FORWARD);
-        if(mapOfEdges.Contains(e))
+        if (mapOfEdges.Contains(e))
         {
             continue;
         }
@@ -100,7 +100,7 @@ bool CADSystemOCE::LoadCAD()
         if (curve.GetType() != 7)
         {
             int i = mapOfEdges.Add(e);
-            AddCurve(i,e);
+            AddCurve(i, e);
         }
     }
 
@@ -183,7 +183,7 @@ void CADSystemOCE::AddSurf(int i, TopoDS_Shape in)
     CADSurfSharedPtr newSurf = GetCADSurfFactory().CreateInstance(key);
     boost::static_pointer_cast<CADSurfOCE>(newSurf)->Initialise(i, in);
 
-    //do the exploration on forward oriented
+    // do the exploration on forward oriented
     TopoDS_Shape face = in.Oriented(TopAbs_FORWARD);
     TopTools_IndexedMapOfShape mapOfWires;
     TopExp::MapShapes(face, TopAbs_WIRE, mapOfWires);
@@ -207,9 +207,9 @@ void CADSystemOCE::AddSurf(int i, TopoDS_Shape in)
             {
                 int e = mapOfEdges.FindIndex(edge);
                 edgeloop->edges.push_back(m_curves[e]);
-                edgeloop->edgeo.push_back(
-                    exp.Orientation() == TopAbs_FORWARD ? eForwards
-                                                        : eBackwards);
+                edgeloop->edgeo.push_back(exp.Orientation() == TopAbs_FORWARD
+                                              ? eForwards
+                                              : eBackwards);
             }
             exp.Next();
         }
@@ -232,7 +232,8 @@ void CADSystemOCE::AddSurf(int i, TopoDS_Shape in)
     {
         for (int j = 0; j < edgeloops[k]->edges.size(); j++)
         {
-            edgeloops[k]->edges[j]->SetAdjSurf(make_pair(newSurf, edgeloops[k]->edgeo[j]));
+            edgeloops[k]->edges[j]->SetAdjSurf(
+                make_pair(newSurf, edgeloops[k]->edgeo[j]));
         }
     }
 
@@ -384,7 +385,11 @@ TopoDS_Shape CADSystemOCE::BuildNACA(string naca)
     BRepBuilderAPI_MakeFace face(domWire, true);
     face.Add(aeroWire);
 
-    return face.Face();
+    ShapeFix_Face sf(face.Face());
+    sf.FixOrientation();
+
+
+    return sf.Face();
 }
 }
 }

library/NekMeshUtils/CADSystem/OCE/OpenCascade.h

@@ -71,6 +71,8 @@
 #include <ShapeAnalysis_Curve.hxx>
 #include <BRepLProp_SLProps.hxx>
 #include <Standard_Macro.hxx>
+#include <ShapeFix_Face.hxx>
+
 
 #include <GeomAPI_PointsToBSpline.hxx>
 #include <Geom_BSplineCurve.hxx>

library/NekMeshUtils/MeshElements/Element.h

@@ -358,13 +358,7 @@ public:
             }
         }
     }
-
-    NEKMESHUTILS_EXPORT void Flip()
-    {
-        ASSERTL0(m_face.size() == 0,"only works in 2D")
-        std::swap(m_vertex[0],m_vertex[1]);
-    }
-
+    
     CADObjectSharedPtr m_parentCAD;
 
 protected:

library/NekMeshUtils/SurfaceMeshing/FaceMesh.cpp

@@ -891,11 +891,6 @@ void FaceMesh::BuildLocalMesh()
     {
         ElmtConfig conf(LibUtilities::eTriangle, 1, false, false);
 
-        if(!m_cadsurf->IsReversedNormal())
-        {
-            swap(m_connec[i][0],m_connec[i][1]);
-        }
-
         vector<int> tags;
         tags.push_back(m_compId);
         ElementSharedPtr E = GetElementFactory().CreateInstance(

library/NekMeshUtils/VolumeMeshing/BLMeshing/BLMesh.cpp

@@ -1030,14 +1030,6 @@ void BLMesh::Setup()
     // but do not care about triangles which are not in the bl
     for (int i = 0; i < m_mesh->m_element[2].size(); i++)
     {
-        // orientate the triangle
-        if (!m_mesh->m_cad
-                 ->GetSurf(m_mesh->m_element[2][i]->m_parentCAD->GetId())
-                 ->IsReversedNormal())
-        {
-            m_mesh->m_element[2][i]->Flip();
-        }
-
         vector<unsigned int>::iterator f =
             find(m_blsurfs.begin(), m_blsurfs.end(),
                  m_mesh->m_element[2][i]->m_parentCAD->GetId());

utilities/NekMesh/InputModules/InputMCF.cpp

@@ -333,7 +333,7 @@ void InputMCF::Process()
     }
 
     ////*** VARIATIONAL OPTIMISATION ****////
-    /*if(m_varopti)
+    if(m_varopti)
     {
         unsigned int np = boost::thread::physical_concurrency();
         if(m_mesh->m_verbose)