Merge branch 'feature/2D_periodic_BCs' into 'master'

2D Periodic Boundary Conditions

Adds the ability to  make pairs of curves periodic

See merge request !733

CHANGELOG.md

@@ -4,6 +4,8 @@ Changelog
 v4.5.0
 ------
 **NekMesh**:
+
+- Add periodic boundary condition meshing in 2D (!733)
 - Adjust boundary layer thickness in corners in 2D (!739)
 
 **Documentation**:
@@ -108,9 +110,9 @@ v4.4.0
 - Change variable names in mcf file to make more sense (!736)
 - Fix issues in varopti module so that in can be compiled without meshgen on
   (!736)
-- Replace LAPACK Eigenvalue calculation with handwritten function in 
+- Replace LAPACK Eigenvalue calculation with handwritten function in
   varopti (!738)
-- Improved node-colouring algorithm for better load-balancing 
+- Improved node-colouring algorithm for better load-balancing
   in varopti (!738)
 - Simplified calculation of the energy functional in varopti for improved
   performance (!738)

library/NekMeshUtils/2DGenerator/2DGenerator.cpp

@@ -40,6 +40,8 @@
 #include <LibUtilities/BasicUtils/ParseUtils.hpp>
 #include <LibUtilities/BasicUtils/Progressbar.hpp>
 
+#include <boost/algorithm/string.hpp>
+
 using namespace std;
 namespace Nektar
 {
@@ -56,6 +58,8 @@ Generator2D::Generator2D(MeshSharedPtr m) : ProcessModule(m)
         ConfigOption(false, "", "Generate parallelograms on these curves");
     m_config["blthick"] =
         ConfigOption(false, "0.0", "Parallelogram layer thickness");
+    m_config["periodic"] =
+        ConfigOption(false, "", "Set of pairs of periodic curves");
     m_config["bltadjust"] =
         ConfigOption(false, "2.0", "Boundary layer thickness adjustment");
     m_config["adjustblteverywhere"] =
@@ -73,12 +77,9 @@ void Generator2D::Process()
         cout << endl << "2D meshing" << endl;
         cout << endl << "\tCurve meshing:" << endl << endl;
     }
-
     m_mesh->m_numNodes = m_mesh->m_cad->GetNumVerts();
-
     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);
 
@@ -90,10 +91,8 @@ void Generator2D::Process()
             LibUtilities::PrintProgressbar(i, m_mesh->m_cad->GetNumCurve(),
                                            "Curve progress");
         }
-
         vector<unsigned int>::iterator f =
             find(m_blCurves.begin(), m_blCurves.end(), i);
-
         if (f == m_blCurves.end())
         {
             m_curvemeshes[i] =
@@ -104,17 +103,24 @@ void Generator2D::Process()
             m_curvemeshes[i] = MemoryManager<CurveMesh>::AllocateSharedPtr(
                 i, m_mesh, m_config["blthick"].as<string>());
         }
-
         m_curvemeshes[i]->Mesh();
     }
 
+    ////////
+    // consider periodic curves
+
+    if (m_config["periodic"].beenSet)
+    {
+        PeriodicPrep();
+        MakePeriodic();
+    }
+
     ////////////////////////////////////////
 
     if (m_config["blcurves"].beenSet)
     {
         // we need to do the boundary layer generation in a face by face basis
         MakeBLPrep();
-
         for (int i = 1; i <= m_mesh->m_cad->GetNumSurf(); i++)
         {
             MakeBL(i);
@@ -125,7 +131,6 @@ void Generator2D::Process()
     {
         cout << endl << "\tFace meshing:" << endl << endl;
     }
-
     // linear mesh all surfaces
     for (int i = 1; i <= m_mesh->m_cad->GetNumSurf(); i++)
     {
@@ -148,28 +153,25 @@ void Generator2D::Process()
         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);
     }
 
+    // m_mesh->m_expDim = 1;
+    // m_mesh->m_element[2].clear();
+
     ProcessVertices();
     ProcessEdges();
     ProcessFaces();
     ProcessElements();
     ProcessComposites();
-
     Report();
 }
 
