Merge branch 'fix/3DbndLayer' into 'master'

fix/3d-boundary-layers

Recent merge of feature/mesh_spacing broke a number of parts of the mesh generation, mostly do with orientation.

The CAD orientation routines have now been rewritten in this branch and the relevant code updated.

3 more regression tests, which didnt work prior to this fix, have been added

A small amount of wall time optimisation has been done to the OCC cad interface and the reduction of some L0 asserts to L1.

See merge request !744

library/NekMeshUtils/2DGenerator/2DGenerator.cpp

@@ -52,7 +52,7 @@ ModuleKey Generator2D::className = GetModuleFactory().RegisterCreatorFunction(
 Generator2D::Generator2D(MeshSharedPtr m) : ProcessModule(m)
 {
     m_config["blcurves"] =
-        ConfigOption(false, "0", "Generate parallelograms on these curves");
+        ConfigOption(false, "", "Generate parallelograms on these curves");
     m_config["blthick"] =
         ConfigOption(false, "0.0", "Parallelogram layer thickness");
 }
@@ -74,6 +74,9 @@ void Generator2D::Process()
     m_thickness_ID =
         m_thickness.DefineFunction("x y z", m_config["blthick"].as<string>());
 
+    ParseUtils::GenerateSeqVector(m_config["blcurves"].as<string>().c_str(),
+                                  m_blCurves);
+
     // linear mesh all curves
     for (int i = 1; i <= m_mesh->m_cad->GetNumCurve(); i++)
     {
@@ -102,63 +105,24 @@ void Generator2D::Process()
 
     ////////////////////////////////////////
 
-    EdgeSet::iterator it;
-    for (it = m_mesh->m_edgeSet.begin(); it != m_mesh->m_edgeSet.end(); it++)
-    {
-        vector<NodeSharedPtr> ns;
-        ns.push_back((*it)->m_n1);
-        ns.push_back((*it)->m_n2);
-
-        // for each iterator create a LibUtilities::eSegement
-        // push segment into m_mesh->m_element[1]
-        // tag for the elements shoudl be the CAD number of the curves
-
-        ElmtConfig conf(LibUtilities::eSegment, 1, false, false);
-
-        vector<int> tags;
-        tags.push_back((*it)->m_parentCAD->GetId());
-
-        ElementSharedPtr E2 = GetElementFactory().CreateInstance(
-            LibUtilities::eSegment, conf, ns, tags);
-
-        m_mesh->m_element[1].push_back(E2);
-    }
-
-    for (int i = 1; i <= m_mesh->m_cad->GetNumSurf(); i++)
-    {
-        m_facemeshes[i] = MemoryManager<FaceMesh>::AllocateSharedPtr(
-            i, m_mesh, m_curvemeshes, 100);
-
-        m_facemeshes[i]->OrientateCurves();
-    }
-
     if (m_config["blcurves"].beenSet)
     {
         // we need to do the boundary layer generation in a face by face basis
         MakeBLPrep();
 
-        // Im going to do a horrendous trick to get the edge orientaion.
-        // Going to activate the first routine of facemeshing without actually
-        // face meshing, this will orientate the edgeloop objects (hopefully);
-        // which can be used by the makebl command to know the normal
-        // orienation
         for (int i = 1; i <= m_mesh->m_cad->GetNumSurf(); i++)
         {
-            MakeBL(i, m_facemeshes[i]->GetEdges());
+            MakeBL(i);
         }
     }
 
-    // m_mesh->m_element[1].clear();
-
     if (m_mesh->m_verbose)
     {
         cout << endl << "\tFace meshing:" << endl << endl;
     }
 
     // linear mesh all surfaces
-    map<int, FaceMeshSharedPtr>::iterator fit;
-    int i = 1;
-    for (fit = m_facemeshes.begin(); fit != m_facemeshes.end(); fit++)
+    for (int i = 1; i <= m_mesh->m_cad->GetNumSurf(); i++)
     {
         if (m_mesh->m_verbose)
         {
@@ -166,16 +130,34 @@ void Generator2D::Process()
                                            "Face progress");
         }
 
-        if (m_config["blcurves"].beenSet)
-        {
-            // for bl surfaces orientate curves needs to be run again to
-            // push nodes to the edges of the system
-            fit->second->ResetCurvemeshes(m_curvemeshes);
-            fit->second->OrientateCurves();
-        }
+        m_facemeshes[i] =
+            MemoryManager<FaceMesh>::AllocateSharedPtr(i,m_mesh,
+                m_curvemeshes, 99+i);
+        m_facemeshes[i]->Mesh();
+    }
+
+    ////////////////////////////////////
 
-        fit->second->Mesh();
-        i++;
+    EdgeSet::iterator it;
+    for (it = m_mesh->m_edgeSet.begin(); it != m_mesh->m_edgeSet.end(); it++)
+    {
+        vector<NodeSharedPtr> ns;
+        ns.push_back((*it)->m_n1);
+        ns.push_back((*it)->m_n2);
+
+        // for each iterator create a LibUtilities::eSegement
+        // push segment into m_mesh->m_element[1]
+        // tag for the elements shoudl be the CAD number of the curves
+
+        ElmtConfig conf(LibUtilities::eSegment, 1, false, false);
+
+        vector<int> tags;
+        tags.push_back((*it)->m_parentCAD->GetId());
+
+        ElementSharedPtr E2 = GetElementFactory().CreateInstance(
+            LibUtilities::eSegment, conf, ns, tags);
+
+        m_mesh->m_element[1].push_back(E2);
     }
 
     ProcessVertices();
@@ -195,8 +177,6 @@ void Generator2D::MakeBLPrep()
     }
 
     // identify the nodes which will become the boundary layer.
