inital rework

library/NekMeshUtils/CADSystem/CADCurve.h

@@ -71,7 +71,7 @@ public:
      *
      * @return Array of two entries, min and max parametric coordinate.
      */
-    virtual Array<OneD, NekDouble> Bounds() = 0;
+    virtual Array<OneD, NekDouble> GetBounds() = 0;
 
     /**
      * @brief Calculates the arclength between the two paremetric points \p ti
@@ -120,15 +120,15 @@ public:
     /**
      * @brief set the ids of the surfaces either side of the curve
      */
-    void SetAdjSurf(std::vector<CADSurfSharedPtr> i)
+    void SetAdjSurf(std::pair<CADSurfSharedPtr, Orientation> i)
     {
-        m_adjSurfs = i;
+        m_adjSurfs.push_back(i);
     }
 
     /*
      * @brief returns the ids of neigbouring surfaces
      */
-    std::vector<CADSurfSharedPtr> GetAdjSurf()
+    std::vector<std::pair<CADSurfSharedPtr, Orientation> > GetAdjSurf()
     {
         return m_adjSurfs;
     }
@@ -168,7 +168,7 @@ protected:
     /// Length of edge
     NekDouble m_length;
     /// List of surfaces which this curve belongs to.
-    std::vector<CADSurfSharedPtr> m_adjSurfs;
+    std::vector<std::pair<CADSurfSharedPtr, Orientation> > m_adjSurfs;
     /// list of end vertices
     std::vector<CADVertSharedPtr> m_mainVerts;
 };

library/NekMeshUtils/CADSystem/CADSurf.cpp