@@ -197,17 +199,13 @@ void Generator2D::MakeBLPrep()
 void Generator2D::MakeBL(int faceid)
 {
     map<int, Array<OneD, NekDouble> > edgeNormals;
-
     int eid = 0;
-
     for (vector<unsigned>::iterator it = m_blCurves.begin();
          it != m_blCurves.end(); ++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)
@@ -218,30 +216,27 @@ void Generator2D::MakeBL(int faceid)
             Array<OneD, NekDouble> p1, p2;
             p1 = es[j]->m_n1->GetCADSurfInfo(faceid);
             p2 = es[j]->m_n2->GetCADSurfInfo(faceid);
+            Array<OneD, NekDouble> n(2);
+            n[0] = p1[1] - p2[1];
+            n[1] = p2[0] - p1[0];
             if (edgeo == CADOrientation::eBackwards)
             {
-                swap(p1, p2);
+                n[0] *= -1.0;
+                n[1] *= -1.0;
             }
-            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> np(2);
+            np[0] = p1[0] + n[0];
+            np[1] = p1[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;
         }
     }
@@ -260,21 +255,18 @@ void Generator2D::MakeBL(int faceid)
     map<NodeSharedPtr, vector<EdgeSharedPtr> >::iterator it;
     for (it = m_nodesToEdge.begin(); it != m_nodesToEdge.end(); it++)
     {
-        Array<OneD, NekDouble> n(3);
+        Array<OneD, NekDouble> n(3, 0.0);
         ASSERTL0(it->second.size() == 2,
                  "wierdness, most likely bl_surfs are incorrect");
         Array<OneD, NekDouble> n1 = edgeNormals[it->second[0]->m_id];
         Array<OneD, NekDouble> n2 = edgeNormals[it->second[1]->m_id];
-
         n[0]          = (n1[0] + n2[0]) / 2.0;
         n[1]          = (n1[1] + n2[1]) / 2.0;
         NekDouble mag = sqrt(n[0] * n[0] + n[1] * n[1]);
         n[0] /= mag;
         n[1] /= mag;
-
         NekDouble t = m_thickness.Evaluate(m_thickness_ID, it->first->m_x,
                                            it->first->m_y, 0.0, 0.0);
-
         // Adjust thickness according to angle between normals
         if (adjust)
         {
@@ -286,10 +278,8 @@ void Generator2D::MakeBL(int faceid)
             }
         }
 
-        n[0] = n[0] * t + it->first->m_x;
-        n[1] = n[1] * t + it->first->m_y;
-        n[2] = 0.0;
-
+        n[0]             = n[0] * t + it->first->m_x;
+        n[1]             = n[1] * t + it->first->m_y;
         NodeSharedPtr nn = boost::shared_ptr<Node>(
             new Node(m_mesh->m_numNodes++, n[0], n[1], 0.0));
         CADSurfSharedPtr s = m_mesh->m_cad->GetSurf(faceid);
@@ -303,7 +293,6 @@ void Generator2D::MakeBL(int faceid)
     {
         CADOrientation::Orientation edgeo =
             m_mesh->m_cad->GetCurve(*it)->GetOrienationWRT(faceid);
-
         vector<NodeSharedPtr> ns = m_curvemeshes[*it]->GetMeshPoints();
         vector<NodeSharedPtr> newNs;
         for (int i = 0; i < ns.size(); i++)
@@ -312,7 +301,6 @@ void Generator2D::MakeBL(int faceid)
         }
         m_curvemeshes[*it] =
             MemoryManager<CurveMesh>::AllocateSharedPtr(*it, m_mesh, newNs);