-    ParseUtils::GenerateSeqVector(m_config["blcurves"].as<string>().c_str(),
-                                  m_blCurves);
 
     for (vector<unsigned>::iterator it = m_blCurves.begin();
          it != m_blCurves.end(); ++it)
@@ -210,17 +190,8 @@ void Generator2D::MakeBLPrep()
     }
 }
 
-void Generator2D::MakeBL(int faceid, vector<EdgeLoopSharedPtr> e)
+void Generator2D::MakeBL(int faceid)
 {
-    map<int, int> edgeToOrient;
-    for (int j = 0; j < e.size(); j++)
-    {
-        for (int i = 0; i < e[j]->edges.size(); ++i)
-        {
-            edgeToOrient[e[j]->edges[i]->GetId()] = e[j]->edgeo[i];
-        }
-    }
-
     map<int, Array<OneD, NekDouble> > edgeNormals;
 
     int eid = 0;
@@ -228,25 +199,45 @@ void Generator2D::MakeBL(int faceid, vector<EdgeLoopSharedPtr> e)
     for (vector<unsigned>::iterator it = m_blCurves.begin();
          it != m_blCurves.end(); ++it)
     {
-        int edgeo = edgeToOrient[*it];
+        CADOrientation::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 =
-                (edgeo == 0) ? es[j]->m_n1->GetLoc() : es[j]->m_n2->GetLoc();
-            Array<OneD, NekDouble> p2 =
-                (edgeo == 0) ? es[j]->m_n2->GetLoc() : es[j]->m_n1->GetLoc();
+            Array<OneD, NekDouble> p1, p2;
+            p1 = es[j]->m_n1->GetCADSurfInfo(faceid);
+            p2 = es[j]->m_n2->GetCADSurfInfo(faceid);
+            if (edgeo == CADOrientation::eBackwards)
+            {
+                swap(p1, p2);
+            }
             Array<OneD, NekDouble> n(2);
             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;
         }
     }