+////////////////////////////////////////////////////////////////////////////////
+//
+//  File: CADSystem.cpp
+//
+//  For more information, please see: http://www.nektar.info/
+//
+//  The MIT License
+//
+//  Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
+//  Department of Aeronautics, Imperial College London (UK), and Scientific
+//  Computing and Imaging Institute, University of Utah (USA).
+//
+//  License for the specific language governing rights and limitations under
+//  Permission is hereby granted, free of charge, to any person obtaining a
+//  copy of this software and associated documentation files (the "Software"),
+//  to deal in the Software without restriction, including without limitation
+//  the rights to use, copy, modify, merge, publish, distribute, sublicense,
+//  and/or sell copies of the Software, and to permit persons to whom the
+//  Software is furnished to do so, subject to the following conditions:
+//
+//  The above copyright notice and this permission notice shall be included
+//  in all copies or substantial portions of the Software.
+//
+//  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+//  OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+//  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+//  THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+//  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+//  FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+//  DEALINGS IN THE SOFTWARE.
+//
+//  Description: cad object methods.
+//
+////////////////////////////////////////////////////////////////////////////////
+
+#include "CADSurf.h"
+#include "CADCurve.h"
+
+using namespace std;
+
+namespace Nektar
+{
+namespace NekMeshUtils
+{
+
+void CADSurf::OrientateEdges(CADSurfSharedPtr surf,
+                             vector<EdgeLoopSharedPtr> &ein)
+{
+    // this piece of code orientates the surface,
+    // it used to be face mesh but its easier to have it here
+    vector<vector<Array<OneD, NekDouble> > > loopt;
+    for (int i = 0; i < ein.size(); i++)
+    {
+        vector<Array<OneD, NekDouble> > loop;
+        for (int j = 0; j < ein[i]->edges.size(); j++)
+        {
+            Array<OneD, NekDouble> bnds = ein[i]->edges[j]->GetBounds();
+            NekDouble dt = (bnds[1] - bnds[0]) / (20 - 1);
+            if (ein[i]->edgeo[j] == eForwards)
+            {
+                for (int k = 0; k < 20; k++)
+                {
+                    NekDouble t = bnds[0] + dt * k;
+                    Array<OneD, NekDouble> l  = ein[i]->edges[j]->P(t);
+                    Array<OneD, NekDouble> uv = surf->locuv(l);
+                    loop.push_back(uv);
+                }
+            }
+            else
+            {
+                for (int k = 19; k >= 0; k--)
+                {
+                    NekDouble t = bnds[0] + dt * k;
+                    Array<OneD, NekDouble> l  = ein[i]->edges[j]->P(t);
+                    Array<OneD, NekDouble> uv = surf->locuv(l);
+                    loop.push_back(uv);
+                }
+            }
+            loopt.push_back(loop);
+        }
+    }
+
+    for (int i = 0; i < loopt.size(); i++)
+    {
+        NekDouble area = 0.0;
+        NekDouble mn   = numeric_limits<double>::max();
+
+        for (int j = 0; j < loopt[i].size(); j++)
+        {
+            mn = min(loopt[i][j][1], mn);
+        }
+
+        for (int j = 0; j < loopt[i].size(); j++)
+        {
+            loopt[i][j][1] -= mn;
+        }
+
+        for (int j = 0; j < loopt[i].size() - 1; j++)
+        {
+            area += (loopt[i][j + 1][0] - loopt[i][j][0]) *
+                    (loopt[i][j][1] + loopt[i][j + 1][1]) / 2.0;
+        }
+        area += (loopt[i][0][0] - loopt[i][loopt[i].size() - 1][0]) *
+                (loopt[i][loopt[i].size() - 1][1] + loopt[i][0][1]) / 2.0;
+
+        ein[i]->area = area;
+    }
+
+    // order by absoulte area
+    int ct = 0;
+    do
+    {
+        ct = 0;
+        for (int i = 0; i < ein.size() - 1; i++)
+        {
+            if (fabs(ein[i]->area) < fabs(ein[i + 1]->area))
+            {
+                // swap
+                swap(ein[i], ein[i + 1]);
+                ct += 1;
+            }
+        }
+
+    } while (ct > 0);
+
+    for (int i = 0; i < ein.size(); i++)
+    {
+        Array<OneD, NekDouble> n1info, n2info;
+        n1info = loopt[i][0];
+        n2info = loopt[i][1];
+
+        Array<OneD, NekDouble> N(2);
+        NekDouble mag = sqrt((n1info[0] - n2info[0]) * (n1info[0] - n2info[0]) +
+                             (n1info[1] - n2info[1]) * (n1info[1] - n2info[1]));
+        ASSERTL0(mag > 1e-30, "infinity");
+        N[0] = -1.0 * (n2info[1] - n1info[1]) / mag;
+        N[1] = (n2info[0] - n1info[0]) / mag;
+
+        Array<OneD, NekDouble> P(2);
+        P[0] = (n1info[0] + n2info[0]) / 2.0 + 1e-8 * N[0];
+        P[1] = (n1info[1] + n2info[1]) / 2.0 + 1e-8 * N[1];
+
+        // now test to see if p is inside or outside the shape
+        // vector to the right
+        int intercepts = 0;
+        for (int j = 0; j < loopt[i].size() - 1; j++)
+        {
+            Array<OneD, NekDouble> nt1, nt2;
+            nt1 = loopt[i][j];
+            nt2 = loopt[i][j + 1];
+
+            if (fabs(nt2[1] - nt1[1]) < 1e-30)
+            {
+                continue;
+            }
+
+            NekDouble lam = (P[1] - nt1[1]) / (nt2[1] - nt1[1]);
+            NekDouble S   = nt1[0] - P[0] + (nt2[0] - nt1[0]) * lam;
+
+            if (!(lam < 0) && !(lam > 1) && S > 0)
+            {
+                intercepts++;
+            }
+        }
+        {
+            Array<OneD, NekDouble> nt1, nt2;
+            nt1 = loopt[i].back();
+            nt2 = loopt[i][0];
+
+            if (fabs(nt2[1] - nt1[1]) < 1e-30)
+            {
+                continue;
+            }
+
+            NekDouble lam = (P[1] - nt1[1]) / (nt2[1] - nt1[1]);
+            NekDouble S   = nt1[0] - P[0] + (nt2[0] - nt1[0]) * lam;
+
+            if (!(lam < 0) && !(lam > 1) && S > 0)
+            {
+                intercepts++;
+            }
+        }
+        if (intercepts % 2 == 0)
+        {
+            P[0]       = (n1info[0] + n2info[0]) / 2.0 - 1e-6 * N[0];
+            P[1]       = (n1info[1] + n2info[1]) / 2.0 - 1e-6 * N[1];
+            intercepts = 0;
+            for (int j = 0; j < loopt[i].size() - 1; j++)
+            {
+                Array<OneD, NekDouble> nt1, nt2;
+                nt1 = loopt[i][j];
+                nt2 = loopt[i][j + 1];
+
+                if (fabs(nt2[1] - nt1[1]) < 1e-30)
+                {
+                    continue;
+                }
+
+                NekDouble lam = (P[1] - nt1[1]) / (nt2[1] - nt1[1]);
+                NekDouble S   = nt1[0] - P[0] + (nt2[0] - nt1[0]) * lam;
+
+                if (!(lam < 0) && !(lam > 1) && S > 0)
+                {
+                    intercepts++;
+                }
+            }
+            {
+                Array<OneD, NekDouble> nt1, nt2;
+                nt1 = loopt[i].back();
+                nt2 = loopt[i][0];
+
+                if (fabs(nt2[1] - nt1[1]) < 1e-30)
+                {
+                    continue;
+                }
+
+                NekDouble lam = (P[1] - nt1[1]) / (nt2[1] - nt1[1]);
+                NekDouble S   = nt1[0] - P[0] + (nt2[0] - nt1[0]) * lam;
+
+                if (!(lam < 0) && !(lam > 1) && S > 0)
+                {
+                    intercepts++;
+                }
+            }
+            if (intercepts % 2 == 0)
+            {
+                cerr << "still failed to find point inside loop" << endl;
+            }
+        }
+
+        ein[i]->center = P;
+    }
+
+    if (ein[0]->area < 0) // reverse the first uvLoop
+    {
+        reverse(ein[0]->edgeo.begin(), ein[0]->edgeo.end());
+        reverse(ein[0]->edges.begin(), ein[0]->edges.end());
+        // need to flip edgeo
+        for (int i = 0; i < ein[0]->edgeo.size(); i++)
+        {
+            if (ein[0]->edgeo[i] == eForwards)
+            {
+                ein[0]->edgeo[i] = eBackwards;
+            }
+            else
+            {
+                ein[0]->edgeo[i] = 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());
+
+            // need to flip edgeo
+            for (int j = 0; j < ein[i]->edgeo.size(); j++)
+            {
+                if (ein[i]->edgeo[j] == eForwards)
+                {
+                    ein[i]->edgeo[j] = eBackwards;
+                }
+                else
+                {
+                    ein[i]->edgeo[j] = eForwards;
+                }
+            }
+        }
+    }
+}
+
+}
+}

