Add Mesh::MakeOrder function to convert mesh to uniform polynomial order,...

Add Mesh::MakeOrder function to convert mesh to uniform polynomial order, working high-order triangle output from Gmsh

library/NekMeshUtils/MeshElements/Edge.h

@@ -38,6 +38,7 @@
 
 #include <iomanip>
 
+#include <LibUtilities/Foundations/ManagerAccess.h>
 #include <SpatialDomains/SegGeom.h>
 
 #include <NekMeshUtils/NekMeshUtilsDeclspec.h>
@@ -169,6 +170,44 @@ public:
         return ret;
     }
 
+    void MakeOrder(int                                order,
+                   SpatialDomains::GeometrySharedPtr  geom,
+                   LibUtilities::PointsType           edgeType,
+                   int                                coordDim,
+                   int                               &id)
+    {
+        int nPoints = order + 1;
+        StdRegions::StdExpansionSharedPtr xmap = geom->GetXmap();
+
+        Array<OneD, NekDouble> edgePoints;
+        LibUtilities::PointsKey edgeKey(nPoints, edgeType);
+        LibUtilities::PointsManager()[edgeKey]->GetPoints(edgePoints);
+
+        Array<OneD, Array<OneD, NekDouble> > phys(coordDim);
+
+        for (int i = 0; i < coordDim; ++i)
+        {
+            phys[i] = Array<OneD, NekDouble>(xmap->GetTotPoints());
+            xmap->BwdTrans(geom->GetCoeffs(i), phys[i]);
+        }
+
+        m_edgeNodes.resize(nPoints - 2);
+
+        for (int i = 1; i < nPoints - 1; ++i)
+        {
+            Array<OneD, NekDouble> x(3, 0.0);
+            for (int j = 0; j < coordDim; ++j)
+            {
+                x[j] = xmap->PhysEvaluate(edgePoints + i, phys[j]);
+            }
+
+            m_edgeNodes[i-1] = boost::shared_ptr<Node>(
+                new Node(id++, x[0], x[1], x[2]));
+        }
+
+        m_curveType = edgeType;
+    }
+
     /// ID of edge.
     unsigned int m_id;
     /// First vertex node.

library/NekMeshUtils/MeshElements/Element.h

@@ -448,13 +448,29 @@ public:
                  "This function should be implemented at a shape level.");
         return boost::shared_ptr<SpatialDomains::Geometry>();
     }
+
     NEKMESHUTILS_EXPORT int GetMaxOrder();
+
     /// Complete this object.
-    NEKMESHUTILS_EXPORT virtual void Complete(int order)
+    NEKMESHUTILS_EXPORT virtual void MakeOrder(
+        int                                order,
+        SpatialDomains::GeometrySharedPtr  geom,
+        LibUtilities::PointsType           edgeType,
+        int                                coordDim,
+        int                               &id)
+    {
+        ASSERTL0(false,
+                 "This function should be implemented at a shape level.");
+    }
+
+    NEKMESHUTILS_EXPORT virtual StdRegions::Orientation GetEdgeOrient(
+        int edgeId, EdgeSharedPtr edge)
     {
         ASSERTL0(false,
                  "This function should be implemented at a shape level.");
+        return StdRegions::eNoOrientation;
     }
+
     NEKMESHUTILS_EXPORT void Print()
     {
         int i, j;

library/NekMeshUtils/MeshElements/Mesh.cpp

@@ -34,6 +34,7 @@
 ////////////////////////////////////////////////////////////////////////////////
 
 #include <NekMeshUtils/MeshElements/Mesh.h>
+#include <LibUtilities/Foundations/ManagerAccess.h>
 
 using namespace std;
 
@@ -79,5 +80,88 @@ unsigned int Mesh::GetNumEntities()
 
     return nEnt;
 }