-
         if (edgeo == CADOrientation::eBackwards)
         {
             reverse(ns.begin(), ns.end());
@@ -320,22 +308,16 @@ void Generator2D::MakeBL(int faceid)
         for (int i = 0; i < ns.size() - 1; ++i)
         {
             vector<NodeSharedPtr> qns;
-
             qns.push_back(ns[i]);
             qns.push_back(ns[i + 1]);
             qns.push_back(nodeNormals[ns[i + 1]]);
             qns.push_back(nodeNormals[ns[i]]);
-
             ElmtConfig conf(LibUtilities::eQuadrilateral, 1, false, false);
-
             vector<int> tags;
             tags.push_back(101);
-
             ElementSharedPtr E = GetElementFactory().CreateInstance(
                 LibUtilities::eQuadrilateral, conf, qns, tags);
-
             E->m_parentCAD = m_mesh->m_cad->GetSurf(faceid);
-
             for (int j = 0; j < E->GetEdgeCount(); ++j)
             {
                 pair<EdgeSet::iterator, bool> testIns;
@@ -353,6 +335,60 @@ void Generator2D::MakeBL(int faceid)
     }
 }
 
+void Generator2D::PeriodicPrep()
+{
+    m_periodicPairs.clear();
+    set<unsigned> periodic;
+
+    // Build periodic curve pairs
+    string s = m_config["periodic"].as<string>();
+    vector<string> lines;
+
+    boost::split(lines, s, boost::is_any_of(":"));
+
+    for (vector<string>::iterator il = lines.begin(); il != lines.end(); ++il)
+    {
+        vector<string> tmp;
+        boost::split(tmp, *il, boost::is_any_of(","));
+
+        ASSERTL0(tmp.size() == 2, "periodic pairs ill-defined");
+
+        vector<unsigned> data(2);
+        data[0] = boost::lexical_cast<unsigned>(tmp[0]);
+        data[1] = boost::lexical_cast<unsigned>(tmp[1]);
+
+        ASSERTL0(!periodic.count(data[0]), "curve already periodic");
+        ASSERTL0(!periodic.count(data[1]), "curve already periodic");
+
+        m_periodicPairs[data[0]] = data[1];
+        periodic.insert(data[0]);
+        periodic.insert(data[1]);
+    }
+}
+
+void Generator2D::MakePeriodic()
+{
+    // Override slave curves
+
+    for (map<unsigned, unsigned>::iterator ip = m_periodicPairs.begin();
+         ip != m_periodicPairs.end(); ++ip)
+    {
+        m_curvemeshes[ip->second]->PeriodicOverwrite(m_curvemeshes[ip->first]);
+    }
+
+    if (m_mesh->m_verbose)
+    {
+        cout << "\t\tPeriodic boundary conditions" << endl;
+        for (map<unsigned, unsigned>::iterator it = m_periodicPairs.begin();
+             it != m_periodicPairs.end(); ++it)
+        {
+            cout << "\t\t\tCurves " << it->first << " => " << it->second
+                 << endl;
+        }
+        cout << endl;
+    }
+}
+
 void Generator2D::Report()
 {
     if (m_mesh->m_verbose)

library/NekMeshUtils/2DGenerator/2DGenerator.h

@@ -61,14 +61,18 @@ public:
     static ModuleKey className;
 
     Generator2D(MeshSharedPtr m);
+
     virtual ~Generator2D();
 
     virtual void Process();
 
 private:
-
     void MakeBLPrep();
 
+    void PeriodicPrep();
+
+    void MakePeriodic();
+
     void MakeBL(int faceid);
 
     void Report();
@@ -76,6 +80,8 @@ private:
     std::map<int, FaceMeshSharedPtr> m_facemeshes;
     /// map of individual curve meshes of the curves in the domain
     std::map<int, CurveMeshSharedPtr> m_curvemeshes;
+    /// map of periodic curve pairs
+    std::map<unsigned, unsigned> m_periodicPairs;
 
     std::vector<unsigned int> m_blCurves;
     LibUtilities::AnalyticExpressionEvaluator m_thickness;

library/NekMeshUtils/CADSystem/CADSystem.cpp

@@ -33,10 +33,10 @@
 //
 ////////////////////////////////////////////////////////////////////////////////
 
-#include <NekMeshUtils/CADSystem/CADSystem.h>
-#include <NekMeshUtils/CADSystem/CADVert.h>
 #include <NekMeshUtils/CADSystem/CADCurve.h>
 #include <NekMeshUtils/CADSystem/CADSurf.h>
+#include <NekMeshUtils/CADSystem/CADSystem.h>
+#include <NekMeshUtils/CADSystem/CADVert.h>
 
 using namespace std;
 
@@ -45,7 +45,7 @@ namespace Nektar
 namespace NekMeshUtils
 {
 
-EngineFactory& GetEngineFactory()
+EngineFactory &GetEngineFactory()
 {
     typedef Loki::SingletonHolder<EngineFactory, Loki::CreateUsingNew,
                                   Loki::NoDestroy, Loki::SingleThreaded>
@@ -53,7 +53,7 @@ EngineFactory& GetEngineFactory()
     return Type::Instance();
 }
 
-CADVertFactory& GetCADVertFactory()
+CADVertFactory &GetCADVertFactory()
 {
     typedef Loki::SingletonHolder<CADVertFactory, Loki::CreateUsingNew,
                                   Loki::NoDestroy, Loki::SingleThreaded>
@@ -61,7 +61,7 @@ CADVertFactory& GetCADVertFactory()
     return Type::Instance();
 }
 
-CADCurveFactory& GetCADCurveFactory()
+CADCurveFactory &GetCADCurveFactory()
 {
     typedef Loki::SingletonHolder<CADCurveFactory, Loki::CreateUsingNew,
                                   Loki::NoDestroy, Loki::SingleThreaded>
@@ -69,7 +69,7 @@ CADCurveFactory& GetCADCurveFactory()
     return Type::Instance();
 }
 
-CADSurfFactory& GetCADSurfFactory()
+CADSurfFactory &GetCADSurfFactory()
 {
     typedef Loki::SingletonHolder<CADSurfFactory, Loki::CreateUsingNew,
                                   Loki::NoDestroy, Loki::SingleThreaded>
@@ -77,5 +77,37 @@ CADSurfFactory& GetCADSurfFactory()
     return Type::Instance();
 }
 
+Array<OneD, NekDouble> CADSystem::GetPeriodicTranslationVector(int first,
+                                                               int second)
+{
+    ASSERTL0(GetNumSurf() == 1, "wont work for multi surfaces yet");
+
+    CADCurveSharedPtr c1 = GetCurve(first);
+    CADCurveSharedPtr c2 = GetCurve(second);
+
+    NekDouble tst = c1->GetTotLength() - c2->GetTotLength();
+    ASSERTL0(fabs(tst) < 1e-6, "periodic curves not same length");
+
+    vector<CADVertSharedPtr> v1 = c1->GetVertex();
+    Array<OneD, NekDouble> p1 = v1[0]->GetLoc();
+
+    Array<OneD, NekDouble> p2;
+    vector<CADVertSharedPtr> v2 = c2->GetVertex();
+    if (c1->GetOrienationWRT(1) == c2->GetOrienationWRT(1))
+    {
+        p2 = v2[1]->GetLoc();
+    }
+    else
+    {
+        p2 = v2[0]->GetLoc();
+    }
+
+    Array<OneD, NekDouble> ret(3);
+    ret[0] = p2[0] - p1[0];
+    ret[1] = p2[1] - p1[1];
+    ret[2] = p2[2] - p1[2];
+
+    return ret;
+}
 }
 }

library/NekMeshUtils/CADSystem/CADSystem.h

@@ -43,6 +43,8 @@
 
 #include <LibUtilities/BasicUtils/NekFactory.hpp>
 
+#include <NekMeshUtils/NekMeshUtilsDeclspec.h>
+
 #include "CADObject.h"
 
 namespace Nektar
@@ -201,6 +203,9 @@ public:
         return m_verts.size();
     }
 