library/NekMeshUtils/CADSystem/CADSurf.h

@@ -37,8 +37,8 @@
 #define NekMeshUtils_CADSYSTEM_CADSURF
 
 #include <NekMeshUtils/CADSystem/CADObject.h>
-#include <NekMeshUtils/CADSystem/CADVert.h>
 #include <NekMeshUtils/CADSystem/CADSystem.h>
+#include <NekMeshUtils/CADSystem/CADVert.h>
 
 namespace Nektar
 {
@@ -69,6 +69,9 @@ public:
     {
     }
 
+    static void OrientateEdges(
+        CADSurfSharedPtr surf, std::vector<EdgeLoopSharedPtr> &ein);
+
     /**
      * @brief Get the loop structures which bound the cad surface
      */
@@ -77,6 +80,11 @@ public:
         return m_edges;
     }
 
+    void SetEdges(std::vector<EdgeLoopSharedPtr> ein)
+    {
+        m_edges = ein;
+    }
+
     /**
      * @brief Get the limits of the parametric space for the surface.
      *
@@ -90,7 +98,7 @@ public:
      * @param uv Array of u and v parametric coords.
      * @return Array of xyz components of normal vector.
      */
-    virtual Array<OneD, NekDouble> N    (Array<OneD, NekDouble> uv) = 0;
+    virtual Array<OneD, NekDouble> N(Array<OneD, NekDouble> uv) = 0;
 
     /**
      * @brief Get the set of first derivatives at parametric point u,v
@@ -98,7 +106,7 @@ public:
      * @param uv Array of u and v parametric coords.
      * @return Array of xyz copmonents of first derivatives.
      */