@@ -285,7 +276,8 @@ void Generator2D::MakeBL(int faceid, vector<EdgeLoopSharedPtr> e)
     for (vector<unsigned>::iterator it = m_blCurves.begin();
          it != m_blCurves.end(); ++it)
     {
-        int edgeo = edgeToOrient[*it];
+        CADOrientation::Orientation edgeo =
+            m_mesh->m_cad->GetCurve(*it)->GetOrienationWRT(faceid);
 
         vector<NodeSharedPtr> ns = m_curvemeshes[*it]->GetMeshPoints();
         vector<NodeSharedPtr> newNs;
@@ -296,7 +288,7 @@ void Generator2D::MakeBL(int faceid, vector<EdgeLoopSharedPtr> e)
         m_curvemeshes[*it] =
             MemoryManager<CurveMesh>::AllocateSharedPtr(*it, m_mesh, newNs);
 
-        if (edgeo == 1)
+        if (edgeo == CADOrientation::eBackwards)
         {
             reverse(ns.begin(), ns.end());
         }

library/NekMeshUtils/2DGenerator/2DGenerator.h

@@ -69,7 +69,7 @@ private:
 
     void MakeBLPrep();
 
-    void MakeBL(int faceid, std::vector<EdgeLoopSharedPtr> e);
+    void MakeBL(int faceid);
 
     void Report();
     /// map of individual surface meshes from parametric surfaces

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, CADOrientation::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, CADOrientation::Orientation> > GetAdjSurf()
     {
         return m_adjSurfs;
     }
@@ -163,12 +163,30 @@ public:
      */
     virtual NekDouble loct(Array<OneD, NekDouble> xyz) = 0;
 
+    CADOrientation::Orientation GetOrienationWRT(int surf)
+    {
+        for(int i = 0; i < m_adjSurfs.size(); i++)
+        {
+            if(m_adjSurfs[i].first->GetId() == surf)
+            {
+                return m_adjSurfs[i].second;
+            }
+        }
+
+        ASSERTL0(false,"surf not in adjecency list");
+        return CADOrientation::eUnknown;
+    }
+
+    virtual Array<OneD, NekDouble> NormalWRT(NekDouble t, int surf)=0;
+    virtual Array<OneD, NekDouble> N(NekDouble t)=0;
+
+
 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, CADOrientation::Orientation> > m_adjSurfs;
     /// list of end vertices
     std::vector<CADVertSharedPtr> m_mainVerts;
 };

library/NekMeshUtils/CADSystem/CADObject.h

@@ -55,6 +55,16 @@ enum cadType
 };
 }
 
