more updates

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++)
     {
@@ -98,6 +101,7 @@ void Generator2D::Process()
         }
 
         m_curvemeshes[i]->Mesh();
+
     }
 
     ////////////////////////////////////////
@@ -107,14 +111,9 @@ void Generator2D::Process()
         // 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);
+            //MakeBL(i);
         }
     }
 
@@ -179,8 +178,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)

library/NekMeshUtils/CADSystem/OCE/CADCurveOCE.cpp

@@ -138,8 +138,8 @@ Array<OneD, NekDouble> CADCurveOCE::NormalWRT(NekDouble t, int surf)
     Array<OneD, NekDouble> uv = surface.first->locuv(p);
     Array<OneD, NekDouble> d1 = surface.first->D1(uv);
 
-    NekDouble t1 = t - 1e-6;
-    NekDouble t2 = t + 1e-6;
+    NekDouble t1 = t - 1e-8;
+    NekDouble t2 = t + 1e-8;
 
     if(surface.second == CADSystem::eBackwards)
     {
@@ -149,8 +149,8 @@ Array<OneD, NekDouble> CADCurveOCE::NormalWRT(NekDouble t, int surf)
     Array<OneD, NekDouble> uv1 = surface.first->locuv(P(t1));
     Array<OneD, NekDouble> uv2 = surface.first->locuv(P(t2));
 
-    NekDouble du = uv2[0] - uv1[0];
-    NekDouble dv = uv2[1] - uv1[1];
+    NekDouble du = uv2[1] - uv1[1];
+    NekDouble dv = -1.0*(uv2[0] - uv1[0]);
 
     Array<OneD, NekDouble> N(3,0.0);
     N[0] = (d1[3] * du + d1[6] * dv) / 2.0;
@@ -168,7 +168,7 @@ Array<OneD, NekDouble> CADCurveOCE::NormalWRT(NekDouble t, int surf)
 Array<OneD, NekDouble> CADCurveOCE::N(NekDouble t)
 {
     GeomLProp_CLProps d(m_c,2,1e-8);
-    d.SetParameter(t);
+    d.SetParameter(t+1e-8);
 
     gp_Vec d2 = d.D2();
     if(d2.Magnitude() < 1e-8)

library/NekMeshUtils/SurfaceMeshing/CurveMesh.cpp

@@ -49,7 +49,7 @@ void CurveMesh::Mesh()
     m_bounds      = m_cadcurve->GetBounds();
     m_curvelength = m_cadcurve->GetTotLength();
     m_numSamplePoints =
-        int(m_curvelength / m_mesh->m_octree->GetMinDelta()) + 5;
+        int(m_curvelength / m_mesh->m_octree->GetMinDelta()) * 2;
     ds = m_curvelength / (m_numSamplePoints - 1);
 
     GetSampleFunction();
@@ -214,6 +214,8 @@ NekDouble CurveMesh::EvaluateDS(NekDouble s)
     int a = 0;
     int b = 0;
 
+    ASSERTL1(!(s < 0) && !(s > m_curvelength),"s out of bounds");
+
     if (s == 0)
     {
         return m_dst[0][0];
@@ -251,6 +253,8 @@ NekDouble CurveMesh::EvaluatePS(NekDouble s)
     int a = 0;
     int b = 0;
 
+    ASSERTL1(!(s < 0) && !(s > m_curvelength),"s out of bounds");
+
     if (s == 0)
     {
         return m_ps[0][0];
@@ -294,63 +298,55 @@ NekDouble CurveMesh::EvaluatePS(NekDouble s)
 void CurveMesh::GetSampleFunction()
 {
     m_dst.resize(m_numSamplePoints);
-    Array<OneD, NekDouble> loc = m_cadcurve->P(m_bounds[0]);
 
     vector<NekDouble> dsti;
     dsti.resize(3);
 
-    dsti[0] = m_mesh->m_octree->Query(loc);
-    dsti[1] = 0.0;
-    dsti[2] = m_bounds[0];
-
-    m_dst[0] = dsti;
-
-    for (int i = 1; i < m_numSamplePoints; i++)
+    for (int i = 0; i < m_numSamplePoints; i++)
     {
         dsti[1]     = i * ds;
         NekDouble t = m_cadcurve->tAtArcLength(dsti[1]);
 
-        loc = m_cadcurve->P(t);
+        Array<OneD, NekDouble> loc = m_cadcurve->P(t);
 
-        NekDouble ts =
+        /*NekDouble ts =
             m_bl.Evaluate(m_blID, loc[0], loc[1], loc[2], 0.0);
 
         if (ts > 0.0)
         {
             Array<OneD, NekDouble> N = m_cadcurve->N(t);
+            Array<OneD, NekDouble> Nwrt = m_cadcurve->NormalWRT(t, 0);
+
             if(N[0]*N[0] + N[1]*N[1] + N[2]*N[2] < 1e-6)
             {
                 dsti[0] = m_mesh->m_octree->Query(loc);
             }
-
-            Array<OneD, NekDouble> Nwrt = m_cadcurve->NormalWRT(t, 0);
-            NekDouble R = 1.0 / m_cadcurve->Curvature(t);
-            NekDouble dot = N[0]*Nwrt[0] + N[1]*Nwrt[1] + N[2]*Nwrt[2];
-            bool concave = false;
-            if(dot > 0)
-            {
-                concave = true;
-            }
-            Array<OneD, NekDouble> tloc(3);
-            tloc[0] = loc[0] + ts * Nwrt[0];
-            tloc[1] = loc[1] + ts * Nwrt[1];
-            tloc[2] = loc[2] + ts * Nwrt[2];
-
-            NekDouble d = m_mesh->m_octree->Query(tloc);
-
-            if(concave)
+            else if ( N[0]*Nwrt[0] + N[1]*Nwrt[1] + N[2]*Nwrt[2] > 0)
             {
-                dsti[0] = d * R / (R - ts);
+                //concave
+                dsti[0] = m_mesh->m_octree->Query(loc);
             }
             else
             {
+                NekDouble R = 1.0 / m_cadcurve->Curvature(t);
+                if(R > 2.0*t)
+                {
+                    R = 2.0*t;
+                }
+                Array<OneD, NekDouble> tloc(3);
+                tloc[0] = loc[0] + ts * Nwrt[0];
+                tloc[1] = loc[1] + ts * Nwrt[1];
+                tloc[2] = loc[2] + ts * Nwrt[2];
+
+                NekDouble d = m_mesh->m_octree->Query(tloc);
+
                 dsti[0] = d * R / (R + ts);
             }
         }
         else
-        {
+        {*/
             dsti[0] = m_mesh->m_octree->Query(loc);
-        }
+        //}
 
         dsti[2] = t;