Fix prism output up to order 4 (but no higher right now), fix the use of face...

Fix prism output up to order 4 (but no higher right now), fix the use of face IDs for generic 3D elements

library/NekMeshUtils/MeshElements/Element.h

@@ -471,6 +471,13 @@ public:
         return StdRegions::eNoOrientation;
     }
 
+    NEKMESHUTILS_EXPORT virtual int GetFaceVertex(int i, int j)
+    {
+        ASSERTL0(false,
+                 "This function should be implemented at a shape level.");
+        return 0;
+    }
+
     NEKMESHUTILS_EXPORT void Print()
     {
         int i, j;

library/NekMeshUtils/MeshElements/Hexahedron.cpp

@@ -47,6 +47,12 @@ LibUtilities::ShapeType Hexahedron::m_type =
     GetElementFactory().RegisterCreatorFunction(
         LibUtilities::eHexahedron, Hexahedron::create, "Hexahedron");
 
+/// Vertex IDs that make up hexahedron faces.
+int Hexahedron::m_faceIds[6][4] = {
+    {0, 1, 2, 3}, {0, 1, 5, 4}, {1, 2, 6, 5},
+    {3, 2, 6, 7}, {0, 3, 7, 4}, {4, 5, 6, 7}
+};
+
 /**
  * @brief Create a hexahedral element.
  */
@@ -102,12 +108,6 @@ Hexahedron::Hexahedron(ElmtConfig pConf,
 
     // Create faces
     int face_edges[6][4];
-    int face_ids[6][4] = {{0, 1, 2, 3},
-                          {0, 1, 5, 4},
-                          {1, 2, 6, 5},
-                          {3, 2, 6, 7},
-                          {0, 3, 7, 4},
-                          {4, 5, 6, 7}};
     for (int j = 0; j < 6; ++j)
     {
         vector<NodeSharedPtr> faceVertices;
@@ -115,9 +115,9 @@ Hexahedron::Hexahedron(ElmtConfig pConf,
         vector<NodeSharedPtr> faceNodes;
         for (int k = 0; k < 4; ++k)
         {
-            faceVertices.push_back(m_vertex[face_ids[j][k]]);
-            NodeSharedPtr a = m_vertex[face_ids[j][k]];
-            NodeSharedPtr b = m_vertex[face_ids[j][(k + 1) % 4]];
+            faceVertices.push_back(m_vertex[m_faceIds[j][k]]);
+            NodeSharedPtr a = m_vertex[m_faceIds[j][k]];
+            NodeSharedPtr b = m_vertex[m_faceIds[j][(k + 1) % 4]];
             for (unsigned int i = 0; i < m_edge.size(); ++i)
             {
                 if (((*(m_edge[i]->m_n1) == *a) &&

library/NekMeshUtils/MeshElements/Hexahedron.h

@@ -87,6 +87,13 @@ public:
         int                               &id);
 
     NEKMESHUTILS_EXPORT static unsigned int GetNumNodes(ElmtConfig pConf);
+    NEKMESHUTILS_EXPORT virtual int GetFaceVertex(int i, int j)
+    {
+        return m_faceIds[i][j];
+    }
+
+private:
+    static int m_faceIds[6][4];
 };
 }
 }

library/NekMeshUtils/MeshElements/Mesh.cpp

@@ -109,12 +109,12 @@ void Mesh::MakeOrder(int                      order,
     }
     else if (distType == LibUtilities::eGaussLobattoLegendre)
     {
-        pTypes[LibUtilities::eSegment]       = LibUtilities::ePolyEvenlySpaced;
-        pTypes[LibUtilities::eTriangle]      = LibUtilities::eNodalTriEvenlySpaced;
-        pTypes[LibUtilities::eQuadrilateral] = LibUtilities::ePolyEvenlySpaced;
-        pTypes[LibUtilities::eTetrahedron]   = LibUtilities::eNodalTetEvenlySpaced;
-        pTypes[LibUtilities::ePrism]         = LibUtilities::eNodalPrismEvenlySpaced;
-        pTypes[LibUtilities::eHexahedron]    = LibUtilities::ePolyEvenlySpaced;
+        pTypes[LibUtilities::eSegment]       = LibUtilities::eGaussLobattoLegendre;
+        pTypes[LibUtilities::eTriangle]      = LibUtilities::eNodalTriElec;
+        pTypes[LibUtilities::eQuadrilateral] = LibUtilities::eGaussLobattoLegendre;
+        pTypes[LibUtilities::eTetrahedron]   = LibUtilities::eNodalTetElec;
+        // Prism still to do.
+        pTypes[LibUtilities::eHexahedron]    = LibUtilities::eGaussLobattoLegendre;
     }
 
     for(eit = m_edgeSet.begin(); eit != m_edgeSet.end(); eit++)

library/NekMeshUtils/MeshElements/Prism.cpp

@@ -49,6 +49,11 @@ LibUtilities::ShapeType Prism::m_type =
     GetElementFactory().RegisterCreatorFunction(
         LibUtilities::ePrism, Prism::create, "Prism");
 
+/// Vertex IDs that make up prism faces.
+int Prism::m_faceIds[5][4] = {
+    {0, 1, 2, 3}, {0, 1, 4, -1}, {1, 2, 5, 4}, {3, 2, 5, -1}, {0, 3, 5, 4}
+};
+
 /**
  * @brief Create a prism element.
  */
@@ -114,8 +119,6 @@ Prism::Prism(ElmtConfig pConf,
     }
 
     // Create faces
-    int face_ids[5][4] = {
-        {0, 1, 2, 3}, {0, 1, 4, -1}, {1, 2, 5, 4}, {3, 2, 5, -1}, {0, 3, 5, 4}};
     int face_edges[5][4];
 
     int face_offset[5];
@@ -135,9 +138,9 @@ Prism::Prism(ElmtConfig pConf,
 
         for (int k = 0; k < nEdge; ++k)
         {
-            faceVertices.push_back(m_vertex[face_ids[j][k]]);
-            NodeSharedPtr a = m_vertex[face_ids[j][k]];
-            NodeSharedPtr b = m_vertex[face_ids[j][(k + 1) % nEdge]];
+            faceVertices.push_back(m_vertex[m_faceIds[j][k]]);
+            NodeSharedPtr a = m_vertex[m_faceIds[j][k]];
+            NodeSharedPtr b = m_vertex[m_faceIds[j][(k + 1) % nEdge]];
             unsigned int i;
             for (i = 0; i < m_edge.size(); ++i)
             {
@@ -319,6 +322,8 @@ void Prism::MakeOrder(int                                order,
     const int nPrismIntPts = (nPoints - 2) * (nPoints - 3) * (nPoints - 2) / 2;
     m_volumeNodes.resize(nPrismIntPts);
 
+    cout << "numpoints = " << nPrismIntPts << endl;
+
     for (int i = nPrismPts - nPrismIntPts, cnt = 0; i < nPrismPts; ++i, ++cnt)
     {
         Array<OneD, NekDouble> xp(3);
@@ -326,6 +331,8 @@ void Prism::MakeOrder(int                                order,
         xp[1] = py[i];
         xp[2] = pz[i];
 
+        cout << xp[0] << " " << xp[1] << " " << xp[2] << endl;
+
         Array<OneD, NekDouble> x(3, 0.0);
         for (int j = 0; j < coordDim; ++j)
         {

library/NekMeshUtils/MeshElements/Prism.h

@@ -86,6 +86,10 @@ public:
         int                                coordDim,
         int                               &id);
     NEKMESHUTILS_EXPORT static unsigned int GetNumNodes(ElmtConfig pConf);
+    NEKMESHUTILS_EXPORT virtual int GetFaceVertex(int i, int j)
+    {
+        return m_faceIds[i][j];
+    }
 
     /**
      * Orientation of prism; unchanged = 0; clockwise = 1;
@@ -95,6 +99,9 @@ public:
 
 protected:
     void OrientPrism();
+
+private:
+    static int m_faceIds[5][4];
 };
 }
 }

library/NekMeshUtils/MeshElements/Pyramid.cpp

@@ -47,6 +47,11 @@ LibUtilities::ShapeType Pyramid::type =
     GetElementFactory().RegisterCreatorFunction(
         LibUtilities::ePyramid, Pyramid::create, "Pyramid");
 
+/// Vertex IDs that make up pyramid faces.
+int Pyramid::m_faceIds[5][4] = {
+    {0, 1, 2, 3}, {0, 1, 4, -1}, {1, 2, 4, -1}, {3, 2, 4, -1}, {0, 3, 4, -1}
+};
+
 /**
  * @brief Create a pyramidic element.
  */
@@ -98,11 +103,6 @@ Pyramid::Pyramid(ElmtConfig pConf,
     }
 
     // Create faces
-    int face_ids[5][4] = {{0, 1, 2, 3},
-                          {0, 1, 4, -1},
-                          {1, 2, 4, -1},
-                          {3, 2, 4, -1},
-                          {0, 3, 4, -1}};
     int face_edges[5][4];
     int faceoffset = 5 + 8 * n;
     for (int j = 0; j < 5; ++j)
@@ -114,9 +114,9 @@ Pyramid::Pyramid(ElmtConfig pConf,
 
         for (int k = 0; k < nEdge; ++k)
         {
-            faceVertices.push_back(m_vertex[face_ids[j][k]]);
-            NodeSharedPtr a = m_vertex[face_ids[j][k]];
-            NodeSharedPtr b = m_vertex[face_ids[j][(k + 1) % nEdge]];
+            faceVertices.push_back(m_vertex[m_faceIds[j][k]]);
+            NodeSharedPtr a = m_vertex[m_faceIds[j][k]];
+            NodeSharedPtr b = m_vertex[m_faceIds[j][(k + 1) % nEdge]];
             for (unsigned int i = 0; i < m_edge.size(); ++i)
             {
                 if ((m_edge[i]->m_n1 == a && m_edge[i]->m_n2 == b) ||

library/NekMeshUtils/MeshElements/Pyramid.h

@@ -77,11 +77,18 @@ public:
     NEKMESHUTILS_EXPORT virtual SpatialDomains::GeometrySharedPtr GetGeom(
         int coordDim);
     NEKMESHUTILS_EXPORT static unsigned int GetNumNodes(ElmtConfig pConf);
+    NEKMESHUTILS_EXPORT virtual int GetFaceVertex(int i, int j)
+    {
+        return m_faceIds[i][j];
+    }
 
     /**
      * Orientation of pyramid.
      */
     int orientationMap[5];
+
+private:
+    static int m_faceIds[5][4];
 };
 }
 }

library/NekMeshUtils/MeshElements/Tetrahedron.cpp

@@ -51,6 +51,11 @@ LibUtilities::ShapeType Tetrahedron::m_type =
     GetElementFactory().RegisterCreatorFunction(
         LibUtilities::eTetrahedron, Tetrahedron::create, "Tetrahedron");
 
+/// Vertex IDs that make up tetrahedron faces.
+int Tetrahedron::m_faceIds[4][3] = {
+    {0, 1, 2}, {0, 1, 3}, {1, 2, 3}, {0, 2, 3}
+};
+
 /**
  * @brief Create a tetrahedron element.
  */
@@ -116,9 +121,6 @@ Tetrahedron::Tetrahedron(ElmtConfig pConf,
         }
     }
 
-    // Face-vertex IDs
-    int face_ids[4][3] = {{0, 1, 2}, {0, 1, 3}, {1, 2, 3}, {0, 2, 3}};
-
     // Face-edge IDs
     int face_edges[4][3];
 
@@ -134,9 +136,9 @@ Tetrahedron::Tetrahedron(ElmtConfig pConf,
         // additional note below).
         for (int k = 0; k < 3; ++k)
         {
-            faceVertices.push_back(m_vertex[face_ids[j][k]]);
-            NodeSharedPtr a = m_vertex[face_ids[j][k]];
-            NodeSharedPtr b = m_vertex[face_ids[j][(k + 1) % 3]];
+            faceVertices.push_back(m_vertex[m_faceIds[j][k]]);
+            NodeSharedPtr a = m_vertex[m_faceIds[j][k]];
+            NodeSharedPtr b = m_vertex[m_faceIds[j][(k + 1) % 3]];
             for (unsigned int i = 0; i < m_edge.size(); ++i)
             {
                 if (((*(m_edge[i]->m_n1) == *a) &&
@@ -159,10 +161,10 @@ Tetrahedron::Tetrahedron(ElmtConfig pConf,
             const int nFaceNodes = n * (n - 1) / 2;
 
             // Get the vertex IDs of whatever face we are processing.
-            vector<int> faceIds(3);
-            faceIds[0] = faceVertices[0]->m_id;
-            faceIds[1] = faceVertices[1]->m_id;
-            faceIds[2] = faceVertices[2]->m_id;
+            vector<int> faceVertIds(3);
+            faceVertIds[0] = faceVertices[0]->m_id;
+            faceVertIds[1] = faceVertices[1]->m_id;
+            faceVertIds[2] = faceVertices[2]->m_id;
 
             // Find out the original face number as we were given it in
             // the constructor using the orientation map.
@@ -187,14 +189,14 @@ Tetrahedron::Tetrahedron(ElmtConfig pConf,
 
             // Find the original face vertex IDs.
             vector<int> origFaceIds(3);
-            origFaceIds[0] = pNodeList[face_ids[origFace][0]]->m_id;
-            origFaceIds[1] = pNodeList[face_ids[origFace][1]]->m_id;
-            origFaceIds[2] = pNodeList[face_ids[origFace][2]]->m_id;
+            origFaceIds[0] = pNodeList[m_faceIds[origFace][0]]->m_id;
+            origFaceIds[1] = pNodeList[m_faceIds[origFace][1]]->m_id;
+            origFaceIds[2] = pNodeList[m_faceIds[origFace][2]]->m_id;
 
             // Construct a HOTriangle object which performs the
             // orientation magically for us.
             HOTriangle<NodeSharedPtr> hoTri(origFaceIds, faceNodes);
-            hoTri.Align(faceIds);
+            hoTri.Align(faceVertIds);
 
             // Copy the face nodes back again.
             faceNodes = hoTri.surfVerts;
@@ -384,8 +386,6 @@ bool operator==(const struct TetOrient &a, const struct TetOrient &b)
  */
 void Tetrahedron::OrientTet()
 {
-    static int face_ids[4][3] = {{0, 1, 2}, {0, 1, 3}, {1, 2, 3}, {0, 2, 3}};
-
     TetOrientSet faces;
 
     // Create a copy of the original vertex ordering. This is used to
@@ -395,9 +395,9 @@ void Tetrahedron::OrientTet()
     {
         vector<int> nodes(3);
 
-        nodes[0] = m_vertex[face_ids[i][0]]->m_id;
-        nodes[1] = m_vertex[face_ids[i][1]]->m_id;
-        nodes[2] = m_vertex[face_ids[i][2]]->m_id;
+        nodes[0] = m_vertex[m_faceIds[i][0]]->m_id;
+        nodes[1] = m_vertex[m_faceIds[i][1]]->m_id;
+        nodes[2] = m_vertex[m_faceIds[i][2]]->m_id;
 
         sort(nodes.begin(), nodes.end());
         struct TetOrient faceNodes(nodes, i);
@@ -490,9 +490,9 @@ void Tetrahedron::OrientTet()
     {
         vector<int> nodes(3);
 
-        nodes[0] = m_vertex[face_ids[i][0]]->m_id;
-        nodes[1] = m_vertex[face_ids[i][1]]->m_id;
-        nodes[2] = m_vertex[face_ids[i][2]]->m_id;
+        nodes[0] = m_vertex[m_faceIds[i][0]]->m_id;
+        nodes[1] = m_vertex[m_faceIds[i][1]]->m_id;
+        nodes[2] = m_vertex[m_faceIds[i][2]]->m_id;
         sort(nodes.begin(), nodes.end());
 
         struct TetOrient faceNodes(nodes, 0);

library/NekMeshUtils/MeshElements/Tetrahedron.h

@@ -82,12 +82,19 @@ public:
         int                               &id);
 
     NEKMESHUTILS_EXPORT static unsigned int GetNumNodes(ElmtConfig pConf);
+    NEKMESHUTILS_EXPORT virtual int GetFaceVertex(int i, int j)
+    {
+        return m_faceIds[i][j];
+    }
 
     int m_orientationMap[4];
     int m_origVertMap[4];
 
 protected:
     void OrientTet();
+
+private:
+    static int m_faceIds[4][3];
 };
 }
 }

utilities/NekMesh/InputModules/InputGmsh.cpp

@@ -166,10 +166,40 @@ std::vector<int> triTensorNodeOrdering(const std::vector<int> &nodes, int n)
     return nodeList;
 }
 
+typedef boost::tuple<int, int, int> Mode;
+struct cmpop
+{
+    bool operator()(Mode const &a, Mode const &b) const
+    {
+        if (a.get<0>() < b.get<0>())
+        {
+            return true;
+        }
+        if (a.get<0>() > b.get<0>())
+        {
+            return false;
+        }
+        if (a.get<1>() < b.get<1>())
+        {
+            return true;
+        }
+        if (a.get<1>() > b.get<1>())
+        {
+            return false;
+        }
+        if (a.get<2>() < b.get<2>())
+        {
+            return true;
+        }
+
+        return false;
+    }
+};
+
+
 std::vector<int> tetTensorNodeOrdering(const std::vector<int> &nodes, int n)
 {
     std::vector<int> nodeList;
-    int cnt2;
     int nTri = n*(n+1)/2;
     int nTet = n*(n+1)*(n+2)/6;
 
@@ -193,35 +223,6 @@ std::vector<int> tetTensorNodeOrdering(const std::vector<int> &nodes, int n)
 
     // Set up a map that takes (a,b,c) -> m to help us figure out where things
     // are inside the tetrahedron.
-    typedef boost::tuple<int, int, int> Mode;
-    struct cmpop
-    {
-        bool operator()(Mode const &a, Mode const &b) const
-        {
-            if (a.get<0>() < b.get<0>())
-            {
-                return true;
-            }
-            if (a.get<0>() > b.get<0>())
-            {
-                return false;
-            }
-            if (a.get<1>() < b.get<1>())
-            {
-                return true;
-            }
-            if (a.get<1>() > b.get<1>())
-            {
-                return false;
-            }
-            if (a.get<2>() < b.get<2>())
-            {
-                return true;
-            }
-
-            return false;
-        }
-    };
     std::map<Mode, int, cmpop> tmp;
 
     for (int k = 0, cnt = 0; k < n; ++k)
@@ -344,9 +345,6 @@ std::vector<int> tetTensorNodeOrdering(const std::vector<int> &nodes, int n)
 std::vector<int> prismTensorNodeOrdering(const std::vector<int> &nodes, int n)
 {
     std::vector<int> nodeList;
-    int cnt2;
-    int nTri = n*(n+1)/2;
-    int nPrism = n*n*(n+1)/2;
 
     if (n == 0)
     {
@@ -355,34 +353,35 @@ std::vector<int> prismTensorNodeOrdering(const std::vector<int> &nodes, int n)
 
     nodeList.resize(nodes.size());
 
-    if (n == 1)
+    if (n == 2)
     {
         nodeList[0] = nodes[1];
         nodeList[1] = nodes[0];
         return nodeList;
     }
 
-    static int nekToGmshVerts[6] = {3, 4, 1, 0, 5, 2};
-
-    // Vertices
-    nodeList[0] = nodes[nekToGmshVerts[0]];
-    nodeList[1] = nodes[nekToGmshVerts[1]];
-    nodeList[2] = nodes[nekToGmshVerts[2]];
-    nodeList[3] = nodes[nekToGmshVerts[3]];
-    nodeList[4] = nodes[nekToGmshVerts[4]];
-    nodeList[5] = nodes[nekToGmshVerts[5]];
-    // nodeList[0]          = nodes[nekToGmshVerts[0]];
-    // nodeList[n - 1]      = nodes[nekToGmshVerts[1]];
-    // nodeList[n*n - 1]    = nodes[nekToGmshVerts[2]];
-    // nodeList[n*(n-1)]    = nodes[nekToGmshVerts[3]];
-    // nodeList[nPrism-1-n] = nodes[nekToGmshVerts[4]];
-    // nodeList[nPrism-1]   = nodes[nekToGmshVerts[5]];
-
-    if (n == 2)
+    // For some reason, this ordering is different. Whereas Gmsh usually orders
+    // vertices first, followed by edges, this ordering goes VVE-VVE-VVE for the
+    // three edges that contain edge-interior information, but also vertex
+    // orientation is not the same as the original Gmsh prism. The below is done
+    // by looking at the ordering by hand and is hence why we don't support
+    // order > 4 right now.
+    if (n == 3)
     {
-        return nodeList;
+        nodeList[0] = nodes[1];
+        nodeList[1] = nodes[4];
+        nodeList[2] = nodes[2];
+        nodeList[3] = nodes[5];
+        nodeList[4] = nodes[0];
+        nodeList[5] = nodes[3];
+        nodeList[6] = nodes[7];
+        nodeList[7] = nodes[8];
+        nodeList[8] = nodes[6];
     }
 
+    ASSERTL0(n < 4, "Prism Gmsh input and output is incomplete for orders "
+                    "larger than 4");
+
     return nodeList;
 }
 
@@ -1359,8 +1358,8 @@ vector<int> InputGmsh::PrismReordering(ElmtConfig conf)
         intPoints.push_back(i);
     }
 
-    // Reorder this stuff
-    tmp = prismTensorNodeOrdering(intPoints, order - 2);
+    // Reorder interior points
+    tmp = prismTensorNodeOrdering(intPoints, order - 1);
     mapping.insert(mapping.end(), tmp.begin(), tmp.end());
 
     return mapping;

utilities/NekMesh/OutputModules/OutputGmsh.cpp

@@ -56,7 +56,6 @@ ModuleKey OutputGmsh::className =
 OutputGmsh::OutputGmsh(MeshSharedPtr m) : OutputModule(m)
 {
     map<unsigned int, ElmtConfig> igelmap = InputGmsh::GenElmMap();
-
     map<unsigned int, ElmtConfig>::iterator it;
 
     // Populate #InputGmsh::elmMap and use this to construct an
@@ -238,25 +237,20 @@ void OutputGmsh::Process()
             for (int j = 0; j < faceList.size(); ++j)
             {
                 nodeList = faceList[j]->m_faceNodes;
+                int nFaceVerts = faceList[j]->m_vertexList.size();
+                vector<int> faceIds(nFaceVerts), volFaceIds(nFaceVerts);
 
-                if (faceList[j]->m_vertexList.size() == 3)
+                for (int k = 0; k < nFaceVerts; ++k)
                 {
-                    // TODO: move face_ids into a member function of Element
-                    int face_ids[5][4] = {
-                        {0, 1, 2, 3}, {0, 1, 4, -1}, {1, 2, 5, 4}, {3, 2, 5, -1}, {0, 3, 5, 4}};
-                    //int face_ids[4][3] = {
-                    //    {0, 1, 2}, {0, 1, 3}, {1, 2, 3}, {0, 2, 3}};
-                    vector<int> faceIds(3), volFaceIds(3);
-
-                    for (int k = 0; k < 3; ++k)
-                    {
-                        faceIds   [k] = faceList[j]->m_vertexList[k]->m_id;
-                        volFaceIds[k] = e->GetVertexList()[face_ids[j][k]]->m_id;
-                    }
+                    faceIds   [k] = faceList[j]->m_vertexList[k]->m_id;
+                    volFaceIds[k] =
+                        e->GetVertexList()[e->GetFaceVertex(j, k)]->m_id;
+                }
 
+                if (nFaceVerts == 3)
+                {
                     HOTriangle<NodeSharedPtr> hoTri(faceIds, nodeList);
                     hoTri.Align(volFaceIds);
-
                     for (int k = 0; k < hoTri.surfVerts.size(); ++k)
                     {
                         tags.push_back(hoTri.surfVerts[k]->m_id);
@@ -264,20 +258,6 @@ void OutputGmsh::Process()
                 }
                 else
                 {
-                    // TODO: move face_ids into a member function of Element
-                    int face_ids[5][4] = {
-                        {0, 1, 2, 3}, {0, 1, 4, -1}, {1, 2, 5, 4}, {3, 2, 5, -1}, {0, 3, 5, 4}};
-                    //int face_ids[6][4] = {{0, 1, 2, 3}, {0, 1, 5, 4},
-                    //                      {1, 2, 6, 5}, {3, 2, 6, 7},
-                    //                      {0, 3, 7, 4}, {4, 5, 6, 7}};
-                    vector<int> faceIds(4), volFaceIds(4);
-
-                    for (int k = 0; k < 4; ++k)
-                    {
-                        faceIds   [k] = faceList[j]->m_vertexList[k]->m_id;
-                        volFaceIds[k] = e->GetVertexList()[face_ids[j][k]]->m_id;
-                    }
-
                     HOQuadrilateral<NodeSharedPtr> hoQuad(faceIds, nodeList);
                     hoQuad.Align(volFaceIds);