+
+/**
+ * @brief Convert this mesh into a high-order mesh.
+ */
+void Mesh::MakeOrder(int                      order,
+                     LibUtilities::PointsType distType)
+{
+    int nq = order + 1;
+    int id = m_vertexSet.size();
+
+    EdgeSet::iterator eit;
+    FaceSet::iterator fit;
+
+    boost::unordered_map<int, SpatialDomains::Geometry1DSharedPtr> edgeGeoms;
+    boost::unordered_map<int, SpatialDomains::Geometry2DSharedPtr> faceGeoms;
+    boost::unordered_map<int, SpatialDomains::GeometrySharedPtr> volGeoms;
+
+    // Decide on distribution of points to use for each shape type.
+    std::map<LibUtilities::ShapeType, LibUtilities::PointsType> pTypes;
+    if (distType == LibUtilities::ePolyEvenlySpaced)
+    {
+        pTypes[LibUtilities::eSegment]  = LibUtilities::ePolyEvenlySpaced;
+        pTypes[LibUtilities::eTriangle] = LibUtilities::eNodalTriEvenlySpaced;
+    }
+
+    for(eit = m_edgeSet.begin(); eit != m_edgeSet.end(); eit++)
+    {
+        SpatialDomains::Geometry1DSharedPtr geom =
+            (*eit)->GetGeom(m_spaceDim);
+        geom->FillGeom();
+        edgeGeoms[(*eit)->m_id] = geom;
+    }
+
+    for(fit = m_faceSet.begin(); fit != m_faceSet.end(); fit++)
+    {
+        SpatialDomains::Geometry2DSharedPtr geom =
+            (*fit)->GetGeom(m_spaceDim);
+        geom->FillGeom();
+        faceGeoms[(*fit)->m_id] = geom;
+    }
+
+    for(int i = 0; i < m_element[m_expDim].size(); i++)
+    {
+        ElementSharedPtr el = m_element[m_expDim][i];
+        SpatialDomains::GeometrySharedPtr geom =
+            el->GetGeom(m_spaceDim);
+        geom->FillGeom();
+        volGeoms[el->GetId()] = geom;
+    }
+
+    boost::unordered_set<int> processedEdges, processedFaces, processedVolumes;
+
+    for(eit = m_edgeSet.begin(); eit != m_edgeSet.end(); eit++)
+    {
+        int edgeId = (*eit)->m_id;
+
+        if (processedEdges.find(edgeId) != processedEdges.end())
+        {
+            continue;
+        }
+
+        (*eit)->MakeOrder(order, edgeGeoms[edgeId],
+                          pTypes[LibUtilities::eSegment], m_spaceDim, id);
+        processedEdges.insert(edgeId);
+    }
+
+    const int nElmt = m_element[m_expDim].size();
+    for (int i = 0; i < nElmt; ++i)
+    {
+        ElementSharedPtr el = m_element[m_expDim][i];
+        int elmtId = el->GetId();
+
+        if (processedVolumes.find(elmtId) != processedVolumes.end())
+        {
+            continue;
+        }
+
+        el->MakeOrder(order, volGeoms[elmtId], pTypes[el->GetConf().m_e],
+                      m_spaceDim, id);
+        processedEdges.insert(elmtId);
+    }
+}
+
 }
 }

library/NekMeshUtils/MeshElements/Mesh.h

@@ -132,6 +132,9 @@ public:
     NEKMESHUTILS_EXPORT unsigned int GetNumBndryElements();
     /// Returns the total number of entities in the mesh.
     NEKMESHUTILS_EXPORT unsigned int GetNumEntities();
+
+    void MakeOrder(int                      order,
+                   LibUtilities::PointsType distType);
 };
 /// Shared pointer to a mesh.
 typedef boost::shared_ptr<Mesh> MeshSharedPtr;

library/NekMeshUtils/MeshElements/Triangle.cpp

@@ -138,6 +138,25 @@ SpatialDomains::GeometrySharedPtr Triangle::GetGeom(int coordDim)
     return ret;
 }
 
+StdRegions::Orientation Triangle::GetEdgeOrient(int edgeId, EdgeSharedPtr edge)
+{
+    int locVert = edgeId;
+    if (edge->m_n1 == m_vertex[locVert])
+    {
+        return StdRegions::eForwards;
+    }
+    else if (edge->m_n2 == m_vertex[locVert])
+    {
+        return StdRegions::eBackwards;
+    }
+    else
+    {
+        ASSERTL1(false, "Edge is not connected to this triangle.");
+    }
+
+    return StdRegions::eNoOrientation;
+}
+
 /**
  * @brief Return the number of nodes defining a triangle.
  */