+    NEKMESHUTILS_EXPORT Array<OneD, NekDouble> GetPeriodicTranslationVector(
+                int first, int second);
+
 protected:
     /// Name of cad file
     std::string m_name;

library/NekMeshUtils/CADSystem/OCE/CADSurfOCE.cpp

@@ -61,8 +61,8 @@ void CADSurfOCE::Initialise(int i, TopoDS_Shape in)
     gp_Pnt ori(0.0, 0.0, 0.0);
     transform.SetScale(ori, 1.0 / 1000.0);
     TopLoc_Location mv(transform);
-
     in.Move(mv);
+
     m_occSurface = BRepAdaptor_Surface(TopoDS::Face(in));
     m_id         = i;
 
@@ -85,14 +85,14 @@ Array<OneD, NekDouble> CADSurfOCE::locuv(Array<OneD, NekDouble> p)
 
     Array<OneD, NekDouble> uvr(2);
 
-    gp_Pnt2d p2 = m_sas->ValueOfUV(loc, 1e-3);
+    gp_Pnt2d p2 = m_sas->ValueOfUV(loc, Precision::Confusion());
     uvr[0]      = p2.X();
     uvr[1]      = p2.Y();
 
     gp_Pnt p3 = m_sas->Value(p2);
-    if (p3.Distance(loc) > 1.0)
+    if (p3.Distance(loc) > 1e-6)
     {
-        cout << "large locuv distance " << p3.Distance(loc) << " " << m_id
+        cout << "large locuv distance " << p3.Distance(loc)/1000.0 << " " << m_id
              << endl;
     }
 
@@ -173,13 +173,9 @@ NekDouble CADSurfOCE::DistanceTo(Array<OneD, NekDouble> p)
 {
     gp_Pnt loc(p[0] * 1000.0, p[1] * 1000.0, p[2] * 1000.0);
 
-    // alternative locuv methods
-    ShapeAnalysis_Surface sas(m_s);
-    sas.SetDomain(m_bounds[0], m_bounds[1], m_bounds[2], m_bounds[3]);
-
-    gp_Pnt2d p2 = sas.ValueOfUV(loc, 1e-7);
+    gp_Pnt2d p2 = m_sas->ValueOfUV(loc, Precision::Confusion());
 
-    gp_Pnt p3 = sas.Value(p2);
+    gp_Pnt p3 = m_sas->Value(p2);
 
     return p3.Distance(loc);
 }