-    virtual Array<OneD, NekDouble> D1   (Array<OneD, NekDouble> uv) = 0;
+    virtual Array<OneD, NekDouble> D1(Array<OneD, NekDouble> uv) = 0;
 
     /**
      * @brief Get the set of second derivatives at parametric point u,v
@@ -106,7 +114,7 @@ public:
      * @param uv array of u and v parametric coords
      * @return array of xyz copmonents of second derivatives
      */
-    virtual Array<OneD, NekDouble> D2   (Array<OneD, NekDouble> uv) = 0;
+    virtual Array<OneD, NekDouble> D2(Array<OneD, NekDouble> uv) = 0;
 
     /**
      * @brief Get the x,y,z at parametric point u,v.
@@ -114,7 +122,7 @@ public:
      * @param uv Array of u and v parametric coords.
      * @return Array of xyz location.
      */
-    virtual Array<OneD, NekDouble> P    (Array<OneD, NekDouble> uv) = 0;
+    virtual Array<OneD, NekDouble> P(Array<OneD, NekDouble> uv) = 0;
 
     /**
      * @brief Performs a reverse look up to find u,v and x,y,z.
@@ -126,7 +134,8 @@ public:
 
     /**
      * @brief does unconstrained locuv to project point from anywhere
-     * and calculate the distance between the orthonormal projection to the surface
+     * and calculate the distance between the orthonormal projection to the
+     * surface
      * and the point
      */
     virtual NekDouble DistanceTo(Array<OneD, NekDouble> p) = 0;
@@ -135,7 +144,8 @@ public:
      * @brief takes a point from anywhere find the nearest surface point and its
      * uv
      */