@@ -150,92 +169,59 @@ unsigned int Triangle::GetNumNodes(ElmtConfig pConf)
         return (n + 1) * (n + 2) / 2;
 }
 
-void Triangle::Complete(int order)
+void Triangle::MakeOrder(int                                order,
+                         SpatialDomains::GeometrySharedPtr  geom,
+                         LibUtilities::PointsType           pType,
+                         int                                coordDim,
+                         int                               &id)
 {
-    int i, j;
-
-    // Create basis key for a nodal tetrahedron.
-    LibUtilities::BasisKey B0(
-        LibUtilities::eOrtho_A,
-        order + 1,
-        LibUtilities::PointsKey(order + 1,
-                                LibUtilities::eGaussLobattoLegendre));
-    LibUtilities::BasisKey B1(
-        LibUtilities::eOrtho_B,
-        order + 1,
-        LibUtilities::PointsKey(order + 1,
-                                LibUtilities::eGaussRadauMAlpha1Beta0));
-
-    // Create a standard nodal triangle in order to get the
-    // Vandermonde matrix to perform interpolation to nodal points.
-    StdRegions::StdNodalTriExpSharedPtr nodalTri =
-        MemoryManager<StdRegions::StdNodalTriExp>::AllocateSharedPtr(
-            B0, B1, LibUtilities::eNodalTriEvenlySpaced);
-
-    SpatialDomains::TriGeomSharedPtr geom =
-        boost::dynamic_pointer_cast<SpatialDomains::TriGeom>(this->GetGeom(3));
-
-    // Create basis key for a triangle.
-    LibUtilities::BasisKey C0(
-        LibUtilities::eOrtho_A,
-        order + 1,
-        LibUtilities::PointsKey(order + 1,
-                                LibUtilities::eGaussLobattoLegendre));
-    LibUtilities::BasisKey C1(
-        LibUtilities::eOrtho_B,
-        order + 1,
-        LibUtilities::PointsKey(order + 1,
-                                LibUtilities::eGaussRadauMAlpha1Beta0));
-
-    // Create a triangle.
-    LocalRegions::TriExpSharedPtr tri =
-        MemoryManager<LocalRegions::TriExp>::AllocateSharedPtr(C0, C1, geom);
-
-    // Get coordinate array for tetrahedron.
-    int nqtot = tri->GetTotPoints();
-    Array<OneD, NekDouble> alloc(6 * nqtot);
-    Array<OneD, NekDouble> xi(alloc);
-    Array<OneD, NekDouble> yi(alloc + nqtot);
-    Array<OneD, NekDouble> zi(alloc + 2 * nqtot);
-    Array<OneD, NekDouble> xo(alloc + 3 * nqtot);
-    Array<OneD, NekDouble> yo(alloc + 4 * nqtot);
-    Array<OneD, NekDouble> zo(alloc + 5 * nqtot);
-    Array<OneD, NekDouble> tmp;
-
-    tri->GetCoords(xi, yi, zi);
-
-    for (i = 0; i < 3; ++i)
+    // Triangles of order < 3 have no interior volume points.
+    if (order < 3)
     {
-        Array<OneD, NekDouble> coeffs(nodalTri->GetNcoeffs());
-        tri->FwdTrans(alloc + i * nqtot, coeffs);
-        // Apply Vandermonde matrix to project onto nodal space.
-        nodalTri->ModalToNodal(coeffs, tmp = alloc + (i + 3) * nqtot);
+        m_volumeNodes.clear();
+        return;
     }
 
-    // Now extract points from the co-ordinate arrays into the
-    // edge/face/volume nodes. First, extract edge-interior nodes.
-    for (i = 0; i < 3; ++i)
+    int nPoints = order + 1;
+    StdRegions::StdExpansionSharedPtr xmap = geom->GetXmap();
+
+    Array<OneD, NekDouble> px, py;
+    LibUtilities::PointsKey pKey(nPoints, pType);
+    ASSERTL1(pKey.GetPointsDim() == 2, "Points distribution must be 2D");
+    LibUtilities::PointsManager()[pKey]->GetPoints(px, py);
+
+    Array<OneD, Array<OneD, NekDouble> > phys(coordDim);
+
+    for (int i = 0; i < coordDim; ++i)
     {
-        int pos = 3 + i * (order - 1);
-        m_edge[i]->m_edgeNodes.clear();
-        for (j = 0; j < order - 1; ++j)
-        {
-            m_edge[i]->m_edgeNodes.push_back(NodeSharedPtr(
-                new Node(0, xo[pos + j], yo[pos + j], zo[pos + j])));
-        }
+        phys[i] = Array<OneD, NekDouble>(xmap->GetTotPoints());
+        xmap->BwdTrans(geom->GetCoeffs(i), phys[i]);
     }
 
-    // Extract face-interior nodes.
-    int pos = 3 + 3 * (order - 1);
-    for (i = pos; i < (order + 1) * (order + 2) / 2; ++i)
+    const int nTriPts = nPoints * (nPoints + 1) / 2;
+    const int nTriIntPts = (nPoints - 3) * (nPoints - 2) / 2;
+    m_volumeNodes.resize(nTriIntPts);
+
+    for (int i = 3 + 3*(nPoints-2), cnt = 0; i < nTriPts; ++i, ++cnt)
     {
-        m_volumeNodes.push_back(
-            NodeSharedPtr(new Node(0, xo[i], yo[i], zo[i])));
+        Array<OneD, NekDouble> xp(2);
+        xp[0] = px[i];
+        xp[1] = py[i];
+
+        Array<OneD, NekDouble> x(3, 0.0);
+        for (int j = 0; j < coordDim; ++j)
+        {
+            x[j] = xmap->PhysEvaluate(xp, phys[j]);
+        }
+
+        m_volumeNodes[cnt] = boost::shared_ptr<Node>(
+            new Node(id++, x[0], x[1], x[2]));
     }
 
+    m_curveType          = pType;
     m_conf.m_order       = order;
     m_conf.m_faceNodes   = true;
-    m_conf.m_volumeNodes = true;
+    m_conf.m_volumeNodes = false;
 }
 }
 }