@@ -189,13 +185,9 @@ void CADSurfOCE::ProjectTo(Array<OneD, NekDouble> &tp,
 {
     gp_Pnt loc(tp[0] * 1000.0, tp[1] * 1000.0, tp[2] * 1000.0);
 
-    // alternative locuv methods
-    ShapeAnalysis_Surface sas(m_s);
-    sas.SetDomain(m_bounds[0], m_bounds[1], m_bounds[2], m_bounds[3]);
-
-    gp_Pnt2d p2 = sas.ValueOfUV(loc, 1e-7);
+    gp_Pnt2d p2 = m_sas->ValueOfUV(loc, Precision::Confusion());
 
-    gp_Pnt p3 = sas.Value(p2);
+    gp_Pnt p3 = m_sas->Value(p2);
 
     tp[0] = p3.X() / 1000.0;
     tp[1] = p3.Y() / 1000.0;
@@ -226,7 +218,7 @@ Array<OneD, NekDouble> CADSurfOCE::N(Array<OneD, NekDouble> uv)
     Test(uv);
 #endif
 
-    BRepLProp_SLProps slp(m_occSurface, 2, 1e-6);
+    BRepLProp_SLProps slp(m_occSurface, 2, 1e-8);
     slp.SetParameters(uv[0], uv[1]);
 
     if (!slp.IsNormalDefined())

library/NekMeshUtils/CADSystem/OCE/OpenCascade.h

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

library/NekMeshUtils/MeshElements/Quadrilateral.cpp

@@ -105,6 +105,7 @@ Quadrilateral::Quadrilateral(ElmtConfig pConf,
         if (sum > 0.0)
         {
             reverse(m_edge.begin(), m_edge.end());
+            swap(m_vertex[1], m_vertex[3]);
         }
     }
 

library/NekMeshUtils/SurfaceMeshing/CurveMesh.cpp

@@ -353,5 +353,71 @@ void CurveMesh::GetSampleFunction()
         m_dst[i] = dsti;
     }
 }
+
+void CurveMesh::PeriodicOverwrite(CurveMeshSharedPtr from)
+{
+    //clear current mesh points and remove edges from edgeset
+    m_meshpoints.clear();
+    for (int i = 0; i < m_meshedges.size(); i++)
+    {
+        m_mesh->m_edgeSet.erase(m_meshedges[i]);
+    }
+    m_meshedges.clear();
+
+    ///////
+
+    int tid = from->GetId();
+    Array<OneD, NekDouble> T =
+        m_mesh->m_cad->GetPeriodicTranslationVector(tid, m_id);
+
+    CADCurveSharedPtr c1 = m_mesh->m_cad->GetCurve(tid);
+
+    bool reversed = c1->GetOrienationWRT(1) == m_cadcurve->GetOrienationWRT(1);
+
+    vector<NodeSharedPtr> nodes = from->GetMeshPoints();
+
+    vector<pair<CADSurfSharedPtr, CADOrientation::Orientation> > surfs =
+        m_cadcurve->GetAdjSurf();
+
+    for (int i = 1; i < nodes.size() - 1; i++)
+    {
+        Array<OneD, NekDouble> loc = nodes[i]->GetLoc();
+        NodeSharedPtr nn = NodeSharedPtr(new Node(
+            m_mesh->m_numNodes++, loc[0] + T[0], loc[1] + T[1], 0.0));
+
+        for (int j = 0; j < surfs.size(); j++)
+        {
+            nn->SetCADSurf(surfs[j].first->GetId(), surfs[j].first,
+                           surfs[j].first->locuv(nn->GetLoc()));
+        }
+
+        nn->SetCADCurve(m_id, m_cadcurve, m_cadcurve->loct(nn->GetLoc()));
+
+        m_meshpoints.push_back(nn);
+    }
+
+    // Reverse internal nodes of the vector if necessary
+    if (reversed)
+    {
+        reverse(m_meshpoints.begin(), m_meshpoints.end());
+    }
+
+    vector<CADVertSharedPtr> verts = m_cadcurve->GetVertex();
+
+    m_meshpoints.insert(m_meshpoints.begin(), verts[0]->GetNode());
+    m_meshpoints.push_back(verts[1]->GetNode());
+    //dont need to realign cad for vertices
+
+    // make edges and add them to the edgeset for the face mesher to use
+    for (int i = 0; i < m_meshpoints.size() - 1; i++)
+    {
+        EdgeSharedPtr e = boost::shared_ptr<Edge>(
+            new Edge(m_meshpoints[i], m_meshpoints