+namespace CADOrientation
+{
+enum Orientation
+{
+    eUnknown,
+    eForwards,
+    eBackwards
+};
+}
+
 class CADObject
 {
 public:
@@ -67,7 +77,9 @@ public:
     {
     }
 
-    virtual ~CADObject(){}
+    virtual ~CADObject()
+    {
+    }
 
     /**
      * @brief Return ID of the vertex
@@ -82,11 +94,19 @@ public:
         return m_type;
     }
 
+    virtual CADOrientation::Orientation Orientation()
+    {
+        ASSERTL0(false,"must be implemented at the cad object level");
+        return CADOrientation::eUnknown;
+    }
+
 protected:
     /// ID of the vert.
     int m_id;
     /// type of the cad object
     CADType::cadType m_type;
+    /// orientation of the CADObject
+    CADOrientation::Orientation m_orientation;
 };
 
 typedef boost::shared_ptr<CADObject> CADObjectSharedPtr;

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"
+
+#include <boost/geometry.hpp>
+#include <boost/geometry/algorithms/assign.hpp>
+#include <boost/geometry/geometries/point_xy.hpp>
+#include <boost/geometry/geometries/polygon.hpp>
+
+namespace bg = boost::geometry;
+typedef bg::model::d2::point_xy<double> point_xy;
+
+using namespace std;
+
+namespace Nektar
+{
+namespace NekMeshUtils
+{
+
+void CADSurf::OrientateEdges(CADSurfSharedPtr surf,
+                             vector<CADSystem::EdgeLoopSharedPtr> &ein)
+{
+    // this piece of code orientates the surface,
+    // it used to be face mesh but its easier to have it here
+    int np = 20;
+    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]) / (np - 1);
+            if (ein[i]->edgeo[j] == CADOrientation::eForwards)
+            {
+                for (int k = 0; k < np - 1; 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 = np - 1; 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);
+    }
+
+    vector<bg::model::polygon<point_xy, false, true> > polygons;
+
+    for (int i = 0; i < loopt.size(); i++)
+    {
+        bg::model::polygon<point_xy, false, true> polygon;
+        vector<point_xy> points;
+        for (int j = 0; j < loopt[i].size(); j++)
+        {
+            points.push_back(point_xy(loopt[i][j][0], loopt[i][j][1]));
+        }
+        // boost requires for closed polygons (last point == first point)
+        points.push_back(point_xy(loopt[i][0][0], loopt[i][0][1]));
+
+        bg::assign_points(polygon, points);
+
+        NekDouble area = bg::area(polygon);
+
+        ein[i]->area = area;
+
+        point_xy cen;
+        bg::centroid(polygon, cen);
+
+        ein[i]->center    = Array<OneD, NekDouble>(2);
+        ein[i]->center[0] = cen.x();
+        ein[i]->center[1] = cen.y();
+
+        polygons.push_back(polygon);
+    }
+
+    // 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]);
+                swap(loopt[i], loopt[i + 1]);
+                swap(polygons[i], polygons[i + 1]);
+                ct += 1;
+            }
+        }
+
+    } while (ct > 0);
+
+    // only need center points for inner loops
+    for (int i = 1; i < ein.size(); i++)
+    {
+        point_xy p(ein[i]->center[0], ein[i]->center[1]);
+
+        if (!bg::within(p, polygons[i]))
+        {
+            Array<OneD, NekDouble> n1 = loopt[i][0];
+            Array<OneD, NekDouble> n2 = loopt[i][1];
+
+            Array<OneD, NekDouble> N(2);
+            NekDouble mag = sqrt((n1[0] - n2[0]) * (n1[0] - n2[0]) +
+                                 (n1[1] - n2[1]) * (n1[1] - n2[1]));
+
+            N[0] = (n2[1] - n1[1]) / mag;
+            N[1] = -1.0 * (n2[0] - n1[0]) / mag;
+
+            Array<OneD, NekDouble> P(2);
+            P[0] = n1[0] + N[0];
+            P[1] = n1[1] + N[1];
+
+            ein[i]->center = P;
+
+            p = point_xy(P[0], P[1]);
+
+            ASSERTL0(boost::geometry::within(p, polygons[i]), "point is not side loop");
+        }
+    }
+}
+}
+}

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
 {
@@ -63,20 +63,29 @@ public:
     CADSurf()
     {
         m_type = CADType::eSurf;
+        m_orientation = CADOrientation::eForwards;
     }
 
     ~CADSurf()
     {
     }
 
+    static void OrientateEdges(
+        CADSurfSharedPtr surf, std::vector<CADSystem::EdgeLoopSharedPtr> &ein);
+
     /**
      * @brief Get the loop structures which bound the cad surface
      */
-    std::vector<EdgeLoopSharedPtr> GetEdges()
+    std::vector<CADSystem::EdgeLoopSharedPtr> GetEdges()
     {
         return m_edges;
     }
 
+    void SetEdges(std::vector<CADSystem::EdgeLoopSharedPtr> ein)
+    {
+        m_edges = ein;
+    }
+
     /**
      * @brief Get the limits of the parametric space for the surface.
      *
@@ -90,7 +99,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 +107,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 +115,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 +123,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 +135,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 +145,8 @@ public:
      * @brief takes a point from anywhere find the nearest surface point and its