library/NekMeshUtils/MeshElements/Triangle.h

@@ -207,7 +207,14 @@ public:
 
     NEKMESHUTILS_EXPORT virtual SpatialDomains::GeometrySharedPtr GetGeom(
         int coordDim);
-    NEKMESHUTILS_EXPORT virtual void Complete(int order);
+    NEKMESHUTILS_EXPORT virtual StdRegions::Orientation GetEdgeOrient(
+        int edgeId, EdgeSharedPtr edge);
+    NEKMESHUTILS_EXPORT virtual void MakeOrder(
+        int                                order,
+        SpatialDomains::GeometrySharedPtr  geom,
+        LibUtilities::PointsType           pType,
+        int                                coordDim,
+        int                               &id);
 
     NEKMESHUTILS_EXPORT static unsigned int GetNumNodes(ElmtConfig pConf);
 };

utilities/NekMesh/InputModules/InputGmsh.h

@@ -66,15 +66,15 @@ public:
      * Element map; takes a msh id to an %ElmtConfig object.
      */
     static std::map<unsigned int, ElmtConfig> elmMap;
+    static std::vector<int> CreateReordering(unsigned int InputGmshEntity);
 
 private:
     int GetNnodes(unsigned int InputGmshEntity);
-    std::vector<int> CreateReordering(unsigned int InputGmshEntity);
-    std::vector<int> TriReordering(ElmtConfig conf);
-    std::vector<int> QuadReordering(ElmtConfig conf);
-    std::vector<int> HexReordering(ElmtConfig conf);
-    std::vector<int> PrismReordering(ElmtConfig conf);
-    std::vector<int> TetReordering(ElmtConfig conf);
+    static std::vector<int> TriReordering(ElmtConfig conf);
+    static std::vector<int> QuadReordering(ElmtConfig conf);
+    static std::vector<int> HexReordering(ElmtConfig conf);
+    static std::vector<int> PrismReordering(ElmtConfig conf);
+    static std::vector<int> TetReordering(ElmtConfig conf);
 };
 }
 }