-    virtual void ProjectTo(Array<OneD, NekDouble> &tp, Array<OneD, NekDouble> &uv) = 0;
+    virtual void ProjectTo(Array<OneD, NekDouble> &tp,
+                           Array<OneD, NekDouble> &uv) = 0;
 
     /**
      * @brief returns curvature at point uv
@@ -151,7 +161,6 @@ public:
     }
 
 protected:
-
     /**
      * @brief sets the flag to reverse the normal for this suface,
      * this is determined in cadsystem and ensures all surface normals,
@@ -161,7 +170,7 @@ protected:
     {
         m_correctNormal = false;
     }
-    
+
     /// normal
     bool m_correctNormal;
     /// List of bounding edges in loops with orientation.
@@ -175,8 +184,7 @@ typedef boost::shared_ptr<CADSurf> CADSurfSharedPtr;
 
 typedef LibUtilities::NekFactory<std::string, CADSurf> CADSurfFactory;
 
-CADSurfFactory& GetCADSurfFactory();
-
+CADSurfFactory &GetCADSurfFactory();
 }
 }
 

library/NekMeshUtils/CADSystem/CADSystem.h

@@ -58,6 +58,12 @@ typedef boost::shared_ptr<CADCurve> CADCurveSharedPtr;
 class CADSurf;
 typedef boost::shared_ptr<CADSurf> CADSurfSharedPtr;
 
+enum Orientation
+{
+    eForwards,
+    eBackwards
+};
+
 /**
  * @brief struct which descibes a collection of cad edges which are a
  *        loop on the cad surface
@@ -65,7 +71,7 @@ typedef boost::shared_ptr<CADSurf> CADSurfSharedPtr;
 struct EdgeLoop
 {
     std::vector<CADCurveSharedPtr> edges;
-    std::vector<int> edgeo; //0 is forward 1 is backward
+    std::vector<Orientation> edgeo;
     Array<OneD, NekDouble> center;
     NekDouble area;
 };

library/NekMeshUtils/CADSystem/OCE/CADCurveOCE.cpp

@@ -70,7 +70,7 @@ NekDouble CADCurveOCE::tAtArcLength(NekDouble s)
 
 NekDouble CADCurveOCE::Length(NekDouble ti, NekDouble tf)
 {
-    Array<OneD, NekDouble> b = Bounds();
+    Array<OneD, NekDouble> b = GetBounds();
     Handle(Geom_Curve) NewCurve = new Geom_TrimmedCurve(m_c, ti, tf);
     TopoDS_Edge NewEdge = BRepBuilderAPI_MakeEdge(NewCurve);
     GProp_GProps System;
@@ -81,7 +81,7 @@ NekDouble CADCurveOCE::Length(NekDouble ti, NekDouble tf)
 NekDouble CADCurveOCE::loct(Array<OneD, NekDouble> xyz)
 {
     NekDouble t = 0.0;
-    Array<OneD, NekDouble> b = Bounds();
+    Array<OneD, NekDouble> b = GetBounds();
 
     gp_Pnt loc(xyz[0]*1000.0, xyz[1]*1000.0, xyz[2]*1000.0);
 
@@ -136,7 +136,7 @@ NekDouble CADCurveOCE::Curvature(NekDouble t)
     return d.Curvature() * 1000.0;
 }
 
-Array<OneD, NekDouble> CADCurveOCE::Bounds()
+Array<OneD, NekDouble> CADCurveOCE::GetBounds()
 {
     Array<OneD, NekDouble> t(2);
     t[0] = m_occCurve.FirstParameter();

library/NekMeshUtils/CADSystem/OCE/CADCurveOCE.h

@@ -63,7 +63,7 @@ public:
     {
     }
 
-    virtual Array<OneD, NekDouble> Bounds();
+    virtual Array<OneD, NekDouble> GetBounds();
     virtual NekDouble Length(NekDouble ti, NekDouble tf);
     virtual Array<OneD, NekDouble> P(NekDouble t);
     virtual Array<OneD, NekDouble> D2(NekDouble t);
@@ -88,7 +88,7 @@ public:
         BRepGProp::LinearProperties(m_occEdge, System);
         m_length = System.Mass();
 
-        Array<OneD, NekDouble> b = Bounds();
+        Array<OneD, NekDouble> b = GetBounds();
         m_c = BRep_Tool::Curve(TopoDS::Edge(cp), b[0], b[1]);
 
         m_id   = i;

library/NekMeshUtils/CADSystem/OCE/CADSurfOCE.cpp

@@ -45,8 +45,7 @@ namespace NekMeshUtils
 std::string CADSurfOCE::key = GetCADSurfFactory().RegisterCreatorFunction(
     "oce", CADSurfOCE::create, "CADSurfOCE");
 
-void CADSurfOCE::Initialise(int i, TopoDS_Shape in,
-                            vector<EdgeLoopSharedPtr> ein)
+void CADSurfOCE::Initialise(int i, TopoDS_Shape in)
 {
     // this bit of code changes the units of the cad from mm opencascade
     // defualt to m
@@ -63,7 +62,6 @@ void CADSurfOCE::Initialise(int i, TopoDS_Shape in,
         SetReverseNomral();
     }
 
-    m_edges = ein;
     gp_Trsf transform;
     gp_Pnt ori(0.0, 0.0, 0.0);
     transform.SetScale(ori, 1.0 / 1000.0);

library/NekMeshUtils/CADSystem/OCE/CADSurfOCE.h

@@ -63,7 +63,7 @@ public:
     {
     }
 
-    void Initialise(int i, TopoDS_Shape in, std::vector<EdgeLoopSharedPtr> ein);
+    void Initialise(int i, TopoDS_Shape in);
 
     virtual Array<OneD, NekDouble> GetBounds();
     virtual Array<OneD, NekDouble> N    (Array<OneD, NekDouble> uv);

library/NekMeshUtils/CADSystem/OCE/CADSystemOCE.cpp

@@ -34,11 +34,11 @@
 ////////////////////////////////////////////////////////////////////////////////
 #include <LibUtilities/BasicUtils/ParseUtils.hpp>
 
-#include <NekMeshUtils/CADSystem/OCE/CADSystemOCE.h>
-#include <NekMeshUtils/CADSystem/OCE/CADVertOCE.h>
+#include <NekMeshUtils/CADSystem/CADSurf.h>
 #include <NekMeshUtils/CADSystem/OCE/CADCurveOCE.h>
 #include <NekMeshUtils/CADSystem/OCE/CADSurfOCE.h>
-#include <NekMeshUtils/CADSystem/CADSurf.h>
+#include <NekMeshUtils/CADSystem/OCE/CADSystemOCE.h>
+#include <NekMeshUtils/CADSystem/OCE/CADVertOCE.h>
 
 using namespace std;
 
@@ -54,7 +54,7 @@ bool CADSystemOCE::LoadCAD()
 {
     if (m_naca.size() == 0)
     {
-        //not a naca profile behave normally
+        // not a naca profile behave normally
         // Takes step file and makes OpenCascade shape
         STEPControl_Reader reader;
         reader = STEPControl_Reader();
@@ -112,9 +112,6 @@ bool CADSystemOCE::LoadCAD()
         }
     }
 
-    // from id of curve to list of ids of surfs
-    map<int, vector<int> > adjsurfmap;
-
     // Adds edges to our type and map
     for (int i = 1; i <= mapOfEdges.Extent(); i++)
     {
@@ -146,79 +143,15 @@ bool CADSystemOCE::LoadCAD()
         TopTools_IndexedMapOfShape mapOfWires;
         TopExp::MapShapes(face, TopAbs_WIRE, mapOfWires);
 
-        // this pice of code does an idiot check on the loops to make sure
-        // they dont cross or touch
-        if (mapOfWires.Extent() > 1)
-        {
-            TopoDS_Wire ow = BRepTools::OuterWire(TopoDS::Face(face));
-
-            vector<TopoDS_Shape> wirefacecuts;
-            vector<gp_Pnt> centersofcutfaces;
-
-            for (int j = 1; j <= mapOfWires.Extent(); j++)
-            {
-                TopoDS_Shape wire = mapOfWires.FindKey(j);
-
-                if (wire != ow)
-                {
-                    BRepBuilderAPI_MakeFace build(
-                        BRep_Tool::Surface(TopoDS::Face(face)), 1e-7);
-                    build.Add(TopoDS::Wire(wire));
-                    TopoDS_Shape newface = build.Shape();
-                    wirefacecuts.push_back(newface);
-                    BRepAdaptor_Surface b =
-                        BRepAdaptor_Surface(TopoDS::Face(newface));
-                    NekDouble u, v;
-                    u = (b.LastUParameter() - b.FirstUParameter()) / 2.0;
-                    v = (b.LastVParameter() - b.FirstVParameter()) / 2.0;
-                    centersofcutfaces.push_back(b.Value(u, v));
-                }
-            }
-            for (int j = 0; j < wirefacecuts.size(); j++)
-            {
-                for (int k = 0; k < wirefacecuts.size(); k++)
-                {
-                    if (j == k)
-                        continue;
-
-                    /// TODO fix this test
-                    BRepClass_FaceClassifier fc(TopoDS::Face(wirefacecuts[j]),
-                                                centersofcutfaces[k], 1e-7);
-                    // ASSERTL0(fc.State() == 1, "Internal face loops make this
-                    // cad impossible to mesh");