utilities/NekMesh/OutputModules/OutputGmsh.cpp

@@ -110,96 +110,38 @@ void OutputGmsh::Process()
         }
     }
 
-    // maxOrder = 2;
-    for (int d = 1; d <= 3; ++d)
+    // Convert this mesh into a high-order mesh of uniform order.
+    cout << "Mesh order of " << maxOrder << " detected" << endl;
+    maxOrder = 6;
+    m_mesh->MakeOrder(maxOrder, LibUtilities::ePolyEvenlySpaced);
+
+    // Add edge- and face-interior nodes to vertex set.
+    EdgeSet::iterator eIt;
+    FaceSet::iterator fIt;
+
+    for (eIt = m_mesh->m_edgeSet.begin(); eIt != m_mesh->m_edgeSet.end(); ++eIt)
     {
-        for (int i = 0; i < m_mesh->m_element[d].size(); ++i)
-        {
-            ElementSharedPtr e = m_mesh->m_element[d][i];
-            if ((e->GetConf().m_order <= 1 && maxOrder > 1) ||
-                (e->GetConf().m_order == maxOrder &&
-                 e->GetConf().m_faceNodes == false))
-            {
-                toComplete.push_back(e);
-            }
-            // Generate geometry information for this element. This will
-            // be stored locally inside each element.
-            SpatialDomains::GeometrySharedPtr geom =
-                m_mesh->m_element[d][i]->GetGeom(m_mesh->m_spaceDim);
-        }
+        m_mesh->m_vertexSet.insert((*eIt)->m_edgeNodes.begin(),
+                                   (*eIt)->m_edgeNodes.end());
     }
 
-    // Complete these elements.
-    for (int i = 0; i < toComplete.size(); ++i)
+    for (fIt = m_mesh->m_faceSet.begin(); fIt != m_mesh->m_faceSet.end(); ++fIt)
     {
-        toComplete[i]->Complete(maxOrder);
+        m_mesh->m_vertexSet.insert((*fIt)->m_faceNodes.begin(),
+                                   (*fIt)->m_faceNodes.end());
     }
 
-    // Do second pass over elements to enumerate high-order vertices.
+    // Do second pass over elements for volume nodes.
     for (int d = 1; d <= 3; ++d)
     {
         for (int i = 0; i < m_mesh->m_element[d].size(); ++i)
         {
-            // Keep track of faces and edges to ensure that high-order
-            // nodes are only added once on common faces/edges.
-            boost::unordered_set<int> edgesDone;
-            boost::unordered_set<int> facesDone;
             ElementSharedPtr e = m_mesh->m_element[d][i];
+            vector<NodeSharedPtr> volList = e->GetVolumeNodes();
 
-            if (e->GetConf().m_order > 1)
+            for (int j = 0; j < volList.size(); ++j)
             {
-                vector<NodeSharedPtr> tmp;
-                vector<EdgeSharedPtr> edgeList = e->GetEdgeList();
-                vector<FaceSharedPtr> faceList = e->GetFaceList();
-                vector<NodeSharedPtr> volList  = e->GetVolumeNodes();
-
-                for (int j = 0; j < edgeList.size(); ++j)
-                {
-                    boost::unordered_set<int>::iterator it =
-                        edgesDone.find(edgeList[j]->m_id);
-                    if (it == edgesDone.end() || d != 3)
-                    {
-                        tmp.insert(tmp.end(),
-                                   edgeList[j]->m_edgeNodes.begin(),
-                                   edgeList[j]->m_edgeNodes.end());
-                        edgesDone.insert(edgeList[j]->m_id);
-                    }
-                }
-
-                for (int j = 0; j < faceList.size(); ++j)
-                {
-                    boost::unordered_set<int>::iterator it =
-                        facesDone.find(faceList[j]->m_id);
-                    if (it == facesDone.end() || d != 3)
-                    {
-                        tmp.insert(tmp.end(),
-                                   faceList[j]->m_faceNodes.begin(),
-                                   faceList[j]->m_faceNodes.end());
-                        facesDone.insert(faceList[j]->m_id);
-                    }
-                }
-
-                tmp.insert(tmp.end(), volList.begin(), volList.end());
-
-                // Even though faces/edges are at this point unique
-                // across the mesh, still need to test inserts since
-                // high-order nodes may already have been inserted into
-                // the list from an adjoining element or a boundary
-                // element.
-                for (int j = 0; j < tmp.size(); ++j)
-                {
-                    pair<NodeSet::iterator, bool> testIns =
-                        m_mesh->m_vertexSet.insert(tmp[j]);
-
-                    if (testIns.second)
-                    {
-                        (*(testIns.first))->m_id = id++;
-                    }
-                    else
-                    {
-                        tmp[j]->m_id = (*(testIns.first))->m_id;
-                    }
-                }
+                m_mesh->m_vertexSet.insert(volList[j]);
             }
         }
     }
@@ -214,7 +156,7 @@ void OutputGmsh::Process()
 
     for (it = tmp.begin(); it != tmp.end(); ++it)
     {
-        m_mshFile << (*it)->m_id << " " << scientific << setprecision(10)
+        m_mshFile << (*it)->m_id+1 << " " << scientific << setprecision(10)
                   << (*it)->m_x << " " << (*it)->m_y << " " << (*it)->m_z
                   << endl;
     }
@@ -226,7 +168,7 @@ void OutputGmsh::Process()
     m_mshFile << "$Elements" << endl;
     m_mshFile << m_mesh->GetNumEntities() << endl;
 
-    id = 0;
+    id = 1;
 
     for (int d = 1; d <= 3; ++d)
     {
@@ -235,7 +177,8 @@ void OutputGmsh::Process()
             ElementSharedPtr e = m_mesh->m_element[d][i];
 
             // First output element ID and type.
-            m_mshFile << id << " " << elmMap[e->GetConf()] << " ";
+            int elmtType = elmMap[e->GetConf()];
+            m_mshFile << id << " " << elmtType << " ";
 
             // Write out number of element tags and then the tags
             // themselves.
@@ -268,170 +211,54 @@ void OutputGmsh::Process()
                 tags.push_back(nodeList[j]->m_id);
             }
 
-            if (e->GetConf().m_order > 1)
+            // Process edge-interior points
+            for (int j = 0; j < edgeList.size(); ++j)
             {
-                for (int j = 0; j < edgeList.size(); ++j)
+                nodeList = edgeList[j]->m_edgeNodes;
+
+                if (e->GetEdgeOrient(j, edgeList[j]) == StdRegions::eForwards)
                 {
-                    nodeList = edgeList[j]->m_edgeNodes;
                     for (int k = 0; k < nodeList.size(); ++k)
                     {
                         tags.push_back(nodeList[k]->m_id);
-                        // cout << "EDGENODE" << endl;
                     }
                 }
-
-                for (int j = 0; j < faceList.size(); ++j)
+                else
                 {
-                    nodeList = faceList[j]->m_faceNodes;
-                    for (int k = 0; k < nodeList.size(); ++k)
+                    for (int k = nodeList.size() - 1; k >= 0; --k)
                     {
-                        // cout << "FACENODE" << endl;
                         tags.push_back(nodeList[k]->m_id);
                     }
                 }
-
-                for (int j = 0; j < volList.size(); ++j)
-                {
-                    // cout << "VOLNODE" << endl;
-                    tags.push_back(volList[j]->m_id);
-                }
             }
 
-            // Re-order tetrahedral vertices.
-            if (e->GetConf().m_e == LibUtilities::eTetrahedron)
+            // Process face-interior points
+            for (int j = 0; j < faceList.size(); ++j)
             {
-                int order = e->GetConf().m_order;
-                if (order > 4)
+                nodeList = faceList[j]->m_faceNodes;
+                for (int k = 0; k < nodeList.size(); ++k)
                 {
-                    cerr << "Temporary error: Gmsh tets only supported "