From 3cd8302cdeec66e6e86bdbcd7f724861e4d05058 Mon Sep 17 00:00:00 2001 From: David Moxey Date: Fri, 15 Jul 2016 07:57:10 +0100 Subject: [PATCH 01/25] Add Mesh::MakeOrder function to convert mesh to uniform polynomial order, working high-order triangle output from Gmsh --- library/NekMeshUtils/MeshElements/Edge.h | 39 +++ library/NekMeshUtils/MeshElements/Element.h | 18 +- library/NekMeshUtils/MeshElements/Mesh.cpp | 84 ++++++ library/NekMeshUtils/MeshElements/Mesh.h | 3 + .../NekMeshUtils/MeshElements/Triangle.cpp | 134 ++++----- library/NekMeshUtils/MeshElements/Triangle.h | 9 +- utilities/NekMesh/InputModules/InputGmsh.h | 12 +- .../NekMesh/OutputModules/OutputGmsh.cpp | 265 +++--------------- .../NekMesh/OutputModules/OutputNekpp.cpp | 5 +- 9 files changed, 267 insertions(+), 302 deletions(-) diff --git a/library/NekMeshUtils/MeshElements/Edge.h b/library/NekMeshUtils/MeshElements/Edge.h index 063eb9a9a..2c2d6c519 100644 --- a/library/NekMeshUtils/MeshElements/Edge.h +++ b/library/NekMeshUtils/MeshElements/Edge.h @@ -38,6 +38,7 @@ #include +#include #include #include @@ -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 edgePoints; + LibUtilities::PointsKey edgeKey(nPoints, edgeType); + LibUtilities::PointsManager()[edgeKey]->GetPoints(edgePoints); + + Array > phys(coordDim); + + for (int i = 0; i < coordDim; ++i) + { + phys[i] = Array(xmap->GetTotPoints()); + xmap->BwdTrans(geom->GetCoeffs(i), phys[i]); + } + + m_edgeNodes.resize(nPoints - 2); + + for (int i = 1; i < nPoints - 1; ++i) + { + Array 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( + new Node(id++, x[0], x[1], x[2])); + } + + m_curveType = edgeType; + } + /// ID of edge. unsigned int m_id; /// First vertex node. diff --git a/library/NekMeshUtils/MeshElements/Element.h b/library/NekMeshUtils/MeshElements/Element.h index 46cf74865..09048aa88 100644 --- a/library/NekMeshUtils/MeshElements/Element.h +++ b/library/NekMeshUtils/MeshElements/Element.h @@ -448,13 +448,29 @@ public: "This function should be implemented at a shape level."); return boost::shared_ptr(); } + 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; diff --git a/library/NekMeshUtils/MeshElements/Mesh.cpp b/library/NekMeshUtils/MeshElements/Mesh.cpp index fd1dd4aa5..f0784b586 100644 --- a/library/NekMeshUtils/MeshElements/Mesh.cpp +++ b/library/NekMeshUtils/MeshElements/Mesh.cpp @@ -34,6 +34,7 @@ //////////////////////////////////////////////////////////////////////////////// #include +#include 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 edgeGeoms; + boost::unordered_map faceGeoms; + boost::unordered_map volGeoms; + + // Decide on distribution of points to use for each shape type. + std::map 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 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); + } +} + } } diff --git a/library/NekMeshUtils/MeshElements/Mesh.h b/library/NekMeshUtils/MeshElements/Mesh.h index efb325204..6cd633ce6 100644 --- a/library/NekMeshUtils/MeshElements/Mesh.h +++ b/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 MeshSharedPtr; diff --git a/library/NekMeshUtils/MeshElements/Triangle.cpp b/library/NekMeshUtils/MeshElements/Triangle.cpp index 3bf09f50d..e465a45a4 100644 --- a/library/NekMeshUtils/MeshElements/Triangle.cpp +++ b/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::AllocateSharedPtr( - B0, B1, LibUtilities::eNodalTriEvenlySpaced); - - SpatialDomains::TriGeomSharedPtr geom = - boost::dynamic_pointer_cast(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::AllocateSharedPtr(C0, C1, geom); - - // Get coordinate array for tetrahedron. - int nqtot = tri->GetTotPoints(); - Array alloc(6 * nqtot); - Array xi(alloc); - Array yi(alloc + nqtot); - Array zi(alloc + 2 * nqtot); - Array xo(alloc + 3 * nqtot); - Array yo(alloc + 4 * nqtot); - Array zo(alloc + 5 * nqtot); - Array 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 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 px, py; + LibUtilities::PointsKey pKey(nPoints, pType); + ASSERTL1(pKey.GetPointsDim() == 2, "Points distribution must be 2D"); + LibUtilities::PointsManager()[pKey]->GetPoints(px, py); + + Array > 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(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 xp(2); + xp[0] = px[i]; + xp[1] = py[i]; + + Array x(3, 0.0); + for (int j = 0; j < coordDim; ++j) + { + x[j] = xmap->PhysEvaluate(xp, phys[j]); + } + + m_volumeNodes[cnt] = boost::shared_ptr( + 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; } } } diff --git a/library/NekMeshUtils/MeshElements/Triangle.h b/library/NekMeshUtils/MeshElements/Triangle.h index 65742f68e..b8f0009c7 100644 --- a/library/NekMeshUtils/MeshElements/Triangle.h +++ b/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); }; diff --git a/utilities/NekMesh/InputModules/InputGmsh.h b/utilities/NekMesh/InputModules/InputGmsh.h index 4b200f4a2..1ff000b06 100644 --- a/utilities/NekMesh/InputModules/InputGmsh.h +++ b/utilities/NekMesh/InputModules/InputGmsh.h @@ -66,15 +66,15 @@ public: * Element map; takes a msh id to an %ElmtConfig object. */ static std::map elmMap; + static std::vector CreateReordering(unsigned int InputGmshEntity); private: int GetNnodes(unsigned int InputGmshEntity); - std::vector CreateReordering(unsigned int InputGmshEntity); - std::vector TriReordering(ElmtConfig conf); - std::vector QuadReordering(ElmtConfig conf); - std::vector HexReordering(ElmtConfig conf); - std::vector PrismReordering(ElmtConfig conf); - std::vector TetReordering(ElmtConfig conf); + static std::vector TriReordering(ElmtConfig conf); + static std::vector QuadReordering(ElmtConfig conf); + static std::vector HexReordering(ElmtConfig conf); + static std::vector PrismReordering(ElmtConfig conf); + static std::vector TetReordering(ElmtConfig conf); }; } } diff --git a/utilities/NekMesh/OutputModules/OutputGmsh.cpp b/utilities/NekMesh/OutputModules/OutputGmsh.cpp index c3aed28d9..95779b459 100644 --- a/utilities/NekMesh/OutputModules/OutputGmsh.cpp +++ b/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 edgesDone; - boost::unordered_set facesDone; ElementSharedPtr e = m_mesh->m_element[d][i]; + vector volList = e->GetVolumeNodes(); - if (e->GetConf().m_order > 1) + for (int j = 0; j < volList.size(); ++j) { - vector tmp; - vector edgeList = e->GetEdgeList(); - vector faceList = e->GetFaceList(); - vector volList = e->GetVolumeNodes(); - - for (int j = 0; j < edgeList.size(); ++j) - { - boost::unordered_set::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::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 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 " - << "up to 4th order - will fix soon!" << endl; - abort(); - } - int pos = 4; - // Swap edge 1->3 nodes with edge 2->3 nodes. - pos = 4 + 4 * (order - 1); - for (int j = 0; j < order - 1; ++j) - { - swap(tags[j + pos], tags[j + pos + order - 1]); - } - // Reverse ordering of other vertical edge-interior - // nodes. - reverse(tags.begin() + 4 + 3 * (order - 1), - tags.begin() + 4 + 4 * (order - 1)); - reverse(tags.begin() + 4 + 4 * (order - 1), - tags.begin() + 4 + 5 * (order - 1)); - reverse(tags.begin() + 4 + 5 * (order - 1), - tags.begin() + 4 + 6 * (order - 1)); - - // Swap face 2 nodes with face 3. - pos = 4 + 6 * (order - 1) + 2 * (order - 2) * (order - 1) / 2; - for (int j = 0; j < (order - 2) * (order - 1) / 2; ++j) - { - swap(tags[j + pos], - tags[j + pos + (order - 2) * (order - 1) / 2]); - } - - // Re-order face points. Gmsh ordering (node->face) is: - // - // Face 0: 0->2->1 - // Face 1: 0->1->3 - // Face 2: 0->3->2 - // Face 3: 3->1->2 - // - // Therefore need to reorder nodes for faces 0, 2 and - // 3 to match nodal ordering. - - // Re-order face 0: transpose - vector tmp((order - 2) * (order - 1) / 2); - int a = 0; - pos = 4 + 6 * (order - 1); - for (int j = 0; j < order - 2; ++j) - { - for (int k = 0; k < order - 2 - j; ++k, ++a) - { - tmp[a] = - tags[pos + j + k * (2 * (order - 2) + 1 - k) / 2]; - } - } - for (int j = 0; j < (order - 1) * (order - 2) / 2; ++j) - { - tags[pos + j] = tmp[j]; - } - - // Re-order face 2: transpose - pos = 4 + 6 * (order - 1) + 2 * (order - 2) * (order - 1) / 2; - a = 0; - for (int j = 0; j < order - 2; ++j) - { - for (int k = 0; k < order - 2 - j; ++k, ++a) - { - tmp[a] = - tags[pos + j + k * (2 * (order - 2) + 1 - k) / 2]; - } - } - for (int j = 0; j < (order - 1) * (order - 2) / 2; ++j) - { - tags[pos + j] = tmp[j]; - } - - // Re-order face 3: reflect in y direction - pos = 4 + 6 * (order - 1) + 3 * (order - 2) * (order - 1) / 2; - a = 0; - for (int j = 0; j < order - 2; ++j) - { - for (int k = order - 3 - j; k >= 0; --k, ++a) - { - tmp[a] = - tags[pos + j + k * (2 * (order - 2) + 1 - k) / 2]; - } - } - - for (int j = 0; j < (order - 1) * (order - 2) / 2; ++j) - { - tags[pos + j] = tmp[j]; + tags.push_back(nodeList[k]->m_id); } } - // Re-order prism vertices. - else if (e->GetConf().m_e == LibUtilities::ePrism) + + // Process volume nodes + for (int j = 0; j < volList.size(); ++j) { - int order = e->GetConf().m_order; - if (order > 2) - { - cerr << "Temporary error: Gmsh prisms only " - << "supported up to 2nd order!" << endl; - abort(); - } + tags.push_back(volList[j]->m_id); + } - // Swap nodes. - vector temp(18); - temp[0] = tags[0]; - temp[1] = tags[1]; - temp[2] = tags[4]; - temp[3] = tags[3]; - temp[4] = tags[2]; - temp[5] = tags[5]; - temp[6] = tags[6]; - temp[7] = tags[10]; - temp[8] = tags[9]; - temp[9] = tags[11]; - temp[10] = tags[7]; - temp[11] = tags[14]; - temp[12] = tags[8]; - temp[13] = tags[13]; - temp[14] = tags[12]; - temp[15] = tags[15]; - temp[16] = tags[17]; - temp[17] = tags[16]; - for (int k = 0; k < 18; ++k) - { - tags[k] = temp[k]; - } + vector reordering = InputGmsh::CreateReordering(elmtType); + vector inv(tags.size()); + for (int j = 0; j < tags.size(); ++j) + { + inv[reordering[j]] = j; } // Finally write element nodes. for (int j = 0; j < tags.size(); ++j) { - m_mshFile << tags[j] << " "; + m_mshFile << tags[inv[j]] + 1 << " "; } m_mshFile << endl; diff --git a/utilities/NekMesh/OutputModules/OutputNekpp.cpp b/utilities/NekMesh/OutputModules/OutputNekpp.cpp index dd436a01d..4fb448049 100644 --- a/utilities/NekMesh/OutputModules/OutputNekpp.cpp +++ b/utilities/NekMesh/OutputModules/OutputNekpp.cpp @@ -98,6 +98,8 @@ void OutputNekpp::Process() cout << "OutputNekpp: Writing file..." << endl; } + m_mesh->MakeOrder(5, LibUtilities::ePolyEvenlySpaced); + TiXmlDocument doc; TiXmlDeclaration *decl = new TiXmlDeclaration("1.0", "utf-8", ""); doc.LinkEndChild(decl); @@ -679,7 +681,8 @@ void OutputNekpp::WriteXmlCurves(TiXmlElement *pRoot) } } // 2D elements in 3-space, output face curvature information - else if (m_mesh->m_expDim == 2 && m_mesh->m_spaceDim == 3) + else if (m_mesh->m_expDim == 2 && + m_mesh->m_spaceDim >= 2) { vector::iterator it; for (it = m_mesh->m_element[m_mesh->m_expDim].begin(); -- GitLab From 2014123e2c3c46ee1ed8dbe1e38ad265a75f587f Mon Sep 17 00:00:00 2001 From: David Moxey Date: Fri, 15 Jul 2016 08:08:02 +0100 Subject: [PATCH 02/25] Working quadrilateral --- library/NekMeshUtils/MeshElements/Mesh.cpp | 6 +- .../MeshElements/Quadrilateral.cpp | 108 ++++++++++-------- .../NekMeshUtils/MeshElements/Quadrilateral.h | 9 +- 3 files changed, 74 insertions(+), 49 deletions(-) diff --git a/library/NekMeshUtils/MeshElements/Mesh.cpp b/library/NekMeshUtils/MeshElements/Mesh.cpp index f0784b586..9ac53c5b2 100644 --- a/library/NekMeshUtils/MeshElements/Mesh.cpp +++ b/library/NekMeshUtils/MeshElements/Mesh.cpp @@ -87,7 +87,6 @@ unsigned int Mesh::GetNumEntities() void Mesh::MakeOrder(int order, LibUtilities::PointsType distType) { - int nq = order + 1; int id = m_vertexSet.size(); EdgeSet::iterator eit; @@ -101,8 +100,9 @@ void Mesh::MakeOrder(int order, std::map pTypes; if (distType == LibUtilities::ePolyEvenlySpaced) { - pTypes[LibUtilities::eSegment] = LibUtilities::ePolyEvenlySpaced; - pTypes[LibUtilities::eTriangle] = LibUtilities::eNodalTriEvenlySpaced; + pTypes[LibUtilities::eSegment] = LibUtilities::ePolyEvenlySpaced; + pTypes[LibUtilities::eTriangle] = LibUtilities::eNodalTriEvenlySpaced; + pTypes[LibUtilities::eQuadrilateral] = LibUtilities::ePolyEvenlySpaced; } for(eit = m_edgeSet.begin(); eit != m_edgeSet.end(); eit++) diff --git a/library/NekMeshUtils/MeshElements/Quadrilateral.cpp b/library/NekMeshUtils/MeshElements/Quadrilateral.cpp index 8a7fa70f8..f8b9471ae 100644 --- a/library/NekMeshUtils/MeshElements/Quadrilateral.cpp +++ b/library/NekMeshUtils/MeshElements/Quadrilateral.cpp @@ -114,63 +114,81 @@ Quadrilateral::Quadrilateral(ElmtConfig pConf, } } -void Quadrilateral::Complete(int order) +StdRegions::Orientation Quadrilateral::GetEdgeOrient( + int edgeId, EdgeSharedPtr edge) { - LibUtilities::BasisKey C0( - LibUtilities::eOrtho_A, - order + 1, - LibUtilities::PointsKey(order + 1, - LibUtilities::eGaussLobattoLegendre)); - - SpatialDomains::QuadGeomSharedPtr geom = - boost::dynamic_pointer_cast(this->GetGeom(3)); - - // Create a quad. - LocalRegions::QuadExpSharedPtr quad = - MemoryManager::AllocateSharedPtr(C0, C0, geom); - - // Get coordinate array for quadrilateral. - int nqtot = quad->GetTotPoints(); - Array alloc(3 * nqtot); - Array x(alloc); - Array y(alloc + 1 * nqtot); - Array z(alloc + 2 * nqtot); - - quad->GetCoords(x, y, z); - - // Now extract points from the co-ordinate arrays into the edge - // and face nodes. First, extract edge-interior nodes. - int edgeMap[4][2] = {{0, 1}, - {order, order + 1}, - {nqtot - 1, -1}, - {order * (order + 1), -order - 1}}; + 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 quadrilateral."); + } - for (int i = 0; i < 4; ++i) + return StdRegions::eNoOrientation; +} + +void Quadrilateral::MakeOrder(int order, + SpatialDomains::GeometrySharedPtr geom, + LibUtilities::PointsType pType, + int coordDim, + int &id) +{ + // Triangles of order < 2 have no interior volume points. + if (order < 2) { - int pos = edgeMap[i][0] + edgeMap[i][1]; - m_edge[i]->m_edgeNodes.clear(); - for (int j = 1; j < order; ++j, pos += edgeMap[i][1]) - { - m_edge[i]->m_edgeNodes.push_back( - NodeSharedPtr(new Node(0, x[pos], y[pos], z[pos]))); - } + m_volumeNodes.clear(); + return; + } + + int nPoints = order + 1; + StdRegions::StdExpansionSharedPtr xmap = geom->GetXmap(); + + Array px; + LibUtilities::PointsKey pKey(nPoints, pType); + ASSERTL1(pKey.GetPointsDim() == 1, "Points distribution must be 1D"); + LibUtilities::PointsManager()[pKey]->GetPoints(px); + + Array > phys(coordDim); + + for (int i = 0; i < coordDim; ++i) + { + phys[i] = Array(xmap->GetTotPoints()); + xmap->BwdTrans(geom->GetCoeffs(i), phys[i]); } - // Extract face-interior nodes. - m_volumeNodes.clear(); - for (int i = 1; i < order; ++i) + int nQuadIntPts = (nPoints - 2) * (nPoints - 2); + m_volumeNodes.resize(nQuadIntPts); + + for (int i = 1, cnt = 0; i < nPoints-1; ++i) { - int pos = i * (order + 1); - for (int j = 1; j < order; ++j) + for (int j = 1; j < nPoints-1; ++j, ++cnt) { - m_volumeNodes.push_back( - NodeSharedPtr(new Node(0, x[pos + j], y[pos + j], z[pos + j]))); + Array xp(2); + xp[0] = px[j]; + xp[1] = px[i]; + + Array x(3, 0.0); + for (int k = 0; k < coordDim; ++k) + { + x[k] = xmap->PhysEvaluate(xp, phys[k]); + } + + m_volumeNodes[cnt] = boost::shared_ptr( + 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; } SpatialDomains::GeometrySharedPtr Quadrilateral::GetGeom(int coordDim) diff --git a/library/NekMeshUtils/MeshElements/Quadrilateral.h b/library/NekMeshUtils/MeshElements/Quadrilateral.h index 3eaf73645..e0f728f79 100644 --- a/library/NekMeshUtils/MeshElements/Quadrilateral.h +++ b/library/NekMeshUtils/MeshElements/Quadrilateral.h @@ -76,7 +76,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); }; -- GitLab From ce2bd1ef2e3c5a969efefac23afba4969280c3d0 Mon Sep 17 00:00:00 2001 From: David Moxey Date: Fri, 15 Jul 2016 20:44:30 +0100 Subject: [PATCH 03/25] Working tets up to order 6 - tetTensorNodeOrdering needs to be completed --- library/NekMeshUtils/MeshElements/Face.h | 110 ++++++++++ library/NekMeshUtils/MeshElements/Mesh.cpp | 23 ++- .../MeshElements/Quadrilateral.cpp | 2 +- .../NekMeshUtils/MeshElements/Tetrahedron.cpp | 195 ++++++++---------- .../NekMeshUtils/MeshElements/Tetrahedron.h | 9 +- utilities/NekMesh/InputModules/InputGmsh.cpp | 161 ++++++++++++++- .../NekMesh/OutputModules/OutputGmsh.cpp | 30 ++- 7 files changed, 407 insertions(+), 123 deletions(-) diff --git a/library/NekMeshUtils/MeshElements/Face.h b/library/NekMeshUtils/MeshElements/Face.h index d4690b339..97c3ef74e 100644 --- a/library/NekMeshUtils/MeshElements/Face.h +++ b/library/NekMeshUtils/MeshElements/Face.h @@ -215,6 +215,116 @@ public: return s.str(); } + void MakeOrder(int order, + SpatialDomains::GeometrySharedPtr geom, + LibUtilities::PointsType pType, + int coordDim, + int &id) + { + if (m_vertexList.size() == 3) + { + // Triangles of order < 3 have no interior volume points. + if (order < 3) + { + m_faceNodes.clear(); + return; + } + + int nPoints = order + 1; + StdRegions::StdExpansionSharedPtr xmap = geom->GetXmap(); + + Array px, py; + LibUtilities::PointsKey pKey(nPoints, pType); + ASSERTL1(pKey.GetPointsDim() == 2, "Points distribution must be 2D"); + LibUtilities::PointsManager()[pKey]->GetPoints(px, py); + + Array > phys(coordDim); + + for (int i = 0; i < coordDim; ++i) + { + phys[i] = Array(xmap->GetTotPoints()); + xmap->BwdTrans(geom->GetCoeffs(i), phys[i]); + } + + const int nTriPts = nPoints * (nPoints + 1) / 2; + const int nTriIntPts = (nPoints - 3) * (nPoints - 2) / 2; + m_faceNodes.resize(nTriIntPts); + + for (int i = 3 + 3*(nPoints-2), cnt = 0; i < nTriPts; ++i, ++cnt) + { + Array xp(2); + xp[0] = px[i]; + xp[1] = py[i]; + + Array x(3, 0.0); + for (int j = 0; j < coordDim; ++j) + { + x[j] = xmap->PhysEvaluate(xp, phys[j]); + } + + m_faceNodes[cnt] = boost::shared_ptr( + new Node(id++, x[0], x[1], x[2])); + } + std::cout << "NFACEPTS = " << m_faceNodes.size() << std::endl; + + + m_curveType = pType; + } + else if (m_vertexList.size() == 4) + { + // Quads of order < 2 have no interior volume points. + if (order < 2) + { + m_faceNodes.clear(); + return; + } + + int nPoints = order + 1; + StdRegions::StdExpansionSharedPtr xmap = geom->GetXmap(); + + Array px; + LibUtilities::PointsKey pKey(nPoints, pType); + ASSERTL1(pKey.GetPointsDim() == 1, "Points distribution must be 1D"); + LibUtilities::PointsManager()[pKey]->GetPoints(px); + + Array > phys(coordDim); + + for (int i = 0; i < coordDim; ++i) + { + phys[i] = Array(xmap->GetTotPoints()); + xmap->BwdTrans(geom->GetCoeffs(i), phys[i]); + } + + int nQuadIntPts = (nPoints - 2) * (nPoints - 2); + m_faceNodes.resize(nQuadIntPts); + + for (int i = 1, cnt = 0; i < nPoints-1; ++i) + { + for (int j = 1; j < nPoints-1; ++j, ++cnt) + { + Array xp(2); + xp[0] = px[j]; + xp[1] = px[i]; + + Array x(3, 0.0); + for (int k = 0; k < coordDim; ++k) + { + x[k] = xmap->PhysEvaluate(xp, phys[k]); + } + + m_faceNodes[cnt] = boost::shared_ptr( + new Node(id++, x[0], x[1], x[2])); + } + } + + m_curveType = pType; + } + else + { + ASSERTL0(false, "Unknown number of vertices"); + } + } + /// Generate either SpatialDomains::TriGeom or /// SpatialDomains::QuadGeom for this element. SpatialDomains::Geometry2DSharedPtr GetGeom(int coordDim) diff --git a/library/NekMeshUtils/MeshElements/Mesh.cpp b/library/NekMeshUtils/MeshElements/Mesh.cpp index 9ac53c5b2..66e3dd963 100644 --- a/library/NekMeshUtils/MeshElements/Mesh.cpp +++ b/library/NekMeshUtils/MeshElements/Mesh.cpp @@ -103,6 +103,7 @@ void Mesh::MakeOrder(int order, pTypes[LibUtilities::eSegment] = LibUtilities::ePolyEvenlySpaced; pTypes[LibUtilities::eTriangle] = LibUtilities::eNodalTriEvenlySpaced; pTypes[LibUtilities::eQuadrilateral] = LibUtilities::ePolyEvenlySpaced; + pTypes[LibUtilities::eTetrahedron] = LibUtilities::eNodalTetEvenlySpaced; } for(eit = m_edgeSet.begin(); eit != m_edgeSet.end(); eit++) @@ -146,20 +147,28 @@ void Mesh::MakeOrder(int order, processedEdges.insert(edgeId); } - const int nElmt = m_element[m_expDim].size(); - for (int i = 0; i < nElmt; ++i) + for(fit = m_faceSet.begin(); fit != m_faceSet.end(); fit++) { - ElementSharedPtr el = m_element[m_expDim][i]; - int elmtId = el->GetId(); + int faceId = (*fit)->m_id; - if (processedVolumes.find(elmtId) != processedVolumes.end()) + if (processedFaces.find(faceId) != processedFaces.end()) { continue; } - el->MakeOrder(order, volGeoms[elmtId], pTypes[el->GetConf().m_e], + LibUtilities::ShapeType type = (*fit)->m_vertexList.size() == 3 ? + LibUtilities::eTriangle : LibUtilities::eQuadrilateral; + (*fit)->MakeOrder(order, faceGeoms[faceId], pTypes[type], m_spaceDim, + id); + processedFaces.insert(faceId); + } + + const int nElmt = m_element[m_expDim].size(); + for (int i = 0; i < nElmt; ++i) + { + ElementSharedPtr el = m_element[m_expDim][i]; + el->MakeOrder(order, volGeoms[el->GetId()], pTypes[el->GetConf().m_e], m_spaceDim, id); - processedEdges.insert(elmtId); } } diff --git a/library/NekMeshUtils/MeshElements/Quadrilateral.cpp b/library/NekMeshUtils/MeshElements/Quadrilateral.cpp index f8b9471ae..1b1d5b4ac 100644 --- a/library/NekMeshUtils/MeshElements/Quadrilateral.cpp +++ b/library/NekMeshUtils/MeshElements/Quadrilateral.cpp @@ -140,7 +140,7 @@ void Quadrilateral::MakeOrder(int order, int coordDim, int &id) { - // Triangles of order < 2 have no interior volume points. + // Quadrilaterals of order < 2 have no interior volume points. if (order < 2) { m_volumeNodes.clear(); diff --git a/library/NekMeshUtils/MeshElements/Tetrahedron.cpp b/library/NekMeshUtils/MeshElements/Tetrahedron.cpp index 27e66c931..bf44253c5 100644 --- a/library/NekMeshUtils/MeshElements/Tetrahedron.cpp +++ b/library/NekMeshUtils/MeshElements/Tetrahedron.cpp @@ -230,138 +230,111 @@ SpatialDomains::GeometrySharedPtr Tetrahedron::GetGeom(int coordDim) return ret; } -/** - * @brief Return the number of nodes defining a tetrahedron. - */ -unsigned int Tetrahedron::GetNumNodes(ElmtConfig pConf) +StdRegions::Orientation Tetrahedron::GetEdgeOrient( + int edgeId, EdgeSharedPtr edge) { - int n = pConf.m_order; - if (pConf.m_volumeNodes && pConf.m_faceNodes) - return (n + 1) * (n + 2) * (n + 3) / 6; - else if (!pConf.m_volumeNodes && pConf.m_faceNodes) - return 4 * (n + 1) * (n + 2) / 2 - 6 * (n + 1) + 4; + int locVert = edgeId; + int edgeVerts[6][2] = { {0,1}, {1,2}, {0,2}, {0,3}, {1,3}, {2,3} }; + + if (edge->m_n1 == m_vertex[edgeVerts[edgeId][0]]) + { + return StdRegions::eForwards; + } + else if (edge->m_n1 == m_vertex[edgeVerts[edgeId][1]]) + { + return StdRegions::eBackwards; + } else - return 6 * (n + 1) - 8; + { + ASSERTL1(false, "Edge is not connected to this quadrilateral."); + } + + return StdRegions::eNoOrientation; } -/** - * @brief . - */ -void Tetrahedron::Complete(int order) +void Tetrahedron::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)); - LibUtilities::BasisKey B2( - LibUtilities::eOrtho_C, - order + 1, - LibUtilities::PointsKey(order + 1, - LibUtilities::eGaussRadauMAlpha2Beta0)); - - // Create a standard nodal tetrahedron in order to get the - // Vandermonde matrix to perform interpolation to nodal points. - StdRegions::StdNodalTetExpSharedPtr nodalTet = - MemoryManager::AllocateSharedPtr( - B0, B1, B2, LibUtilities::eNodalTetEvenlySpaced); - - Array x, y, z; - - nodalTet->GetNodalPoints(x, y, z); - - SpatialDomains::TetGeomSharedPtr geom = - boost::dynamic_pointer_cast(this->GetGeom(3)); - - // Create basis key for a tetrahedron. - 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)); - LibUtilities::BasisKey C2( - LibUtilities::eOrtho_C, - order + 1, - LibUtilities::PointsKey(order + 1, - LibUtilities::eGaussRadauMAlpha2Beta0)); - - // Create a tet. - LocalRegions::TetExpSharedPtr tet = - MemoryManager::AllocateSharedPtr( - C0, C1, C2, geom); - - // Get coordinate array for tetrahedron. - int nqtot = tet->GetTotPoints(); - Array alloc(6 * nqtot); - Array xi(alloc); - Array yi(alloc + nqtot); - Array zi(alloc + 2 * nqtot); - Array xo(alloc + 3 * nqtot); - Array yo(alloc + 4 * nqtot); - Array zo(alloc + 5 * nqtot); - Array tmp; - - tet->GetCoords(xi, yi, zi); - - for (i = 0; i < 3; ++i) + m_conf.m_order = order; + m_volumeNodes.clear(); + + if (order == 1 || order == 2) + { + m_conf.m_volumeNodes = m_conf.m_faceNodes = false; + return; + } + else if (order == 2) + { + m_conf.m_faceNodes = true; + m_conf.m_volumeNodes = false; + return; + } + else if (order == 3) { - Array coeffs(nodalTet->GetNcoeffs()); - tet->FwdTrans(alloc + i * nqtot, coeffs); - // Apply Vandermonde matrix to project onto nodal space. - nodalTet->ModalToNodal(coeffs, tmp = alloc + (i + 3) * nqtot); + m_conf.m_volumeNodes = false; + return; } - // Now extract points from the co-ordinate arrays into the - // edge/face/volume nodes. First, extract edge-interior nodes. - for (i = 0; i < 6; ++i) + int nPoints = order + 1; + StdRegions::StdExpansionSharedPtr xmap = geom->GetXmap(); + + Array px, py, pz; + LibUtilities::PointsKey pKey(nPoints, pType); + ASSERTL1(pKey.GetPointsDim() == 3, "Points distribution must be 3D"); + LibUtilities::PointsManager()[pKey]->GetPoints(px, py, pz); + + Array > phys(coordDim); + + for (int i = 0; i < coordDim; ++i) { - int pos = 4 + 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(xmap->GetTotPoints()); + xmap->BwdTrans(geom->GetCoeffs(i), phys[i]); } - // Now extract face-interior nodes. - for (i = 0; i < 4; ++i) + const int nTetPts = nPoints * (nPoints + 1) * (nPoints + 2) / 6; + const int nTetIntPts = (nPoints - 4) * (nPoints - 3) * (nPoints - 2) / 6; + m_volumeNodes.resize(nTetIntPts); + + for (int i = 4 + 6*(nPoints-2) + 2*(nPoints-3)*(nPoints-2), cnt = 0; + i < nTetPts; ++i, ++cnt) { - int pos = 4 + 6 * (order - 1) + i * (order - 2) * (order - 1) / 2; - m_face[i]->m_faceNodes.clear(); - for (j = 0; j < (order - 2) * (order - 1) / 2; ++j) + Array xp(3); + xp[0] = px[i]; + xp[1] = py[i]; + xp[2] = pz[i]; + + Array x(3, 0.0); + for (int j = 0; j < coordDim; ++j) { - m_face[i]->m_faceNodes.push_back(NodeSharedPtr( - new Node(0, xo[pos + j], yo[pos + j], zo[pos + j]))); + x[j] = xmap->PhysEvaluate(xp, phys[j]); } - } - // Finally extract volume nodes. - int pos = 4 + 6 * (order - 1) + 4 * (order - 2) * (order - 1) / 2; - for (i = pos; i < (order + 1) * (order + 2) * (order + 3) / 6; ++i) - { - m_volumeNodes.push_back( - NodeSharedPtr(new Node(0, xo[i], yo[i], zo[i]))); + m_volumeNodes[cnt] = boost::shared_ptr( + new Node(id++, x[0], x[1], x[2])); } - m_conf.m_order = order; + m_curveType = pType; m_conf.m_faceNodes = true; m_conf.m_volumeNodes = true; } +/** + * @brief Return the number of nodes defining a tetrahedron. + */ +unsigned int Tetrahedron::GetNumNodes(ElmtConfig pConf) +{ + int n = pConf.m_order; + if (pConf.m_volumeNodes && pConf.m_faceNodes) + return (n + 1) * (n + 2) * (n + 3) / 6; + else if (!pConf.m_volumeNodes && pConf.m_faceNodes) + return 4 * (n + 1) * (n + 2) / 2 - 6 * (n + 1) + 4; + else + return 6 * (n + 1) - 8; +} + struct TetOrient { TetOrient(vector nid, int fid) : nid(nid), fid(fid) diff --git a/library/NekMeshUtils/MeshElements/Tetrahedron.h b/library/NekMeshUtils/MeshElements/Tetrahedron.h index b14437bbe..dbd1ee806 100644 --- a/library/NekMeshUtils/MeshElements/Tetrahedron.h +++ b/library/NekMeshUtils/MeshElements/Tetrahedron.h @@ -77,7 +77,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); diff --git a/utilities/NekMesh/InputModules/InputGmsh.cpp b/utilities/NekMesh/InputModules/InputGmsh.cpp index fd2fbf40e..464b6ce62 100644 --- a/utilities/NekMesh/InputModules/InputGmsh.cpp +++ b/utilities/NekMesh/InputModules/InputGmsh.cpp @@ -37,6 +37,8 @@ #include using namespace std; +#include + #include #include #include @@ -164,6 +166,131 @@ std::vector triTensorNodeOrdering(const std::vector &nodes, int n) return nodeList; } +std::vector tetTensorNodeOrdering(const std::vector &nodes, int n) +{ + std::vector nodeList; + int cnt2; + int nTri = n*(n+1)/2; + int nTet = n*(n+1)*(n+2)/6; + + nodeList.resize(nodes.size()); + nodeList[0] = nodes[0]; + + // Vertices + if (n > 1) + { + nodeList[n - 1] = nodes[1]; + nodeList[nTri - 1] = nodes[2]; + nodeList[nTet - 1] = nodes[3]; + } + + if (n == 1) + { + return nodeList; + } + + // Set up a map that takes (a,b,c) -> m to help us figure out where things + // are inside the tetrahedron. + + typedef boost::tuple 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 tmp; + + for (int k = 0, cnt = 0; k < n; ++k) + { + for (int j = 0; j < n - k; ++j) + { + for (int i = 0; i < n - k - j; ++i) + { + tmp[Mode(i,j,k)] = cnt++; + } + } + } + + for (int i = 1; i < n-1; ++i) + { + int eI = i-1; + nodeList[tmp[Mode(i,0,0)]] = nodes[4 + eI]; + nodeList[tmp[Mode(n-1-i,i,0)]] = nodes[4 + (n-2) + eI]; + nodeList[tmp[Mode(0,n-1-i,0)]] = nodes[4 + 2*(n-2) + eI]; + nodeList[tmp[Mode(0,0,n-1-i)]] = nodes[4 + 3*(n-2) + eI]; + nodeList[tmp[Mode(0,i,n-1-i)]] = nodes[4 + 4*(n-2) + eI]; + nodeList[tmp[Mode(i,0,n-1-i)]] = nodes[4 + 5*(n-2) + eI]; + } + + + // Edges + // if (n > 2) + // { + // int cnt = n, cnt2 = 0; + // for (int i = 1; i < n - 1; ++i) + // { + // nodeList[i] = nodes[4 + i - 1]; + // nodeList[cnt] = nodes[3 + 3 * (n - 2) - i]; + // nodeList[cnt + n - i - 1] = nodes[3 + (n - 2) + i - 1]; + // cnt += n - i; + // } + // } + + // // Faces + + // // Interior (recursion) + + // // Interior (recursion) + // if (n > 3) + // { + // // Reorder interior nodes + // std::vector interior((n - 3) * (n - 2) / 2); + // std::copy( + // nodes.begin() + 3 + 3 * (n - 2), nodes.end(), interior.begin()); + // interior = triTensorNodeOrdering(interior, n - 3); + + // // Copy into full node list + // cnt = n; + // cnt2 = 0; + // for (int j = 1; j < n - 2; ++j) + // { + // for (int i = 0; i < n - j - 2; ++i) + // { + // nodeList[cnt + i + 1] = interior[cnt2 + i]; + // } + // cnt += n - j; + // cnt2 += n - 2 - j; + // } + // } + + return nodeList; +} + /** * @brief Set up InputGmsh object. * @@ -728,6 +855,38 @@ vector InputGmsh::TetReordering(ElmtConfig conf) } } + if (conf.m_volumeNodes == false) + { + return mapping; + } + + const int nInt = (order - 3) * (order - 2) * (order - 1) / 6; + if (nInt <= 0) + { + return mapping; + } + + if (nInt == 1) + { + mapping.push_back(mapping.size()); + return mapping; + } + + int ntot = (order+1)*(order+2)*(order+3)/6; + vector interior; + + for (int i = 4 + 6 * n + 4 * n2; i < ntot; ++i) + { + interior.push_back(i); + } + + if (interior.size() > 0) + { + interior = tetTensorNodeOrdering(interior, order-3); + } + + mapping.insert(mapping.end(), interior.begin(), interior.end()); + return mapping; } @@ -995,7 +1154,7 @@ std::map InputGmsh::GenElmMap() tmp[ 26] = ElmtConfig(eSegment, 3, true, false); tmp[ 27] = ElmtConfig(eSegment, 4, true, false); tmp[ 28] = ElmtConfig(eSegment, 5, true, false); - tmp[ 29] = ElmtConfig(eTetrahedron, 3, true, true); + tmp[ 29] = ElmtConfig(eTetrahedron, 3, true, false); tmp[ 30] = ElmtConfig(eTetrahedron, 4, true, true); tmp[ 31] = ElmtConfig(eTetrahedron, 5, true, true); tmp[ 32] = ElmtConfig(eTetrahedron, 4, true, false); diff --git a/utilities/NekMesh/OutputModules/OutputGmsh.cpp b/utilities/NekMesh/OutputModules/OutputGmsh.cpp index 95779b459..5e36f2248 100644 --- a/utilities/NekMesh/OutputModules/OutputGmsh.cpp +++ b/utilities/NekMesh/OutputModules/OutputGmsh.cpp @@ -34,6 +34,7 @@ //////////////////////////////////////////////////////////////////////////////// #include +#include #include "OutputGmsh.h" #include "../InputModules/InputGmsh.h" @@ -211,6 +212,8 @@ void OutputGmsh::Process() tags.push_back(nodeList[j]->m_id); } + int face_ids[4][3] = {{0, 1, 2}, {0, 1, 3}, {1, 2, 3}, {0, 2, 3}}; + // Process edge-interior points for (int j = 0; j < edgeList.size(); ++j) { @@ -236,9 +239,28 @@ void OutputGmsh::Process() for (int j = 0; j < faceList.size(); ++j) { nodeList = faceList[j]->m_faceNodes; - for (int k = 0; k < nodeList.size(); ++k) + + if (faceList[j]->m_vertexList.size() == 3) { - tags.push_back(nodeList[k]->m_id); + vector 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; + } + + HOTriangle hoTri(faceIds, nodeList); + hoTri.Align(volFaceIds); + + for (int k = 0; k < hoTri.surfVerts.size(); ++k) + { + tags.push_back(hoTri.surfVerts[k]->m_id); + } + } + else + { + // todo } } @@ -248,8 +270,12 @@ void OutputGmsh::Process() tags.push_back(volList[j]->m_id); } + // Construct inverse of input reordering + cout << "ELMT TYPE = " << elmtType << endl; vector reordering = InputGmsh::CreateReordering(elmtType); vector inv(tags.size()); + cout << reordering.size() << " " << tags.size() << endl; + for (int j = 0; j < tags.size(); ++j) { inv[reordering[j]] = j; -- GitLab From 9a19d8465735c88570df63bb367ce981431cef56 Mon Sep 17 00:00:00 2001 From: David Moxey Date: Fri, 29 Jul 2016 12:33:12 +0100 Subject: [PATCH 04/25] Working tets up to order 7 --- utilities/NekMesh/InputModules/InputGmsh.cpp | 125 +++++++++++------- .../NekMesh/OutputModules/OutputGmsh.cpp | 2 +- 2 files changed, 79 insertions(+), 48 deletions(-) diff --git a/utilities/NekMesh/InputModules/InputGmsh.cpp b/utilities/NekMesh/InputModules/InputGmsh.cpp index 464b6ce62..abcb7da66 100644 --- a/utilities/NekMesh/InputModules/InputGmsh.cpp +++ b/utilities/NekMesh/InputModules/InputGmsh.cpp @@ -176,15 +176,17 @@ std::vector tetTensorNodeOrdering(const std::vector &nodes, int n) nodeList.resize(nodes.size()); nodeList[0] = nodes[0]; - // Vertices - if (n > 1) + if (n == 1) { - nodeList[n - 1] = nodes[1]; - nodeList[nTri - 1] = nodes[2]; - nodeList[nTet - 1] = nodes[3]; + return nodeList; } - if (n == 1) + // Vertices + nodeList[n - 1] = nodes[1]; + nodeList[nTri - 1] = nodes[2]; + nodeList[nTet - 1] = nodes[3]; + + if (n == 2) { return nodeList; } @@ -205,7 +207,6 @@ std::vector tetTensorNodeOrdering(const std::vector &nodes, int n) { return false; } - if (a.get<1>() < b.get<1>()) { return true; @@ -214,7 +215,6 @@ std::vector tetTensorNodeOrdering(const std::vector &nodes, int n) { return false; } - if (a.get<2>() < b.get<2>()) { return true; @@ -236,6 +236,7 @@ std::vector tetTensorNodeOrdering(const std::vector &nodes, int n) } } + // Edges for (int i = 1; i < n-1; ++i) { int eI = i-1; @@ -247,46 +248,76 @@ std::vector tetTensorNodeOrdering(const std::vector &nodes, int n) nodeList[tmp[Mode(i,0,n-1-i)]] = nodes[4 + 5*(n-2) + eI]; } + if (n == 3) + { + return nodeList; + } + + // For faces, we use the triTensorNodeOrdering routine to make our lives + // slightly easier. + int nFacePts = (n-3)*(n-2)/2; + cout << "nFacePts = " << nFacePts << endl; - // Edges - // if (n > 2) - // { - // int cnt = n, cnt2 = 0; - // for (int i = 1; i < n - 1; ++i) - // { - // nodeList[i] = nodes[4 + i - 1]; - // nodeList[cnt] = nodes[3 + 3 * (n - 2) - i]; - // nodeList[cnt + n - i - 1] = nodes[3 + (n - 2) + i - 1]; - // cnt += n - i; - // } - // } - - // // Faces - - // // Interior (recursion) - - // // Interior (recursion) - // if (n > 3) - // { - // // Reorder interior nodes - // std::vector interior((n - 3) * (n - 2) / 2); - // std::copy( - // nodes.begin() + 3 + 3 * (n - 2), nodes.end(), interior.begin()); - // interior = triTensorNodeOrdering(interior, n - 3); - - // // Copy into full node list - // cnt = n; - // cnt2 = 0; - // for (int j = 1; j < n - 2; ++j) - // { - // for (int i = 0; i < n - j - 2; ++i) - // { - // nodeList[cnt + i + 1] = interior[cnt2 + i]; - // } - // cnt += n - j; - // cnt2 += n - 2 - j; - // } - // } + // Grab face points and reorder into a tensor-product type format + vector > tmpNodes(4); + int offset = 4 + 6*(n-2); + + for (int i = 0; i < 4; ++i) + { + tmpNodes[i].resize(nFacePts); + for (int j = 0; j < nFacePts; ++j) + { + tmpNodes[i][j] = nodes[offset++]; + } + tmpNodes[i] = triTensorNodeOrdering(tmpNodes[i], n-3); + } + + // Face 0 + for (int j = 1, cnt = 0; j < n-2; ++j) + { + for (int i = 1; i < n-j-1; ++i, ++cnt) + { + cout << "FACE 0 " << i << " " << j << endl; + nodeList[tmp[Mode(i,j,0)]] = tmpNodes[0][cnt]; + } + } + + // Face 1 + for (int j = 1, cnt = 0; j < n-2; ++j) + { + for (int i = 1; i < n-j-1; ++i, ++cnt) + { + cout << "FACE 1 " << i << " " << j << endl; + nodeList[tmp[Mode(i,0,j)]] = tmpNodes[1][cnt]; + } + } + + // Face 2 + for (int j = 1, cnt = 0; j < n-2; ++j) + { + for (int i = 1; i < n-j-1; ++i, ++cnt) + { + cout << "FACE 2 " << i << " " << j << endl; + nodeList[tmp[Mode(i,n-1-i-j,j)]] = tmpNodes[3][cnt]; + } + } + + // Face 3 + for (int j = 1, cnt = 0; j < n-2; ++j) + { + for (int i = 1; i < n-j-1; ++i, ++cnt) + { + cout << "FACE 3 " << i << " " << j << endl; + nodeList[tmp[Mode(0,i,j)]] = tmpNodes[2][cnt]; + } + } + + if (n == 4) + { + return nodeList; + } + + // Finally, recurse on interior volume return nodeList; } diff --git a/utilities/NekMesh/OutputModules/OutputGmsh.cpp b/utilities/NekMesh/OutputModules/OutputGmsh.cpp index 5e36f2248..79de25aa4 100644 --- a/utilities/NekMesh/OutputModules/OutputGmsh.cpp +++ b/utilities/NekMesh/OutputModules/OutputGmsh.cpp @@ -113,7 +113,7 @@ void OutputGmsh::Process() // Convert this mesh into a high-order mesh of uniform order. cout << "Mesh order of " << maxOrder << " detected" << endl; - maxOrder = 6; + maxOrder = 7; m_mesh->MakeOrder(maxOrder, LibUtilities::ePolyEvenlySpaced); // Add edge- and face-interior nodes to vertex set. -- GitLab From e28af50e99ac1121cc22a2b25e8e3f3bbf5f7a0e Mon Sep 17 00:00:00 2001 From: David Moxey Date: Fri, 29 Jul 2016 18:01:25 +0100 Subject: [PATCH 05/25] Tets now work up to order 10 --- utilities/NekMesh/InputModules/InputGmsh.cpp | 79 ++++++++++++------- .../NekMesh/OutputModules/OutputGmsh.cpp | 4 +- 2 files changed, 50 insertions(+), 33 deletions(-) diff --git a/utilities/NekMesh/InputModules/InputGmsh.cpp b/utilities/NekMesh/InputModules/InputGmsh.cpp index abcb7da66..4eff2aa2e 100644 --- a/utilities/NekMesh/InputModules/InputGmsh.cpp +++ b/utilities/NekMesh/InputModules/InputGmsh.cpp @@ -193,7 +193,6 @@ std::vector tetTensorNodeOrdering(const std::vector &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 Mode; struct cmpop { @@ -256,7 +255,6 @@ std::vector tetTensorNodeOrdering(const std::vector &nodes, int n) // For faces, we use the triTensorNodeOrdering routine to make our lives // slightly easier. int nFacePts = (n-3)*(n-2)/2; - cout << "nFacePts = " << nFacePts << endl; // Grab face points and reorder into a tensor-product type format vector > tmpNodes(4); @@ -272,53 +270,74 @@ std::vector tetTensorNodeOrdering(const std::vector &nodes, int n) tmpNodes[i] = triTensorNodeOrdering(tmpNodes[i], n-3); } - // Face 0 - for (int j = 1, cnt = 0; j < n-2; ++j) + if (n > 4) { - for (int i = 1; i < n-j-1; ++i, ++cnt) - { - cout << "FACE 0 " << i << " " << j << endl; - nodeList[tmp[Mode(i,j,0)]] = tmpNodes[0][cnt]; - } + // Now align faces + vector triVertId(3), toAlign(3); + triVertId[0] = 0; + triVertId[1] = 1; + triVertId[2] = 2; + + // Faces 0,2: triangle verts {0,2,1} --> {0,1,2} + HOTriangle hoTri(triVertId, tmpNodes[0]); + toAlign[0] = 0; + toAlign[1] = 2; + toAlign[2] = 1; + + hoTri.Align(toAlign); + tmpNodes[0] = hoTri.surfVerts; + + hoTri.surfVerts = tmpNodes[2]; + hoTri.Align(toAlign); + tmpNodes[2] = hoTri.surfVerts; + + // Face 3: triangle verts {1,2,0} --> {0,1,2} + toAlign[0] = 1; + toAlign[1] = 2; + toAlign[2] = 0; + + hoTri.surfVerts = tmpNodes[3]; + hoTri.Align(toAlign); + tmpNodes[3] = hoTri.surfVerts; } - // Face 1 + // Now apply faces. Note that faces 3 and 2 are swapped between Gmsh and + // Nektar++ order. for (int j = 1, cnt = 0; j < n-2; ++j) { for (int i = 1; i < n-j-1; ++i, ++cnt) { - cout << "FACE 1 " << i << " " << j << endl; - nodeList[tmp[Mode(i,0,j)]] = tmpNodes[1][cnt]; + nodeList[tmp[Mode(i,j,0)]] = tmpNodes[0][cnt]; + nodeList[tmp[Mode(i,0,j)]] = tmpNodes[1][cnt]; + nodeList[tmp[Mode(n-1-i-j,i,j)]] = tmpNodes[3][cnt]; + nodeList[tmp[Mode(0,i,j)]] = tmpNodes[2][cnt]; } } - // Face 2 - for (int j = 1, cnt = 0; j < n-2; ++j) + if (n == 4) { - for (int i = 1; i < n-j-1; ++i, ++cnt) - { - cout << "FACE 2 " << i << " " << j << endl; - nodeList[tmp[Mode(i,n-1-i-j,j)]] = tmpNodes[3][cnt]; - } + return nodeList; } - // Face 3 - for (int j = 1, cnt = 0; j < n-2; ++j) + // Finally, recurse on interior volume + vector intNodes; + for (int i = offset; i < nTet; ++i) { - for (int i = 1; i < n-j-1; ++i, ++cnt) - { - cout << "FACE 3 " << i << " " << j << endl; - nodeList[tmp[Mode(0,i,j)]] = tmpNodes[2][cnt]; - } + intNodes.push_back(nodes[i]); } + intNodes = tetTensorNodeOrdering(intNodes, n-4); - if (n == 4) + for (int k = 1, cnt = 0; k < n - 2; ++k) { - return nodeList; + for (int j = 1; j < n - k - 1; ++j) + { + for (int i = 1; i < n - k - j - 1; ++i) + { + nodeList[tmp[Mode(i,j,k)]] = intNodes[cnt++]; + } + } } - // Finally, recurse on interior volume - return nodeList; } diff --git a/utilities/NekMesh/OutputModules/OutputGmsh.cpp b/utilities/NekMesh/OutputModules/OutputGmsh.cpp index 79de25aa4..b25acaee5 100644 --- a/utilities/NekMesh/OutputModules/OutputGmsh.cpp +++ b/utilities/NekMesh/OutputModules/OutputGmsh.cpp @@ -113,7 +113,7 @@ void OutputGmsh::Process() // Convert this mesh into a high-order mesh of uniform order. cout << "Mesh order of " << maxOrder << " detected" << endl; - maxOrder = 7; + maxOrder = 10; m_mesh->MakeOrder(maxOrder, LibUtilities::ePolyEvenlySpaced); // Add edge- and face-interior nodes to vertex set. @@ -271,10 +271,8 @@ void OutputGmsh::Process() } // Construct inverse of input reordering - cout << "ELMT TYPE = " << elmtType << endl; vector reordering = InputGmsh::CreateReordering(elmtType); vector inv(tags.size()); - cout << reordering.size() << " " << tags.size() << endl; for (int j = 0; j < tags.size(); ++j) { -- GitLab From 1e75d49ae714db6b3a81d6fe081e03c7a4fa5113 Mon Sep 17 00:00:00 2001 From: David Moxey Date: Sat, 30 Jul 2016 09:12:32 +0100 Subject: [PATCH 06/25] Beginnings of a working prism --- utilities/NekMesh/InputModules/InputGmsh.cpp | 78 ++++++++++++++++--- .../NekMesh/OutputModules/OutputGmsh.cpp | 2 +- .../NekMesh/OutputModules/OutputNekpp.cpp | 2 - 3 files changed, 70 insertions(+), 12 deletions(-) diff --git a/utilities/NekMesh/InputModules/InputGmsh.cpp b/utilities/NekMesh/InputModules/InputGmsh.cpp index 4eff2aa2e..32d392275 100644 --- a/utilities/NekMesh/InputModules/InputGmsh.cpp +++ b/utilities/NekMesh/InputModules/InputGmsh.cpp @@ -960,6 +960,7 @@ vector InputGmsh::PrismReordering(ElmtConfig conf) // different; need to mirror in the triangular faces, and then // reorder vertices to make ordering anticlockwise on base quad. static int gmshToNekVerts[6] = {3, 4, 1, 0, 5, 2}; + // 3->0, 4->1, 1->2, 0->3, 5->4, 2->5 for (i = 0; i < 6; ++i) { @@ -1004,18 +1005,77 @@ vector InputGmsh::PrismReordering(ElmtConfig conf) return mapping; } - if (order > 2) + int nTriInt = n * (n - 1) / 2; + int nQuadInt = n * n; + + // Gmsh faces: + // 0 = {0,2,1}, 1 = {3,4,5}, 2 = {0,1,4,3}, 3 = {0,3,5,2}, 4 = {1,2,5,4} + // Apply gmshToNekVerts map: + // 0 = {3,5,2}, 1 = {0,1,4}, 2 = {3,2,1,0}, 3 = {3,0,4,5}, 4 = {2,5,4,1} + // Nektar++ faces: + // 0 = {0,1,2,3}, 1 = {0,1,4}, 2 = {1,2,5,4}, 3 = {3,2,5}, 4 = {0,3,5,4} + // {3,2,1,0} {0,1,4} {2,5,4,1} {3,5,2} {3,0,4,5} + // This gives gmsh -> Nektar++ faces: + static int gmshToNekFace[5] = {2, 1, 4, 0, 3}; + + // Face 0: StdRegions::eDir1FwdDir1_Dir2BwdDir2 + // Face 2: StdRegions::eDir1BwdDir2_Dir2FwdDir1 + // Face 4: StdRegions::eDir1BwdDir1_Dir2FwdDir2 + + // Face offsets + vector offsets(5); + offset = 6 + 9 * n; + offsets[0] = offset; + offsets[1] = offset + nTriInt; + offsets[2] = offset + 2 * nTriInt; + offsets[3] = offset + 2 * nTriInt + nQuadInt; + offsets[4] = offset + 2 * nTriInt + 2 * nQuadInt; + + mapping.resize(6 + 9 * n + 3 * nQuadInt + 2 * nTriInt); + + offset = 6 + 9 * n; + + // Face 0 + int offset2 = offsets[gmshToNekFace[0]]; + int j; + for (i = 0; i < n; ++i) { - cerr << "Gmsh prisms of order > 2 with face curvature " - << "not supported in NekMesh (or indeed Gmsh at" - << "time of writing)." << endl; - abort(); + for (j = 0; j < n; ++j) + { + mapping[offset + i * n + j] = offset2 + (n - i - 1) * n + j; + } + } + offset += nQuadInt; + + // Face 1 + offset2 = offsets[gmshToNekFace[1]]; + for (i = 0; i < nTriInt; ++i) + { + mapping[offset++] = offset2++; + } + + // Face 1 + offset2 = offsets[gmshToNekFace[2]]; + for (i = 0; i < nQuadInt; ++i) + { + mapping[offset++] = offset2++; + cout << offset-1 << " " << offset2-1 << endl; } - mapping.resize(18); - mapping[15] = 15; - mapping[16] = 17; - mapping[17] = 16; + // Face 1 + offset2 = offsets[gmshToNekFace[3]]; + for (i = 0; i < nTriInt; ++i) + { + mapping[offset++] = offset2++; + } + + // Face 1 + offset2 = offsets[gmshToNekFace[4]]; + for (i = 0; i < nQuadInt; ++i) + { + mapping[offset++] = offset2++; + cout << offset-1 << " " << offset2-1 << endl; + } return mapping; } diff --git a/utilities/NekMesh/OutputModules/OutputGmsh.cpp b/utilities/NekMesh/OutputModules/OutputGmsh.cpp index b25acaee5..f3eb69b5a 100644 --- a/utilities/NekMesh/OutputModules/OutputGmsh.cpp +++ b/utilities/NekMesh/OutputModules/OutputGmsh.cpp @@ -113,7 +113,7 @@ void OutputGmsh::Process() // Convert this mesh into a high-order mesh of uniform order. cout << "Mesh order of " << maxOrder << " detected" << endl; - maxOrder = 10; + maxOrder = 3; m_mesh->MakeOrder(maxOrder, LibUtilities::ePolyEvenlySpaced); // Add edge- and face-interior nodes to vertex set. diff --git a/utilities/NekMesh/OutputModules/OutputNekpp.cpp b/utilities/NekMesh/OutputModules/OutputNekpp.cpp index 4fb448049..56a86d26f 100644 --- a/utilities/NekMesh/OutputModules/OutputNekpp.cpp +++ b/utilities/NekMesh/OutputModules/OutputNekpp.cpp @@ -98,8 +98,6 @@ void OutputNekpp::Process() cout << "OutputNekpp: Writing file..." << endl; } - m_mesh->MakeOrder(5, LibUtilities::ePolyEvenlySpaced); - TiXmlDocument doc; TiXmlDeclaration *decl = new TiXmlDeclaration("1.0", "utf-8", ""); doc.LinkEndChild(decl); -- GitLab From 2adf60d58c85e9832b0dff5f422795f4d3b2ff82 Mon Sep 17 00:00:00 2001 From: David Moxey Date: Fri, 5 Aug 2016 10:34:10 +0100 Subject: [PATCH 07/25] Working hex up to order 4 --- .../NekMeshUtils/MeshElements/Hexahedron.cpp | 84 +++++++++ .../NekMeshUtils/MeshElements/Hexahedron.h | 8 + library/NekMeshUtils/MeshElements/Mesh.cpp | 1 + .../MeshElements/Quadrilateral.cpp | 17 +- .../NekMeshUtils/MeshElements/Tetrahedron.cpp | 11 +- utilities/NekMesh/InputModules/InputGmsh.cpp | 162 ++++++++++++++++-- .../NekMesh/OutputModules/OutputGmsh.cpp | 9 +- 7 files changed, 270 insertions(+), 22 deletions(-) diff --git a/library/NekMeshUtils/MeshElements/Hexahedron.cpp b/library/NekMeshUtils/MeshElements/Hexahedron.cpp index d6cf48cde..cde1e7406 100644 --- a/library/NekMeshUtils/MeshElements/Hexahedron.cpp +++ b/library/NekMeshUtils/MeshElements/Hexahedron.cpp @@ -176,6 +176,90 @@ SpatialDomains::GeometrySharedPtr Hexahedron::GetGeom(int coordDim) return ret; } +StdRegions::Orientation Hexahedron::GetEdgeOrient( + int edgeId, EdgeSharedPtr edge) +{ + static int edgeVerts[12][2] = { {0,1}, {1,2}, {2,3}, {3,0}, {0,4}, {1,5}, + {2,6}, {3,7}, {4,5}, {5,6}, {6,7}, {7,4} }; + + if (edge->m_n1 == m_vertex[edgeVerts[edgeId][0]]) + { + return StdRegions::eForwards; + } + else if (edge->m_n1 == m_vertex[edgeVerts[edgeId][1]]) + { + return StdRegions::eBackwards; + } + else + { + ASSERTL1(false, "Edge is not connected to this hexahedron."); + } + + return StdRegions::eNoOrientation; +} + +void Hexahedron::MakeOrder(int order, + SpatialDomains::GeometrySharedPtr geom, + LibUtilities::PointsType pType, + int coordDim, + int &id) +{ + m_conf.m_order = order; + m_volumeNodes.clear(); + + if (order == 1) + { + m_conf.m_volumeNodes = m_conf.m_faceNodes = false; + return; + } + + m_conf.m_faceNodes = true; + m_conf.m_volumeNodes = true; + m_curveType = pType; + + int nPoints = order + 1; + StdRegions::StdExpansionSharedPtr xmap = geom->GetXmap(); + + Array px; + LibUtilities::PointsKey pKey(nPoints, pType); + ASSERTL1(pKey.GetPointsDim() == 1, "Points distribution must be 1D"); + LibUtilities::PointsManager()[pKey]->GetPoints(px); + + Array > phys(coordDim); + + for (int i = 0; i < coordDim; ++i) + { + phys[i] = Array(xmap->GetTotPoints()); + xmap->BwdTrans(geom->GetCoeffs(i), phys[i]); + } + + int nHexIntPts = (nPoints - 2) * (nPoints - 2) * (nPoints - 2); + m_volumeNodes.resize(nHexIntPts); + + for (int i = 1, cnt = 0; i < nPoints-1; ++i) + { + for (int j = 1; j < nPoints-1; ++j) + { + for (int k = 1; k < nPoints-1; ++k, ++cnt) + { + Array xp(3); + xp[0] = px[k]; + xp[1] = px[j]; + xp[2] = px[i]; + + Array x(3, 0.0); + for (int k = 0; k < coordDim; ++k) + { + x[k] = xmap->PhysEvaluate(xp, phys[k]); + } + + m_volumeNodes[cnt] = boost::shared_ptr( + new Node(id++, x[0], x[1], x[2])); + } + } + } +} + /** * @brief Return the number of nodes defining a hexahedron. */ diff --git a/library/NekMeshUtils/MeshElements/Hexahedron.h b/library/NekMeshUtils/MeshElements/Hexahedron.h index 8f9aebade..43e3f421c 100644 --- a/library/NekMeshUtils/MeshElements/Hexahedron.h +++ b/library/NekMeshUtils/MeshElements/Hexahedron.h @@ -77,6 +77,14 @@ public: NEKMESHUTILS_EXPORT virtual SpatialDomains::GeometrySharedPtr GetGeom( int coordDim); + 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); }; diff --git a/library/NekMeshUtils/MeshElements/Mesh.cpp b/library/NekMeshUtils/MeshElements/Mesh.cpp index 66e3dd963..d24768667 100644 --- a/library/NekMeshUtils/MeshElements/Mesh.cpp +++ b/library/NekMeshUtils/MeshElements/Mesh.cpp @@ -104,6 +104,7 @@ void Mesh::MakeOrder(int order, pTypes[LibUtilities::eTriangle] = LibUtilities::eNodalTriEvenlySpaced; pTypes[LibUtilities::eQuadrilateral] = LibUtilities::ePolyEvenlySpaced; pTypes[LibUtilities::eTetrahedron] = LibUtilities::eNodalTetEvenlySpaced; + pTypes[LibUtilities::eHexahedron] = LibUtilities::ePolyEvenlySpaced; } for(eit = m_edgeSet.begin(); eit != m_edgeSet.end(); eit++) diff --git a/library/NekMeshUtils/MeshElements/Quadrilateral.cpp b/library/NekMeshUtils/MeshElements/Quadrilateral.cpp index 1b1d5b4ac..84f74963e 100644 --- a/library/NekMeshUtils/MeshElements/Quadrilateral.cpp +++ b/library/NekMeshUtils/MeshElements/Quadrilateral.cpp @@ -140,6 +140,18 @@ void Quadrilateral::MakeOrder(int order, int coordDim, int &id) { + m_conf.m_order = order; + + if (order == 1) + { + m_conf.m_volumeNodes = m_conf.m_faceNodes = false; + return; + } + + m_conf.m_faceNodes = true; + m_conf.m_volumeNodes = false; + m_curveType = pType; + // Quadrilaterals of order < 2 have no interior volume points. if (order < 2) { @@ -184,11 +196,6 @@ void Quadrilateral::MakeOrder(int order, 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 = false; } SpatialDomains::GeometrySharedPtr Quadrilateral::GetGeom(int coordDim) diff --git a/library/NekMeshUtils/MeshElements/Tetrahedron.cpp b/library/NekMeshUtils/MeshElements/Tetrahedron.cpp index bf44253c5..bdb67ab2c 100644 --- a/library/NekMeshUtils/MeshElements/Tetrahedron.cpp +++ b/library/NekMeshUtils/MeshElements/Tetrahedron.cpp @@ -233,8 +233,7 @@ SpatialDomains::GeometrySharedPtr Tetrahedron::GetGeom(int coordDim) StdRegions::Orientation Tetrahedron::GetEdgeOrient( int edgeId, EdgeSharedPtr edge) { - int locVert = edgeId; - int edgeVerts[6][2] = { {0,1}, {1,2}, {0,2}, {0,3}, {1,3}, {2,3} }; + static int edgeVerts[6][2] = { {0,1}, {1,2}, {0,2}, {0,3}, {1,3}, {2,3} }; if (edge->m_n1 == m_vertex[edgeVerts[edgeId][0]]) { @@ -278,6 +277,10 @@ void Tetrahedron::MakeOrder(int order, return; } + m_curveType = pType; + m_conf.m_faceNodes = true; + m_conf.m_volumeNodes = true; + int nPoints = order + 1; StdRegions::StdExpansionSharedPtr xmap = geom->GetXmap(); @@ -315,10 +318,6 @@ void Tetrahedron::MakeOrder(int order, m_volumeNodes[cnt] = boost::shared_ptr( new Node(id++, x[0], x[1], x[2])); } - - m_curveType = pType; - m_conf.m_faceNodes = true; - m_conf.m_volumeNodes = true; } /** diff --git a/utilities/NekMesh/InputModules/InputGmsh.cpp b/utilities/NekMesh/InputModules/InputGmsh.cpp index 32d392275..4549e6041 100644 --- a/utilities/NekMesh/InputModules/InputGmsh.cpp +++ b/utilities/NekMesh/InputModules/InputGmsh.cpp @@ -341,6 +341,147 @@ std::vector tetTensorNodeOrdering(const std::vector &nodes, int n) return nodeList; } +std::vector hexTensorNodeOrdering(const std::vector &nodes, int n) +{ + std::vector nodeList; + int cnt2; + int nEdge = n-2; + int nFace = nEdge * nEdge; + int nHex = n * n * n; + + nodeList.resize(nodes.size()); + nodeList[0] = nodes[0]; + + if (n == 1) + { + return nodeList; + } + + // Vertices: same order as Nektar++ + nodeList[n - 1] = nodes[1]; + nodeList[n*n -1] = nodes[2]; + nodeList[n*(n-1)] = nodes[3]; + nodeList[n*n*(n-1)] = nodes[4]; + nodeList[n - 1 + n*n*(n-1)] = nodes[5]; + nodeList[n*n -1 + n*n*(n-1)] = nodes[6]; + nodeList[n*(n-1) + n*n*(n-1)] = nodes[7]; + + if (n == 2) + { + return nodeList; + } + + // static int hexEdges[12][2] = + // { { 0, 1 }, { n, n*2 } }; + // static int gmshToNekEdge[12] = {0, -3, 4, 1, 5, 2, 6, 7, 8, -11, 9, 10}; + + // // Edges + // for (int i = 1; i < n-1; ++i) + // { + // int eI = i-1; + // nodeList[tmp[Mode(i,0,0)]] = nodes[4 + eI]; + // nodeList[tmp[Mode(n-1-i,i,0)]] = nodes[4 + (n-2) + eI]; + // nodeList[tmp[Mode(0,n-1-i,0)]] = nodes[4 + 2*(n-2) + eI]; + // nodeList[tmp[Mode(0,0,n-1-i)]] = nodes[4 + 3*(n-2) + eI]; + // nodeList[tmp[Mode(0,i,n-1-i)]] = nodes[4 + 4*(n-2) + eI]; + // nodeList[tmp[Mode(i,0,n-1-i)]] = nodes[4 + 5*(n-2) + eI]; + // } + + // if (n == 3) + // { + // return nodeList; + // } + + // // For faces, we use the triTensorNodeOrdering routine to make our lives + // // slightly easier. + // int nFacePts = (n-3)*(n-2)/2; + + // // Grab face points and reorder into a tensor-product type format + // vector > tmpNodes(4); + // int offset = 4 + 6*(n-2); + + // for (int i = 0; i < 4; ++i) + // { + // tmpNodes[i].resize(nFacePts); + // for (int j = 0; j < nFacePts; ++j) + // { + // tmpNodes[i][j] = nodes[offset++]; + // } + // tmpNodes[i] = triTensorNodeOrdering(tmpNodes[i], n-3); + // } + + // if (n > 4) + // { + // // Now align faces + // vector triVertId(3), toAlign(3); + // triVertId[0] = 0; + // triVertId[1] = 1; + // triVertId[2] = 2; + + // // Faces 0,2: triangle verts {0,2,1} --> {0,1,2} + // HOTriangle hoTri(triVertId, tmpNodes[0]); + // toAlign[0] = 0; + // toAlign[1] = 2; + // toAlign[2] = 1; + + // hoTri.Align(toAlign); + // tmpNodes[0] = hoTri.surfVerts; + + // hoTri.surfVerts = tmpNodes[2]; + // hoTri.Align(toAlign); + // tmpNodes[2] = hoTri.surfVerts; + + // // Face 3: triangle verts {1,2,0} --> {0,1,2} + // toAlign[0] = 1; + // toAlign[1] = 2; + // toAlign[2] = 0; + + // hoTri.surfVerts = tmpNodes[3]; + // hoTri.Align(toAlign); + // tmpNodes[3] = hoTri.surfVerts; + // } + + // // Now apply faces. Note that faces 3 and 2 are swapped between Gmsh and + // // Nektar++ order. + // for (int j = 1, cnt = 0; j < n-2; ++j) + // { + // for (int i = 1; i < n-j-1; ++i, ++cnt) + // { + // nodeList[tmp[Mode(i,j,0)]] = tmpNodes[0][cnt]; + // nodeList[tmp[Mode(i,0,j)]] = tmpNodes[1][cnt]; + // nodeList[tmp[Mode(n-1-i-j,i,j)]] = tmpNodes[3][cnt]; + // nodeList[tmp[Mode(0,i,j)]] = tmpNodes[2][cnt]; + // } + // } + + // if (n == 4) + // { + // return nodeList; + // } + + // // Finally, recurse on interior volume + // vector intNodes; + // for (int i = offset; i < nTet; ++i) + // { + // intNodes.push_back(nodes[i]); + // } + // intNodes = tetTensorNodeOrdering(intNodes, n-4); + + // for (int k = 1, cnt = 0; k < n - 2; ++k) + // { + // for (int j = 1; j < n - k - 1; ++j) + // { + // for (int i = 1; i < n - k - j - 1; ++i) + // { + // nodeList[tmp[Mode(i,j,k)]] = intNodes[cnt++]; + // } + // } + // } + + return nodeList; +} + + /** * @brief Set up InputGmsh object. * @@ -1209,14 +1350,15 @@ vector InputGmsh::HexReordering(ElmtConfig conf) } const int totPoints = (order + 1) * (order + 1) * (order + 1); - mapping.resize(totPoints); - - // TODO: Fix ordering of volume nodes. + vector interior; for (i = 8 + 12 * n + 6 * n2; i < totPoints; ++i) { - mapping[i] = i; + interior.push_back(i); } + interior = hexTensorNodeOrdering(interior, order - 1); + mapping.insert(mapping.end(), interior.begin(), interior.end()); + return mapping; } @@ -1236,12 +1378,12 @@ std::map InputGmsh::GenElmMap() std::map tmp; // Elmt type, order, face, volume - tmp[ 1] = ElmtConfig(eSegment, 1, true, true); - tmp[ 2] = ElmtConfig(eTriangle, 1, true, true); - tmp[ 3] = ElmtConfig(eQuadrilateral, 1, true, true); - tmp[ 4] = ElmtConfig(eTetrahedron, 1, true, true); - tmp[ 5] = ElmtConfig(eHexahedron, 1, true, true); - tmp[ 6] = ElmtConfig(ePrism, 1, true, true); + tmp[ 1] = ElmtConfig(eSegment, 1, false, false); + tmp[ 2] = ElmtConfig(eTriangle, 1, false, false); + tmp[ 3] = ElmtConfig(eQuadrilateral, 1, false, false); + tmp[ 4] = ElmtConfig(eTetrahedron, 1, false, false); + tmp[ 5] = ElmtConfig(eHexahedron, 1, false, false); + tmp[ 6] = ElmtConfig(ePrism, 1, false, false); tmp[ 7] = ElmtConfig(ePyramid, 1, true, true); tmp[ 8] = ElmtConfig(eSegment, 2, true, true); tmp[ 9] = ElmtConfig(eTriangle, 2, true, true); diff --git a/utilities/NekMesh/OutputModules/OutputGmsh.cpp b/utilities/NekMesh/OutputModules/OutputGmsh.cpp index f3eb69b5a..37a74069a 100644 --- a/utilities/NekMesh/OutputModules/OutputGmsh.cpp +++ b/utilities/NekMesh/OutputModules/OutputGmsh.cpp @@ -180,6 +180,7 @@ void OutputGmsh::Process() // First output element ID and type. int elmtType = elmMap[e->GetConf()]; m_mshFile << id << " " << elmtType << " "; + cout << elmtType << endl; // Write out number of element tags and then the tags // themselves. @@ -260,7 +261,11 @@ void OutputGmsh::Process() } else { - // todo + // NEEDS ALIGNMENT + for (int k = 0; k < nodeList.size(); ++k) + { + tags.push_back(nodeList[k]->m_id); + } } } @@ -274,6 +279,8 @@ void OutputGmsh::Process() vector reordering = InputGmsh::CreateReordering(elmtType); vector inv(tags.size()); + cout << tags.size() << endl; + for (int j = 0; j < tags.size(); ++j) { inv[reordering[j]] = j; -- GitLab From 30ec5d1f67dcda014c0e4ce99572f403efbecd61 Mon Sep 17 00:00:00 2001 From: David Moxey Date: Fri, 5 Aug 2016 16:03:11 +0100 Subject: [PATCH 08/25] Actually working hex up to order 4 --- utilities/NekMesh/InputModules/InputGmsh.cpp | 234 ++++++++++-------- .../NekMesh/OutputModules/OutputGmsh.cpp | 181 +++++++++++++- 2 files changed, 296 insertions(+), 119 deletions(-) diff --git a/utilities/NekMesh/InputModules/InputGmsh.cpp b/utilities/NekMesh/InputModules/InputGmsh.cpp index 4549e6041..dc21063e0 100644 --- a/utilities/NekMesh/InputModules/InputGmsh.cpp +++ b/utilities/NekMesh/InputModules/InputGmsh.cpp @@ -343,11 +343,8 @@ std::vector tetTensorNodeOrdering(const std::vector &nodes, int n) std::vector hexTensorNodeOrdering(const std::vector &nodes, int n) { + int i, j, k; std::vector nodeList; - int cnt2; - int nEdge = n-2; - int nFace = nEdge * nEdge; - int nHex = n * n * n; nodeList.resize(nodes.size()); nodeList[0] = nodes[0]; @@ -371,112 +368,129 @@ std::vector hexTensorNodeOrdering(const std::vector &nodes, int n) return nodeList; } - // static int hexEdges[12][2] = - // { { 0, 1 }, { n, n*2 } }; - // static int gmshToNekEdge[12] = {0, -3, 4, 1, 5, 2, 6, 7, 8, -11, 9, 10}; - - // // Edges - // for (int i = 1; i < n-1; ++i) - // { - // int eI = i-1; - // nodeList[tmp[Mode(i,0,0)]] = nodes[4 + eI]; - // nodeList[tmp[Mode(n-1-i,i,0)]] = nodes[4 + (n-2) + eI]; - // nodeList[tmp[Mode(0,n-1-i,0)]] = nodes[4 + 2*(n-2) + eI]; - // nodeList[tmp[Mode(0,0,n-1-i)]] = nodes[4 + 3*(n-2) + eI]; - // nodeList[tmp[Mode(0,i,n-1-i)]] = nodes[4 + 4*(n-2) + eI]; - // nodeList[tmp[Mode(i,0,n-1-i)]] = nodes[4 + 5*(n-2) + eI]; - // } - - // if (n == 3) - // { - // return nodeList; - // } - - // // For faces, we use the triTensorNodeOrdering routine to make our lives - // // slightly easier. - // int nFacePts = (n-3)*(n-2)/2; - - // // Grab face points and reorder into a tensor-product type format - // vector > tmpNodes(4); - // int offset = 4 + 6*(n-2); - - // for (int i = 0; i < 4; ++i) - // { - // tmpNodes[i].resize(nFacePts); - // for (int j = 0; j < nFacePts; ++j) - // { - // tmpNodes[i][j] = nodes[offset++]; - // } - // tmpNodes[i] = triTensorNodeOrdering(tmpNodes[i], n-3); - // } - - // if (n > 4) - // { - // // Now align faces - // vector triVertId(3), toAlign(3); - // triVertId[0] = 0; - // triVertId[1] = 1; - // triVertId[2] = 2; - - // // Faces 0,2: triangle verts {0,2,1} --> {0,1,2} - // HOTriangle hoTri(triVertId, tmpNodes[0]); - // toAlign[0] = 0; - // toAlign[1] = 2; - // toAlign[2] = 1; - - // hoTri.Align(toAlign); - // tmpNodes[0] = hoTri.surfVerts; - - // hoTri.surfVerts = tmpNodes[2]; - // hoTri.Align(toAlign); - // tmpNodes[2] = hoTri.surfVerts; - - // // Face 3: triangle verts {1,2,0} --> {0,1,2} - // toAlign[0] = 1; - // toAlign[1] = 2; - // toAlign[2] = 0; - - // hoTri.surfVerts = tmpNodes[3]; - // hoTri.Align(toAlign); - // tmpNodes[3] = hoTri.surfVerts; - // } - - // // Now apply faces. Note that faces 3 and 2 are swapped between Gmsh and - // // Nektar++ order. - // for (int j = 1, cnt = 0; j < n-2; ++j) - // { - // for (int i = 1; i < n-j-1; ++i, ++cnt) - // { - // nodeList[tmp[Mode(i,j,0)]] = tmpNodes[0][cnt]; - // nodeList[tmp[Mode(i,0,j)]] = tmpNodes[1][cnt]; - // nodeList[tmp[Mode(n-1-i-j,i,j)]] = tmpNodes[3][cnt]; - // nodeList[tmp[Mode(0,i,j)]] = tmpNodes[2][cnt]; - // } - // } - - // if (n == 4) - // { - // return nodeList; - // } - - // // Finally, recurse on interior volume - // vector intNodes; - // for (int i = offset; i < nTet; ++i) - // { - // intNodes.push_back(nodes[i]); - // } - // intNodes = tetTensorNodeOrdering(intNodes, n-4); - - // for (int k = 1, cnt = 0; k < n - 2; ++k) - // { - // for (int j = 1; j < n - k - 1; ++j) - // { - // for (int i = 1; i < n - k - j - 1; ++i) - // { - // nodeList[tmp[Mode(i,j,k)]] = intNodes[cnt++]; - // } - // } - // } + static int hexEdges[12][2] = { + { 0, 1 }, { n-1, n }, { n*n-1, -1 }, { n*(n-1), -n }, + { 0, n*n }, { n-1, n*n }, { n*n - 1, n*n }, { n*(n-1), n*n }, + { n*n*(n-1), 1 }, { n*n*(n-1) + n-1, n }, { n*n*n-1, -1 }, + { n*n*(n-1) + n*(n-1), -n } + }; + static int hexFaces[6][3] = { + { 0, 1, n }, { 0, 1, n*n }, { n-1, n, n*n }, + { n*(n-1), 1, n*n }, { 0, n, n*n }, { n*n*(n-1), 1, n } + }; + static int gmshToNekEdge[12] = {0, -3, 4, 1, 5, 2, 6, 7, 8, -11, 9, 10}; + + // Edges + int offset = 8; + for (int i = 0; i < 12; ++i) + { + int e = abs(gmshToNekEdge[i]); + + if (gmshToNekEdge[i] >= 0) + { + for (int j = 1; j < n-1; ++j) + { + nodeList[hexEdges[e][0] + j*hexEdges[e][1]] = nodes[offset++]; + } + } + else + { + for (int j = 1; j < n-1; ++j) + { + nodeList[hexEdges[e][0] + (n-j-1)*hexEdges[e][1]] = nodes[offset++]; + } + } + } + + // Faces + static int gmsh2NekFace[6] = {0, 1, 4, 2, 3, 5}; + + // Map which defines orientation between Gmsh and Nektar++ faces. + StdRegions::Orientation faceOrient[6] = { + StdRegions::eDir1FwdDir2_Dir2FwdDir1, + StdRegions::eDir1FwdDir1_Dir2FwdDir2, + StdRegions::eDir1FwdDir2_Dir2FwdDir1, + StdRegions::eDir1FwdDir1_Dir2FwdDir2, + StdRegions::eDir1BwdDir1_Dir2FwdDir2, + StdRegions::eDir1FwdDir1_Dir2FwdDir2}; + + for (i = 0; i < 6; ++i) + { + int n2 = (n-2)*(n-2); + int face = gmsh2NekFace[i]; + offset = 8 + 12 * (n-2) + i * n2; + + // Create a list of interior face nodes for this face only. + vector faceNodes(n2); + for (j = 0; j < n2; ++j) + { + faceNodes[j] = nodes[offset + j]; + } + + // Now get the reordering of this face, which puts Gmsh + // recursive ordering into Nektar++ row-by-row order. + faceNodes = quadTensorNodeOrdering(faceNodes, n-2); + vector tmp(n2); + + // Finally reorient the face according to the geometry + // differences. + if (faceOrient[i] == StdRegions::eDir1FwdDir1_Dir2FwdDir2) + { + // Orientation is the same, just copy. + tmp = faceNodes; + } + else if (faceOrient[i] == StdRegions::eDir1FwdDir2_Dir2FwdDir1) + { + // Tranposed faces + for (j = 0; j < n-2; ++j) + { + for (k = 0; k < n-2; ++k) + { + tmp[j * (n-2) + k] = faceNodes[k * (n-2) + j]; + } + } + } + else if (faceOrient[i] == StdRegions::eDir1BwdDir1_Dir2FwdDir2) + { + for (j = 0; j < n-2; ++j) + { + for (k = 0; k < n-2; ++k) + { + tmp[j * (n-2) + k] = faceNodes[j * (n-2) + (n - k - 3)]; + } + } + } + + // Now put this into the right place in the output array + for (k = 1; k < n-1; ++k) + { + for (j = 1; j < n-1; ++j) + { + nodeList[hexFaces[face][0] + j*hexFaces[face][1] + k*hexFaces[face][2]] + = faceNodes[(j-1)*(n-2) + j - 1]; + } + } + } + + // Finally, recurse on interior volume + vector intNodes; + for (int i = 8 + 12 * (n-2) + 6 * (n-2) * (n-2); i < n*n*n; ++i) + { + intNodes.push_back(nodes[i]); + } + intNodes = hexTensorNodeOrdering(intNodes, n-2); + + for (int k = 1, cnt = 0; k < n - 1; ++k) + { + for (int j = 1; j < n - 1; ++j) + { + for (int i = 1; i < n - 1; ++i) + { + cout << "LOLINT " << i + j * n + k * n * n << " " << intNodes[cnt] << endl; + nodeList[i + j * n + k * n * n] = intNodes[cnt++]; + } + } + } return nodeList; } diff --git a/utilities/NekMesh/OutputModules/OutputGmsh.cpp b/utilities/NekMesh/OutputModules/OutputGmsh.cpp index 37a74069a..0ad0baf86 100644 --- a/utilities/NekMesh/OutputModules/OutputGmsh.cpp +++ b/utilities/NekMesh/OutputModules/OutputGmsh.cpp @@ -52,6 +52,157 @@ ModuleKey OutputGmsh::className = OutputGmsh::create, "Writes Gmsh msh file."); +/** + * @brief A lightweight struct for dealing with high-order quadrilateral + * alignment. + */ +template struct HOQuadrilateral +{ + HOQuadrilateral(std::vector pVertId, std::vector pSurfVerts) + : vertId(pVertId), surfVerts(pSurfVerts) + { + } + + HOQuadrilateral(std::vector pVertId) : vertId(pVertId) + { + } + + /// The quadrilateral vertex IDs + std::vector vertId; + + /// The quadrilateral surface vertices -- templated so that this can either + /// be nodes or IDs. + std::vector surfVerts; + + void ReverseX() + { + int np = (int)(sqrt(surfVerts.size()) + 0.5); + for (int i = 0; i < np; ++i) + { + for (int j = 0; j < np/2; ++j) + { + swap(surfVerts[i*np + j], surfVerts[i*np + np-j-1]); + } + } + } + + void ReverseY() + { + int np = (int)(sqrt(surfVerts.size()) + 0.5); + // Reverse y direction + for (int j = 0; j < np; ++j) + { + for (int i = 0; i < np/2; ++i) + { + swap(surfVerts[i*np + j], surfVerts[(np-i-1)*np + j]); + } + } + } + + void Transpose() + { + int np = (int)(sqrt(surfVerts.size()) + 0.5); + vector tmp(surfVerts.size()); + + for (int i = 0; i < np; ++i) + { + for (int j = 0; j < np; ++j) + { + tmp[i*np+j] = surfVerts[j*np+i]; + } + } + + surfVerts = tmp; + } + + /** + * @brief Align this surface to a given vertex ID. + */ + void Align(std::vector vertId) + { + int vmap[4] = {-1, -1, -1, -1}; + + // Determine which vertices map to vertId + for (int i = 0; i < 4; ++i) + { + for (int j = 0; j < 4; ++j) + { + if (this->vertId[j] == vertId[i]) + { + vmap[i] = j; + break; + } + } + + ASSERTL1(vmap[i] != -1, + "Could not determine mapping between vertex IDs"); + } + + StdRegions::Orientation orient = StdRegions::eNoOrientation; + + if (vmap[1] == (vmap[0]+1) % 4) + { + switch (vmap[0]) + { + case 0: + orient = StdRegions::eDir1FwdDir1_Dir2FwdDir2; + break; + case 1: + orient = StdRegions::eDir1BwdDir2_Dir2FwdDir1; + break; + case 2: + orient = StdRegions::eDir1BwdDir1_Dir2BwdDir2; + break; + case 3: + orient = StdRegions::eDir1FwdDir2_Dir2BwdDir1; + break; + } + } + else + { + switch (vmap[0]) + { + case 0: + orient = StdRegions::eDir1FwdDir2_Dir2FwdDir1; + break; + case 1: + orient = StdRegions::eDir1BwdDir1_Dir2FwdDir2; + break; + case 2: + orient = StdRegions::eDir1BwdDir2_Dir2BwdDir1; + break; + case 3: + orient = StdRegions::eDir1FwdDir1_Dir2BwdDir2; + break; + } + } + + if (orient == StdRegions::eDir1FwdDir2_Dir2FwdDir1 || + orient == StdRegions::eDir1BwdDir2_Dir2FwdDir1 || + orient == StdRegions::eDir1FwdDir2_Dir2BwdDir1 || + orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1) + { + Transpose(); + } + + if (orient == StdRegions::eDir1BwdDir1_Dir2FwdDir2 || + orient == StdRegions::eDir1BwdDir1_Dir2BwdDir2 || + orient == StdRegions::eDir1FwdDir2_Dir2BwdDir1 || + orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1) + { + ReverseX(); + } + + if (orient == StdRegions::eDir1FwdDir1_Dir2BwdDir2 || + orient == StdRegions::eDir1BwdDir1_Dir2BwdDir2 || + orient == StdRegions::eDir1BwdDir2_Dir2FwdDir1 || + orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1) + { + ReverseY(); + } + } +}; + OutputGmsh::OutputGmsh(MeshSharedPtr m) : OutputModule(m) { map igelmap = InputGmsh::GenElmMap(); @@ -113,7 +264,7 @@ void OutputGmsh::Process() // Convert this mesh into a high-order mesh of uniform order. cout << "Mesh order of " << maxOrder << " detected" << endl; - maxOrder = 3; + maxOrder = 4; m_mesh->MakeOrder(maxOrder, LibUtilities::ePolyEvenlySpaced); // Add edge- and face-interior nodes to vertex set. @@ -180,7 +331,6 @@ void OutputGmsh::Process() // First output element ID and type. int elmtType = elmMap[e->GetConf()]; m_mshFile << id << " " << elmtType << " "; - cout << elmtType << endl; // Write out number of element tags and then the tags // themselves. @@ -213,8 +363,6 @@ void OutputGmsh::Process() tags.push_back(nodeList[j]->m_id); } - int face_ids[4][3] = {{0, 1, 2}, {0, 1, 3}, {1, 2, 3}, {0, 2, 3}}; - // Process edge-interior points for (int j = 0; j < edgeList.size(); ++j) { @@ -243,6 +391,9 @@ void OutputGmsh::Process() if (faceList[j]->m_vertexList.size() == 3) { + // TODO: move face_ids into a member function of Element + int face_ids[4][3] = { + {0, 1, 2}, {0, 1, 3}, {1, 2, 3}, {0, 2, 3}}; vector faceIds(3), volFaceIds(3); for (int k = 0; k < 3; ++k) @@ -261,10 +412,24 @@ void OutputGmsh::Process() } else { - // NEEDS ALIGNMENT - for (int k = 0; k < nodeList.size(); ++k) + // TODO: move face_ids into a member function of Element + 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 faceIds(4), volFaceIds(4); + + for (int k = 0; k < 4; ++k) { - tags.push_back(nodeList[k]->m_id); + faceIds [k] = faceList[j]->m_vertexList[k]->m_id; + volFaceIds[k] = e->GetVertexList()[face_ids[j][k]]->m_id; + } + + HOQuadrilateral hoQuad(faceIds, nodeList); + hoQuad.Align(volFaceIds); + + for (int k = 0; k < hoQuad.surfVerts.size(); ++k) + { + tags.push_back(hoQuad.surfVerts[k]->m_id); } } } @@ -279,8 +444,6 @@ void OutputGmsh::Process() vector reordering = InputGmsh::CreateReordering(elmtType); vector inv(tags.size()); - cout << tags.size() << endl; - for (int j = 0; j < tags.size(); ++j) { inv[reordering[j]] = j; -- GitLab From 0c89bf2db4ec1ef514ded3ac774dd6e64ec07e80 Mon Sep 17 00:00:00 2001 From: David Moxey Date: Fri, 5 Aug 2016 16:58:28 +0100 Subject: [PATCH 09/25] Hexes to order 4 working, cannot test further without Gmsh fixes --- utilities/NekMesh/InputModules/InputGmsh.cpp | 34 +++++++++++-------- .../NekMesh/OutputModules/OutputGmsh.cpp | 2 +- 2 files changed, 20 insertions(+), 16 deletions(-) diff --git a/utilities/NekMesh/InputModules/InputGmsh.cpp b/utilities/NekMesh/InputModules/InputGmsh.cpp index dc21063e0..031b5790b 100644 --- a/utilities/NekMesh/InputModules/InputGmsh.cpp +++ b/utilities/NekMesh/InputModules/InputGmsh.cpp @@ -320,12 +320,12 @@ std::vector tetTensorNodeOrdering(const std::vector &nodes, int n) } // Finally, recurse on interior volume - vector intNodes; + vector intNodes, tmpInt; for (int i = offset; i < nTet; ++i) { intNodes.push_back(nodes[i]); } - intNodes = tetTensorNodeOrdering(intNodes, n-4); + tmpInt = tetTensorNodeOrdering(intNodes, n-4); for (int k = 1, cnt = 0; k < n - 2; ++k) { @@ -333,7 +333,7 @@ std::vector tetTensorNodeOrdering(const std::vector &nodes, int n) { for (int i = 1; i < n - k - j - 1; ++i) { - nodeList[tmp[Mode(i,j,k)]] = intNodes[cnt++]; + nodeList[tmp[Mode(i,j,k)]] = tmpInt[cnt++]; } } } @@ -349,6 +349,8 @@ std::vector hexTensorNodeOrdering(const std::vector &nodes, int n) nodeList.resize(nodes.size()); nodeList[0] = nodes[0]; + cout << "BEING CALED WITH n = " << n << " nodes.size = " << nodes.size() << endl; + if (n == 1) { return nodeList; @@ -368,17 +370,17 @@ std::vector hexTensorNodeOrdering(const std::vector &nodes, int n) return nodeList; } - static int hexEdges[12][2] = { + int hexEdges[12][2] = { { 0, 1 }, { n-1, n }, { n*n-1, -1 }, { n*(n-1), -n }, { 0, n*n }, { n-1, n*n }, { n*n - 1, n*n }, { n*(n-1), n*n }, { n*n*(n-1), 1 }, { n*n*(n-1) + n-1, n }, { n*n*n-1, -1 }, { n*n*(n-1) + n*(n-1), -n } }; - static int hexFaces[6][3] = { + int hexFaces[6][3] = { { 0, 1, n }, { 0, 1, n*n }, { n-1, n, n*n }, { n*(n-1), 1, n*n }, { 0, n, n*n }, { n*n*(n-1), 1, n } }; - static int gmshToNekEdge[12] = {0, -3, 4, 1, 5, 2, 6, 7, 8, -11, 9, 10}; + int gmshToNekEdge[12] = {0, -3, 4, 1, 5, 2, 6, 7, 8, -11, 9, 10}; // Edges int offset = 8; @@ -403,7 +405,7 @@ std::vector hexTensorNodeOrdering(const std::vector &nodes, int n) } // Faces - static int gmsh2NekFace[6] = {0, 1, 4, 2, 3, 5}; + int gmsh2NekFace[6] = {0, 1, 4, 2, 3, 5}; // Map which defines orientation between Gmsh and Nektar++ faces. StdRegions::Orientation faceOrient[6] = { @@ -467,27 +469,29 @@ std::vector hexTensorNodeOrdering(const std::vector &nodes, int n) for (j = 1; j < n-1; ++j) { nodeList[hexFaces[face][0] + j*hexFaces[face][1] + k*hexFaces[face][2]] - = faceNodes[(j-1)*(n-2) + j - 1]; + = faceNodes[(k-1)*(n-2) + j-1]; } } } // Finally, recurse on interior volume - vector intNodes; + vector intNodes, tmpInt; for (int i = 8 + 12 * (n-2) + 6 * (n-2) * (n-2); i < n*n*n; ++i) { intNodes.push_back(nodes[i]); } - intNodes = hexTensorNodeOrdering(intNodes, n-2); - for (int k = 1, cnt = 0; k < n - 1; ++k) + if (intNodes.size()) { - for (int j = 1; j < n - 1; ++j) + tmpInt = hexTensorNodeOrdering(intNodes, n-2); + for (int k = 1, cnt = 0; k < n - 1; ++k) { - for (int i = 1; i < n - 1; ++i) + for (int j = 1; j < n - 1; ++j) { - cout << "LOLINT " << i + j * n + k * n * n << " " << intNodes[cnt] << endl; - nodeList[i + j * n + k * n * n] = intNodes[cnt++]; + for (int i = 1; i < n - 1; ++i) + { + nodeList[i + j * n + k * n * n] = tmpInt[cnt++]; + } } } } diff --git a/utilities/NekMesh/OutputModules/OutputGmsh.cpp b/utilities/NekMesh/OutputModules/OutputGmsh.cpp index 0ad0baf86..ec945dc2d 100644 --- a/utilities/NekMesh/OutputModules/OutputGmsh.cpp +++ b/utilities/NekMesh/OutputModules/OutputGmsh.cpp @@ -264,7 +264,7 @@ void OutputGmsh::Process() // Convert this mesh into a high-order mesh of uniform order. cout << "Mesh order of " << maxOrder << " detected" << endl; - maxOrder = 4; + maxOrder = 7; m_mesh->MakeOrder(maxOrder, LibUtilities::ePolyEvenlySpaced); // Add edge- and face-interior nodes to vertex set. -- GitLab From e57322bfe0c912365b215e5d7b6ad1c189e325bf Mon Sep 17 00:00:00 2001 From: David Moxey Date: Fri, 5 Aug 2016 18:09:45 +0100 Subject: [PATCH 10/25] Move new HOQuadrilateral class to Quadrilateral.h, prism edge mapping working --- library/NekMeshUtils/MeshElements/Prism.cpp | 13 +- .../NekMeshUtils/MeshElements/Quadrilateral.h | 153 ++++++++++++++++++ utilities/NekMesh/InputModules/InputGmsh.cpp | 144 +++++++++++++---- .../NekMesh/OutputModules/OutputGmsh.cpp | 152 +---------------- 4 files changed, 275 insertions(+), 187 deletions(-) diff --git a/library/NekMeshUtils/MeshElements/Prism.cpp b/library/NekMeshUtils/MeshElements/Prism.cpp index 8ac671a4f..d5b86437b 100644 --- a/library/NekMeshUtils/MeshElements/Prism.cpp +++ b/library/NekMeshUtils/MeshElements/Prism.cpp @@ -185,8 +185,17 @@ Prism::Prism(ElmtConfig pConf, faceNodes.push_back(pNodeList[face_offset[face] + i]); } } - m_face.push_back(FaceSharedPtr(new Face( - faceVertices, faceNodes, faceEdges, m_conf.m_faceCurveType))); + + // Try to translate between common face curve types + LibUtilities::PointsType pType = m_conf.m_faceCurveType; + + if (pType == LibUtilities::ePolyEvenlySpaced && (j == 1 || j == 3)) + { + pType = LibUtilities::eNodalTriEvenlySpaced; + } + + m_face.push_back( + FaceSharedPtr(new Face(faceVertices, faceNodes, faceEdges, pType))); } // Re-order edge array to be consistent with Nektar++ ordering. diff --git a/library/NekMeshUtils/MeshElements/Quadrilateral.h b/library/NekMeshUtils/MeshElements/Quadrilateral.h index e0f728f79..3819ad323 100644 --- a/library/NekMeshUtils/MeshElements/Quadrilateral.h +++ b/library/NekMeshUtils/MeshElements/Quadrilateral.h @@ -43,6 +43,7 @@ namespace Nektar { namespace NekMeshUtils { + /** * @brief A 2-dimensional four-sided element. */ @@ -87,6 +88,158 @@ public: NEKMESHUTILS_EXPORT static unsigned int GetNumNodes(ElmtConfig pConf); }; + +/** + * @brief A lightweight struct for dealing with high-order quadrilateral + * alignment. + */ +template struct HOQuadrilateral +{ + HOQuadrilateral(std::vector pVertId, std::vector pSurfVerts) + : vertId(pVertId), surfVerts(pSurfVerts) + { + } + + HOQuadrilateral(std::vector pVertId) : vertId(pVertId) + { + } + + /// The quadrilateral vertex IDs + std::vector vertId; + + /// The quadrilateral surface vertices -- templated so that this can either + /// be nodes or IDs. + std::vector surfVerts; + + void ReverseX() + { + int np = (int)(sqrt(surfVerts.size()) + 0.5); + for (int i = 0; i < np; ++i) + { + for (int j = 0; j < np/2; ++j) + { + swap(surfVerts[i*np + j], surfVerts[i*np + np-j-1]); + } + } + } + + void ReverseY() + { + int np = (int)(sqrt(surfVerts.size()) + 0.5); + // Reverse y direction + for (int j = 0; j < np; ++j) + { + for (int i = 0; i < np/2; ++i) + { + swap(surfVerts[i*np + j], surfVerts[(np-i-1)*np + j]); + } + } + } + + void Transpose() + { + int np = (int)(sqrt(surfVerts.size()) + 0.5); + std::vector tmp(surfVerts.size()); + + for (int i = 0; i < np; ++i) + { + for (int j = 0; j < np; ++j) + { + tmp[i*np+j] = surfVerts[j*np+i]; + } + } + + surfVerts = tmp; + } + + /** + * @brief Align this surface to a given vertex ID. + */ + void Align(std::vector vertId) + { + int vmap[4] = {-1, -1, -1, -1}; + + // Determine which vertices map to vertId + for (int i = 0; i < 4; ++i) + { + for (int j = 0; j < 4; ++j) + { + if (this->vertId[j] == vertId[i]) + { + vmap[i] = j; + break; + } + } + + ASSERTL1(vmap[i] != -1, + "Could not determine mapping between vertex IDs"); + } + + StdRegions::Orientation orient = StdRegions::eNoOrientation; + + if (vmap[1] == (vmap[0]+1) % 4) + { + switch (vmap[0]) + { + case 0: + orient = StdRegions::eDir1FwdDir1_Dir2FwdDir2; + break; + case 1: + orient = StdRegions::eDir1BwdDir2_Dir2FwdDir1; + break; + case 2: + orient = StdRegions::eDir1BwdDir1_Dir2BwdDir2; + break; + case 3: + orient = StdRegions::eDir1FwdDir2_Dir2BwdDir1; + break; + } + } + else + { + switch (vmap[0]) + { + case 0: + orient = StdRegions::eDir1FwdDir2_Dir2FwdDir1; + break; + case 1: + orient = StdRegions::eDir1BwdDir1_Dir2FwdDir2; + break; + case 2: + orient = StdRegions::eDir1BwdDir2_Dir2BwdDir1; + break; + case 3: + orient = StdRegions::eDir1FwdDir1_Dir2BwdDir2; + break; + } + } + + if (orient == StdRegions::eDir1FwdDir2_Dir2FwdDir1 || + orient == StdRegions::eDir1BwdDir2_Dir2FwdDir1 || + orient == StdRegions::eDir1FwdDir2_Dir2BwdDir1 || + orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1) + { + Transpose(); + } + + if (orient == StdRegions::eDir1BwdDir1_Dir2FwdDir2 || + orient == StdRegions::eDir1BwdDir1_Dir2BwdDir2 || + orient == StdRegions::eDir1FwdDir2_Dir2BwdDir1 || + orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1) + { + ReverseX(); + } + + if (orient == StdRegions::eDir1FwdDir1_Dir2BwdDir2 || + orient == StdRegions::eDir1BwdDir1_Dir2BwdDir2 || + orient == StdRegions::eDir1BwdDir2_Dir2FwdDir1 || + orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1) + { + ReverseY(); + } + } +}; + } } diff --git a/utilities/NekMesh/InputModules/InputGmsh.cpp b/utilities/NekMesh/InputModules/InputGmsh.cpp index 031b5790b..a55a31ea9 100644 --- a/utilities/NekMesh/InputModules/InputGmsh.cpp +++ b/utilities/NekMesh/InputModules/InputGmsh.cpp @@ -1119,7 +1119,7 @@ vector InputGmsh::PrismReordering(ElmtConfig conf) // different; need to mirror in the triangular faces, and then // reorder vertices to make ordering anticlockwise on base quad. static int gmshToNekVerts[6] = {3, 4, 1, 0, 5, 2}; - // 3->0, 4->1, 1->2, 0->3, 5->4, 2->5 + // Inverse (gmsh vert -> Nektar++ vert): 3->0, 4->1, 1->2, 0->3, 5->4, 2->5 for (i = 0; i < 6; ++i) { @@ -1134,8 +1134,18 @@ vector InputGmsh::PrismReordering(ElmtConfig conf) // Curvilinear edges. mapping.resize(6 + 9 * n); + // Nektar++ edges: + // 0 = 1 = 2 = 3 = 4 = 5 = 6 = 7 = 8 = + // {0,1}, {1,2}, {3,2}, {0,3}, {0,4}, {1,4}, {2,5}, {3,5}, {4,5} + + // Gmsh edges: + // {0,1}, {0,2}, {0,3}, {1,2}, {1,4}, {2,5}, {3,4}, {3,5}, {4,5} + // Apply inverse mapping: + // {3,2}, {3,5}, {3,0}, {2,5}, {2,1}, {5,4}, {0,1}, {0,4}, {1,4} + // = 2 = 7 = -3 = 6 = -1 = -8 = 0 = 4 = 5 + static int gmshToNekEdge[9] = {2, 7, 3, 6, 1, 8, 0, 4, 5}; - static int gmshToNekRev[9] = {1, 0, 1, 0, 1, 1, 0, 0, 0}; + static int gmshToNekRev[9] = {0, 0, 1, 0, 1, 1, 0, 0, 0}; // Reorder edges. int offset, cnt = 6; @@ -1182,7 +1192,7 @@ vector InputGmsh::PrismReordering(ElmtConfig conf) // Face 4: StdRegions::eDir1BwdDir1_Dir2FwdDir2 // Face offsets - vector offsets(5); + vector offsets(5), offsets2(5); offset = 6 + 9 * n; offsets[0] = offset; offsets[1] = offset + nTriInt; @@ -1190,52 +1200,118 @@ vector InputGmsh::PrismReordering(ElmtConfig conf) offsets[3] = offset + 2 * nTriInt + nQuadInt; offsets[4] = offset + 2 * nTriInt + 2 * nQuadInt; + offsets2[0] = offset; + offsets2[1] = offset + nQuadInt; + offsets2[2] = offset + nQuadInt + nTriInt; + offsets2[3] = offset + 2 * nQuadInt + nTriInt; + offsets2[4] = offset + 2 * nQuadInt + 2 * nTriInt; + mapping.resize(6 + 9 * n + 3 * nQuadInt + 2 * nTriInt); offset = 6 + 9 * n; - // Face 0 - int offset2 = offsets[gmshToNekFace[0]]; - int j; - for (i = 0; i < n; ++i) + vector triVertId(3); + triVertId[0] = 0; + triVertId[1] = 1; + triVertId[2] = 2; + + vector quadVertId(4); + quadVertId[0] = 0; + quadVertId[1] = 1; + quadVertId[2] = 2; + quadVertId[3] = 3; + + for (i = 0; i < 5; ++i) { - for (j = 0; j < n; ++j) + int face = gmshToNekFace[i]; + int offset2 = offsets[i]; + offset = offsets[face]; + + bool tri = i < 2; + int nFacePts = tri ? nTriInt : nQuadInt; + + if (nFacePts == 0) { - mapping[offset + i * n + j] = offset2 + (n - i - 1) * n + j; + continue; } - } - offset += nQuadInt; - // Face 1 - offset2 = offsets[gmshToNekFace[1]]; - for (i = 0; i < nTriInt; ++i) - { - mapping[offset++] = offset2++; - } + // Create a list of interior face nodes for this face only. + vector faceNodes(nFacePts); + vector toAlign(tri ? 3 : 4); + for (int j = 0; j < nFacePts; ++j) + { + faceNodes[j] = offset2 + j; + } - // Face 1 - offset2 = offsets[gmshToNekFace[2]]; - for (i = 0; i < nQuadInt; ++i) - { - mapping[offset++] = offset2++; - cout << offset-1 << " " << offset2-1 << endl; - } + if (tri) + { + // Now get the reordering of this face, which puts Gmsh + // recursive ordering into Nektar++ row-by-row order. + vector tmp = triTensorNodeOrdering(faceNodes, n - 1); + HOTriangle hoTri(triVertId, tmp); - // Face 1 - offset2 = offsets[gmshToNekFace[3]]; - for (i = 0; i < nTriInt; ++i) - { - mapping[offset++] = offset2++; + // Apply reorientation + if (i == 0) + { + // Triangle verts {0,2,1} --> {0,1,2} + toAlign[0] = 0; + toAlign[1] = 2; + toAlign[2] = 1; + hoTri.Align(toAlign); + } + + // Fill in mapping. + for (int j = 0; j < nTriInt; ++j) + { + mapping[offset + j] = hoTri.surfVerts[j]; + } + } + else + { + vector tmp = quadTensorNodeOrdering(faceNodes, n - 1); + HOQuadrilateral hoQuad(quadVertId, tmp); + + // Apply reorientation + if (i == 2) + { + toAlign[0] = 3; + toAlign[1] = 2; + toAlign[2] = 1; + toAlign[3] = 0; + } + else if (i == 3) + { + toAlign[0] = 1; + toAlign[1] = 0; + toAlign[2] = 3; + toAlign[3] = 2; + } + else if (i == 4) + { + toAlign[0] = 1; + toAlign[1] = 2; + toAlign[2] = 3; + toAlign[3] = 0; + } + + hoQuad.Align(toAlign); + + // Fill in mapping. + for (int j = 0; j < nQuadInt; ++j) + { + mapping[offset + j] = hoQuad.surfVerts[j]; + } + } } - // Face 1 - offset2 = offsets[gmshToNekFace[4]]; - for (i = 0; i < nQuadInt; ++i) + if (conf.m_volumeNodes == false) { - mapping[offset++] = offset2++; - cout << offset-1 << " " << offset2-1 << endl; + return mapping; } + // Interior points + offset = offsets[4]; + return mapping; } diff --git a/utilities/NekMesh/OutputModules/OutputGmsh.cpp b/utilities/NekMesh/OutputModules/OutputGmsh.cpp index ec945dc2d..2dee31787 100644 --- a/utilities/NekMesh/OutputModules/OutputGmsh.cpp +++ b/utilities/NekMesh/OutputModules/OutputGmsh.cpp @@ -35,6 +35,7 @@ #include #include +#include #include "OutputGmsh.h" #include "../InputModules/InputGmsh.h" @@ -52,157 +53,6 @@ ModuleKey OutputGmsh::className = OutputGmsh::create, "Writes Gmsh msh file."); -/** - * @brief A lightweight struct for dealing with high-order quadrilateral - * alignment. - */ -template struct HOQuadrilateral -{ - HOQuadrilateral(std::vector pVertId, std::vector pSurfVerts) - : vertId(pVertId), surfVerts(pSurfVerts) - { - } - - HOQuadrilateral(std::vector pVertId) : vertId(pVertId) - { - } - - /// The quadrilateral vertex IDs - std::vector vertId; - - /// The quadrilateral surface vertices -- templated so that this can either - /// be nodes or IDs. - std::vector surfVerts; - - void ReverseX() - { - int np = (int)(sqrt(surfVerts.size()) + 0.5); - for (int i = 0; i < np; ++i) - { - for (int j = 0; j < np/2; ++j) - { - swap(surfVerts[i*np + j], surfVerts[i*np + np-j-1]); - } - } - } - - void ReverseY() - { - int np = (int)(sqrt(surfVerts.size()) + 0.5); - // Reverse y direction - for (int j = 0; j < np; ++j) - { - for (int i = 0; i < np/2; ++i) - { - swap(surfVerts[i*np + j], surfVerts[(np-i-1)*np + j]); - } - } - } - - void Transpose() - { - int np = (int)(sqrt(surfVerts.size()) + 0.5); - vector tmp(surfVerts.size()); - - for (int i = 0; i < np; ++i) - { - for (int j = 0; j < np; ++j) - { - tmp[i*np+j] = surfVerts[j*np+i]; - } - } - - surfVerts = tmp; - } - - /** - * @brief Align this surface to a given vertex ID. - */ - void Align(std::vector vertId) - { - int vmap[4] = {-1, -1, -1, -1}; - - // Determine which vertices map to vertId - for (int i = 0; i < 4; ++i) - { - for (int j = 0; j < 4; ++j) - { - if (this->vertId[j] == vertId[i]) - { - vmap[i] = j; - break; - } - } - - ASSERTL1(vmap[i] != -1, - "Could not determine mapping between vertex IDs"); - } - - StdRegions::Orientation orient = StdRegions::eNoOrientation; - - if (vmap[1] == (vmap[0]+1) % 4) - { - switch (vmap[0]) - { - case 0: - orient = StdRegions::eDir1FwdDir1_Dir2FwdDir2; - break; - case 1: - orient = StdRegions::eDir1BwdDir2_Dir2FwdDir1; - break; - case 2: - orient = StdRegions::eDir1BwdDir1_Dir2BwdDir2; - break; - case 3: - orient = StdRegions::eDir1FwdDir2_Dir2BwdDir1; - break; - } - } - else - { - switch (vmap[0]) - { - case 0: - orient = StdRegions::eDir1FwdDir2_Dir2FwdDir1; - break; - case 1: - orient = StdRegions::eDir1BwdDir1_Dir2FwdDir2; - break; - case 2: - orient = StdRegions::eDir1BwdDir2_Dir2BwdDir1; - break; - case 3: - orient = StdRegions::eDir1FwdDir1_Dir2BwdDir2; - break; - } - } - - if (orient == StdRegions::eDir1FwdDir2_Dir2FwdDir1 || - orient == StdRegions::eDir1BwdDir2_Dir2FwdDir1 || - orient == StdRegions::eDir1FwdDir2_Dir2BwdDir1 || - orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1) - { - Transpose(); - } - - if (orient == StdRegions::eDir1BwdDir1_Dir2FwdDir2 || - orient == StdRegions::eDir1BwdDir1_Dir2BwdDir2 || - orient == StdRegions::eDir1FwdDir2_Dir2BwdDir1 || - orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1) - { - ReverseX(); - } - - if (orient == StdRegions::eDir1FwdDir1_Dir2BwdDir2 || - orient == StdRegions::eDir1BwdDir1_Dir2BwdDir2 || - orient == StdRegions::eDir1BwdDir2_Dir2FwdDir1 || - orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1) - { - ReverseY(); - } - } -}; - OutputGmsh::OutputGmsh(MeshSharedPtr m) : OutputModule(m) { map igelmap = InputGmsh::GenElmMap(); -- GitLab From 46b3826b2bff7f7d07d4ef65c1b28179223df486 Mon Sep 17 00:00:00 2001 From: David Moxey Date: Fri, 5 Aug 2016 22:50:13 +0100 Subject: [PATCH 11/25] Working prism reordering up to order 5, gmsh does not seem to generate higher --- utilities/NekMesh/InputModules/InputGmsh.cpp | 42 ++++++++++---------- 1 file changed, 20 insertions(+), 22 deletions(-) diff --git a/utilities/NekMesh/InputModules/InputGmsh.cpp b/utilities/NekMesh/InputModules/InputGmsh.cpp index a55a31ea9..3d9b7e725 100644 --- a/utilities/NekMesh/InputModules/InputGmsh.cpp +++ b/utilities/NekMesh/InputModules/InputGmsh.cpp @@ -1118,12 +1118,12 @@ vector InputGmsh::PrismReordering(ElmtConfig conf) // To get from Gmsh -> Nektar++ prism, coordinates axes are // different; need to mirror in the triangular faces, and then // reorder vertices to make ordering anticlockwise on base quad. - static int gmshToNekVerts[6] = {3, 4, 1, 0, 5, 2}; + static int nekToGmshVerts[6] = {3, 4, 1, 0, 5, 2}; // Inverse (gmsh vert -> Nektar++ vert): 3->0, 4->1, 1->2, 0->3, 5->4, 2->5 for (i = 0; i < 6; ++i) { - mapping[i] = gmshToNekVerts[i]; + mapping[i] = nekToGmshVerts[i]; } if (order == 1) @@ -1137,13 +1137,12 @@ vector InputGmsh::PrismReordering(ElmtConfig conf) // Nektar++ edges: // 0 = 1 = 2 = 3 = 4 = 5 = 6 = 7 = 8 = // {0,1}, {1,2}, {3,2}, {0,3}, {0,4}, {1,4}, {2,5}, {3,5}, {4,5} - - // Gmsh edges: + // Gmsh edges (from Geo/MPrism.h of Gmsh source): // {0,1}, {0,2}, {0,3}, {1,2}, {1,4}, {2,5}, {3,4}, {3,5}, {4,5} - // Apply inverse mapping: + // Apply inverse of gmshToNekVerts map: // {3,2}, {3,5}, {3,0}, {2,5}, {2,1}, {5,4}, {0,1}, {0,4}, {1,4} + // Nektar++ mapping (negative indicates reverse orientation): // = 2 = 7 = -3 = 6 = -1 = -8 = 0 = 4 = 5 - static int gmshToNekEdge[9] = {2, 7, 3, 6, 1, 8, 0, 4, 5}; static int gmshToNekRev[9] = {0, 0, 1, 0, 1, 1, 0, 0, 0}; @@ -1177,29 +1176,28 @@ vector InputGmsh::PrismReordering(ElmtConfig conf) int nTriInt = n * (n - 1) / 2; int nQuadInt = n * n; - // Gmsh faces: - // 0 = {0,2,1}, 1 = {3,4,5}, 2 = {0,1,4,3}, 3 = {0,3,5,2}, 4 = {1,2,5,4} - // Apply gmshToNekVerts map: - // 0 = {3,5,2}, 1 = {0,1,4}, 2 = {3,2,1,0}, 3 = {3,0,4,5}, 4 = {2,5,4,1} // Nektar++ faces: // 0 = {0,1,2,3}, 1 = {0,1,4}, 2 = {1,2,5,4}, 3 = {3,2,5}, 4 = {0,3,5,4} - // {3,2,1,0} {0,1,4} {2,5,4,1} {3,5,2} {3,0,4,5} + // Gmsh faces (from Geo/MPrism.h of Gmsh source): + // {0,2,1}, {3,4,5}, {0,1,4,3}, {0,3,5,2}, {1,2,5,4} + // Apply inverse of gmshToNekVerts map: + // {3,5,2}, {0,1,4}, {3,2,1,0}, {3,0,4,5}, {2,5,4,1} + // = 3 = 1 = 0 = 4 = 2 // This gives gmsh -> Nektar++ faces: - static int gmshToNekFace[5] = {2, 1, 4, 0, 3}; - - // Face 0: StdRegions::eDir1FwdDir1_Dir2BwdDir2 - // Face 2: StdRegions::eDir1BwdDir2_Dir2FwdDir1 - // Face 4: StdRegions::eDir1BwdDir1_Dir2FwdDir2 + static int gmshToNekFace[5] = {3, 1, 0, 4, 2}; // Face offsets vector offsets(5), offsets2(5); offset = 6 + 9 * n; + + // Offsets in the gmsh order: TTQQQ offsets[0] = offset; offsets[1] = offset + nTriInt; offsets[2] = offset + 2 * nTriInt; offsets[3] = offset + 2 * nTriInt + nQuadInt; offsets[4] = offset + 2 * nTriInt + 2 * nQuadInt; + // Offsets in the Nektar++ order: QTQTQ offsets2[0] = offset; offsets2[1] = offset + nQuadInt; offsets2[2] = offset + nQuadInt + nTriInt; @@ -1225,7 +1223,7 @@ vector InputGmsh::PrismReordering(ElmtConfig conf) { int face = gmshToNekFace[i]; int offset2 = offsets[i]; - offset = offsets[face]; + offset = offsets2[face]; bool tri = i < 2; int nFacePts = tri ? nTriInt : nQuadInt; @@ -1268,7 +1266,7 @@ vector InputGmsh::PrismReordering(ElmtConfig conf) } else { - vector tmp = quadTensorNodeOrdering(faceNodes, n - 1); + vector tmp = quadTensorNodeOrdering(faceNodes, n); HOQuadrilateral hoQuad(quadVertId, tmp); // Apply reorientation @@ -1288,10 +1286,10 @@ vector InputGmsh::PrismReordering(ElmtConfig conf) } else if (i == 4) { - toAlign[0] = 1; - toAlign[1] = 2; - toAlign[2] = 3; - toAlign[3] = 0; + toAlign[0] = 3; + toAlign[1] = 0; + toAlign[2] = 1; + toAlign[3] = 2; } hoQuad.Align(toAlign); -- GitLab From de97118cea34f4349dac855965caee9bafb15d24 Mon Sep 17 00:00:00 2001 From: David Moxey Date: Tue, 9 Aug 2016 18:48:47 +0100 Subject: [PATCH 12/25] Start to add more support for prisms --- library/NekMeshUtils/MeshElements/Mesh.cpp | 10 + library/NekMeshUtils/MeshElements/Prism.cpp | 176 +++++++----------- library/NekMeshUtils/MeshElements/Prism.h | 10 +- .../NekMeshUtils/MeshElements/Tetrahedron.cpp | 3 +- utilities/NekMesh/InputModules/InputGmsh.cpp | 63 ++++++- .../NekMesh/OutputModules/OutputGmsh.cpp | 16 +- 6 files changed, 158 insertions(+), 120 deletions(-) diff --git a/library/NekMeshUtils/MeshElements/Mesh.cpp b/library/NekMeshUtils/MeshElements/Mesh.cpp index d24768667..50018938f 100644 --- a/library/NekMeshUtils/MeshElements/Mesh.cpp +++ b/library/NekMeshUtils/MeshElements/Mesh.cpp @@ -104,6 +104,16 @@ void Mesh::MakeOrder(int order, 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; + } + 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; } diff --git a/library/NekMeshUtils/MeshElements/Prism.cpp b/library/NekMeshUtils/MeshElements/Prism.cpp index d5b86437b..130da17ec 100644 --- a/library/NekMeshUtils/MeshElements/Prism.cpp +++ b/library/NekMeshUtils/MeshElements/Prism.cpp @@ -251,124 +251,90 @@ SpatialDomains::GeometrySharedPtr Prism::GetGeom(int coordDim) return ret; } -/** - * @brief . - */ -void Prism::Complete(int order) +StdRegions::Orientation Prism::GetEdgeOrient( + int edgeId, EdgeSharedPtr edge) { - int i, j, pos; - - // 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_A, - order + 1, - LibUtilities::PointsKey(order + 1, - LibUtilities::eGaussLobattoLegendre)); - LibUtilities::BasisKey B2( - LibUtilities::eOrtho_B, - order + 1, - LibUtilities::PointsKey(order + 1, - LibUtilities::eGaussRadauMAlpha1Beta0)); - - // Create a standard nodal prism in order to get the Vandermonde - // matrix to perform interpolation to nodal points. - StdRegions::StdNodalPrismExpSharedPtr nodalPrism = - MemoryManager::AllocateSharedPtr( - B0, B1, B2, LibUtilities::eNodalPrismEvenlySpaced); - - Array x, y, z; - nodalPrism->GetNodalPoints(x, y, z); - - SpatialDomains::PrismGeomSharedPtr geom = - boost::dynamic_pointer_cast( - this->GetGeom(3)); - - // Create basis key for a prism. - LibUtilities::BasisKey C0( - LibUtilities::eOrtho_A, - order + 1, - LibUtilities::PointsKey(order + 1, - LibUtilities::eGaussLobattoLegendre)); - LibUtilities::BasisKey C1( - LibUtilities::eOrtho_A, - order + 1, - LibUtilities::PointsKey(order + 1, - LibUtilities::eGaussLobattoLegendre)); - LibUtilities::BasisKey C2( - LibUtilities::eOrtho_B, - order + 1, - LibUtilities::PointsKey(order + 1, - LibUtilities::eGaussRadauMAlpha1Beta0)); - - // Create a prism. - LocalRegions::PrismExpSharedPtr prism = - MemoryManager::AllocateSharedPtr( - C0, C1, C2, geom); - - // Get coordinate array for tetrahedron. - int nqtot = prism->GetTotPoints(); - Array alloc(6 * nqtot); - Array xi(alloc); - Array yi(alloc + nqtot); - Array zi(alloc + 2 * nqtot); - Array xo(alloc + 3 * nqtot); - Array yo(alloc + 4 * nqtot); - Array zo(alloc + 5 * nqtot); - Array tmp; - - prism->GetCoords(xi, yi, zi); - - for (i = 0; i < 3; ++i) + static int edgeVerts[9][2] = { + {0,1}, {1,2}, {3,2}, {0,3}, {0,4}, {1,4}, {2,5}, {3,5}, {4,5} + }; + + if (edge->m_n1 == m_vertex[edgeVerts[edgeId][0]]) { - Array coeffs(nodalPrism->GetNcoeffs()); - prism->FwdTrans(alloc + i * nqtot, coeffs); - // Apply Vandermonde matrix to project onto nodal space. - nodalPrism->ModalToNodal(coeffs, tmp = alloc + (i + 3) * nqtot); + return StdRegions::eForwards; } - - // Now extract points from the co-ordinate arrays into the - // edge/face/volume nodes. First, extract edge-interior nodes. - for (i = 0; i < 9; ++i) + else if (edge->m_n1 == m_vertex[edgeVerts[edgeId][1]]) { - pos = 6 + 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]))); - } + return StdRegions::eBackwards; } - - // Now extract face-interior nodes. - pos = 6 + 9 * (order - 1); - for (i = 0; i < 5; ++i) + else { - int facesize = - i % 2 ? (order - 2) * (order - 1) / 2 : (order - 1) * (order - 1); - m_face[i]->m_faceNodes.clear(); - for (j = 0; j < facesize; ++j) - { - m_face[i]->m_faceNodes.push_back(NodeSharedPtr( - new Node(0, xo[pos + j], yo[pos + j], zo[pos + j]))); - } - pos += facesize; + ASSERTL1(false, "Edge is not connected to this quadrilateral."); } - // Finally extract volume nodes. - for (i = pos; i < (order + 1) * (order + 1) * (order + 2) / 2; ++i) + return StdRegions::eNoOrientation; +} + +void Prism::MakeOrder(int order, + SpatialDomains::GeometrySharedPtr geom, + LibUtilities::PointsType pType, + int coordDim, + int &id) +{ + m_conf.m_order = order; + m_volumeNodes.clear(); + + if (order == 1) + { + m_conf.m_volumeNodes = m_conf.m_faceNodes = false; + return; + } + else if (order == 2) { - m_volumeNodes.push_back( - NodeSharedPtr(new Node(0, xo[i], yo[i], zo[i]))); + m_conf.m_faceNodes = true; + m_conf.m_volumeNodes = false; + return; } - m_conf.m_order = order; m_conf.m_faceNodes = true; m_conf.m_volumeNodes = true; + m_curveType = pType; + + int nPoints = order + 1; + StdRegions::StdExpansionSharedPtr xmap = geom->GetXmap(); + + Array px, py, pz; + LibUtilities::PointsKey pKey(nPoints, pType); + ASSERTL1(pKey.GetPointsDim() == 3, "Points distribution must be 3D"); + LibUtilities::PointsManager()[pKey]->GetPoints(px, py, pz); + + Array > phys(coordDim); + + for (int i = 0; i < coordDim; ++i) + { + phys[i] = Array(xmap->GetTotPoints()); + xmap->BwdTrans(geom->GetCoeffs(i), phys[i]); + } + + const int nPrismPts = nPoints * nPoints * (nPoints + 1) / 2; + const int nPrismIntPts = (nPoints - 2) * (nPoints - 3) * (nPoints - 2) / 2; + m_volumeNodes.resize(nPrismIntPts); + + for (int i = nPrismPts - nPrismIntPts, cnt = 0; i < nPrismPts; ++i, ++cnt) + { + Array xp(3); + xp[0] = px[i]; + xp[1] = py[i]; + xp[2] = pz[i]; + + Array x(3, 0.0); + for (int j = 0; j < coordDim; ++j) + { + x[j] = xmap->PhysEvaluate(xp, phys[j]); + } + + m_volumeNodes[cnt] = boost::shared_ptr( + new Node(id++, x[0], x[1], x[2])); + } } /** diff --git a/library/NekMeshUtils/MeshElements/Prism.h b/library/NekMeshUtils/MeshElements/Prism.h index dbe4b2da4..a0b1dd3bb 100644 --- a/library/NekMeshUtils/MeshElements/Prism.h +++ b/library/NekMeshUtils/MeshElements/Prism.h @@ -77,8 +77,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); /** diff --git a/library/NekMeshUtils/MeshElements/Tetrahedron.cpp b/library/NekMeshUtils/MeshElements/Tetrahedron.cpp index bdb67ab2c..fb632ee6e 100644 --- a/library/NekMeshUtils/MeshElements/Tetrahedron.cpp +++ b/library/NekMeshUtils/MeshElements/Tetrahedron.cpp @@ -301,8 +301,7 @@ void Tetrahedron::MakeOrder(int order, const int nTetIntPts = (nPoints - 4) * (nPoints - 3) * (nPoints - 2) / 6; m_volumeNodes.resize(nTetIntPts); - for (int i = 4 + 6*(nPoints-2) + 2*(nPoints-3)*(nPoints-2), cnt = 0; - i < nTetPts; ++i, ++cnt) + for (int i = nTetPts - nTetIntPts, cnt = 0; i < nTetPts; ++i, ++cnt) { Array xp(3); xp[0] = px[i]; diff --git a/utilities/NekMesh/InputModules/InputGmsh.cpp b/utilities/NekMesh/InputModules/InputGmsh.cpp index 3d9b7e725..db85caee0 100644 --- a/utilities/NekMesh/InputModules/InputGmsh.cpp +++ b/utilities/NekMesh/InputModules/InputGmsh.cpp @@ -341,6 +341,51 @@ std::vector tetTensorNodeOrdering(const std::vector &nodes, int n) return nodeList; } +std::vector prismTensorNodeOrdering(const std::vector &nodes, int n) +{ + std::vector nodeList; + int cnt2; + int nTri = n*(n+1)/2; + int nPrism = n*n*(n+1)/2; + + if (n == 0) + { + return nodeList; + } + + nodeList.resize(nodes.size()); + + if (n == 1) + { + 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) + { + return nodeList; + } + + return nodeList; +} + std::vector hexTensorNodeOrdering(const std::vector &nodes, int n) { int i, j, k; @@ -349,8 +394,6 @@ std::vector hexTensorNodeOrdering(const std::vector &nodes, int n) nodeList.resize(nodes.size()); nodeList[0] = nodes[0]; - cout << "BEING CALED WITH n = " << n << " nodes.size = " << nodes.size() << endl; - if (n == 1) { return nodeList; @@ -1190,14 +1233,14 @@ vector InputGmsh::PrismReordering(ElmtConfig conf) vector offsets(5), offsets2(5); offset = 6 + 9 * n; - // Offsets in the gmsh order: TTQQQ + // Offsets in the gmsh order: face ordering is TTQQQ offsets[0] = offset; offsets[1] = offset + nTriInt; offsets[2] = offset + 2 * nTriInt; offsets[3] = offset + 2 * nTriInt + nQuadInt; offsets[4] = offset + 2 * nTriInt + 2 * nQuadInt; - // Offsets in the Nektar++ order: QTQTQ + // Offsets in the Nektar++ order: face ordering is QTQTQ offsets2[0] = offset; offsets2[1] = offset + nQuadInt; offsets2[2] = offset + nQuadInt + nTriInt; @@ -1308,7 +1351,17 @@ vector InputGmsh::PrismReordering(ElmtConfig conf) } // Interior points - offset = offsets[4]; + offset = offsets[4] + nQuadInt; + vector intPoints, tmp; + + for (int i = offset; i < (order+1) * (order+1) * (order+2) / 2; ++i) + { + intPoints.push_back(i); + } + + // Reorder this stuff + tmp = prismTensorNodeOrdering(intPoints, order - 2); + mapping.insert(mapping.end(), tmp.begin(), tmp.end()); return mapping; } diff --git a/utilities/NekMesh/OutputModules/OutputGmsh.cpp b/utilities/NekMesh/OutputModules/OutputGmsh.cpp index 2dee31787..b7d7f0b4d 100644 --- a/utilities/NekMesh/OutputModules/OutputGmsh.cpp +++ b/utilities/NekMesh/OutputModules/OutputGmsh.cpp @@ -114,7 +114,7 @@ void OutputGmsh::Process() // Convert this mesh into a high-order mesh of uniform order. cout << "Mesh order of " << maxOrder << " detected" << endl; - maxOrder = 7; + maxOrder = 4; m_mesh->MakeOrder(maxOrder, LibUtilities::ePolyEvenlySpaced); // Add edge- and face-interior nodes to vertex set. @@ -242,8 +242,10 @@ void OutputGmsh::Process() if (faceList[j]->m_vertexList.size() == 3) { // TODO: move face_ids into a member function of Element - int face_ids[4][3] = { - {0, 1, 2}, {0, 1, 3}, {1, 2, 3}, {0, 2, 3}}; + 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 faceIds(3), volFaceIds(3); for (int k = 0; k < 3; ++k) @@ -263,9 +265,11 @@ void OutputGmsh::Process() else { // TODO: move face_ids into a member function of Element - 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}}; + 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 faceIds(4), volFaceIds(4); for (int k = 0; k < 4; ++k) -- GitLab From 5be44839a0098bd578cc13cc1cfe89ea15b6ddcb Mon Sep 17 00:00:00 2001 From: David Moxey Date: Thu, 11 Aug 2016 16:53:34 +0100 Subject: [PATCH 13/25] 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 | 7 ++ .../NekMeshUtils/MeshElements/Hexahedron.cpp | 18 +-- .../NekMeshUtils/MeshElements/Hexahedron.h | 7 ++ library/NekMeshUtils/MeshElements/Mesh.cpp | 12 +- library/NekMeshUtils/MeshElements/Prism.cpp | 17 ++- library/NekMeshUtils/MeshElements/Prism.h | 7 ++ library/NekMeshUtils/MeshElements/Pyramid.cpp | 16 +-- library/NekMeshUtils/MeshElements/Pyramid.h | 7 ++ .../NekMeshUtils/MeshElements/Tetrahedron.cpp | 44 +++---- .../NekMeshUtils/MeshElements/Tetrahedron.h | 7 ++ utilities/NekMesh/InputModules/InputGmsh.cpp | 107 +++++++++--------- .../NekMesh/OutputModules/OutputGmsh.cpp | 38 ++----- 12 files changed, 154 insertions(+), 133 deletions(-) diff --git a/library/NekMeshUtils/MeshElements/Element.h b/library/NekMeshUtils/MeshElements/Element.h index b7a36ac63..ec32cada2 100644 --- a/library/NekMeshUtils/MeshElements/Element.h +++ b/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; diff --git a/library/NekMeshUtils/MeshElements/Hexahedron.cpp b/library/NekMeshUtils/MeshElements/Hexahedron.cpp index cde1e7406..a307fa232 100644 --- a/library/NekMeshUtils/MeshElements/Hexahedron.cpp +++ b/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 faceVertices; @@ -115,9 +115,9 @@ Hexahedron::Hexahedron(ElmtConfig pConf, vector 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) && diff --git a/library/NekMeshUtils/MeshElements/Hexahedron.h b/library/NekMeshUtils/MeshElements/Hexahedron.h index 43e3f421c..9b38e7588 100644 --- a/library/NekMeshUtils/MeshElements/Hexahedron.h +++ b/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]; }; } } diff --git a/library/NekMeshUtils/MeshElements/Mesh.cpp b/library/NekMeshUtils/MeshElements/Mesh.cpp index 50018938f..ca5e94fdb 100644 --- a/library/NekMeshUtils/MeshElements/Mesh.cpp +++ b/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++) diff --git a/library/NekMeshUtils/MeshElements/Prism.cpp b/library/NekMeshUtils/MeshElements/Prism.cpp index 130da17ec..25943c094 100644 --- a/library/NekMeshUtils/MeshElements/Prism.cpp +++ b/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 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 x(3, 0.0); for (int j = 0; j < coordDim; ++j) { diff --git a/library/NekMeshUtils/MeshElements/Prism.h b/library/NekMeshUtils/MeshElements/Prism.h index a0b1dd3bb..65b4c506b 100644 --- a/library/NekMeshUtils/MeshElements/Prism.h +++ b/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]; }; } } diff --git a/library/NekMeshUtils/MeshElements/Pyramid.cpp b/library/NekMeshUtils/MeshElements/Pyramid.cpp index 0c15b0878..59ca3f543 100644 --- a/library/NekMeshUtils/MeshElements/Pyramid.cpp +++ b/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) || diff --git a/library/NekMeshUtils/MeshElements/Pyramid.h b/library/NekMeshUtils/MeshElements/Pyramid.h index 815d19bcc..76ba07ee4 100644 --- a/library/NekMeshUtils/MeshElements/Pyramid.h +++ 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]; }; } } diff --git a/library/NekMeshUtils/MeshElements/Tetrahedron.cpp b/library/NekMeshUtils/MeshElements/Tetrahedron.cpp index fb632ee6e..895b4422b 100644 --- a/library/NekMeshUtils/MeshElements/Tetrahedron.cpp +++ b/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 faceIds(3); - faceIds[0] = faceVertices[0]->m_id; - faceIds[1] = faceVertices[1]->m_id; - faceIds[2] = faceVertices[2]->m_id; + vector 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 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 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 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 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); diff --git a/library/NekMeshUtils/MeshElements/Tetrahedron.h b/library/NekMeshUtils/MeshElements/Tetrahedron.h index cc59db66a..d8d1e4cf0 100644 --- a/library/NekMeshUtils/MeshElements/Tetrahedron.h +++ b/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]; }; } } diff --git a/utilities/NekMesh/InputModules/InputGmsh.cpp b/utilities/NekMesh/InputModules/InputGmsh.cpp index db85caee0..5f909662f 100644 --- a/utilities/NekMesh/InputModules/InputGmsh.cpp +++ b/utilities/NekMesh/InputModules/InputGmsh.cpp @@ -166,10 +166,40 @@ std::vector triTensorNodeOrdering(const std::vector &nodes, int n) return nodeList; } +typedef boost::tuple 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 tetTensorNodeOrdering(const std::vector &nodes, int n) { std::vector nodeList; - int cnt2; int nTri = n*(n+1)/2; int nTet = n*(n+1)*(n+2)/6; @@ -193,35 +223,6 @@ std::vector tetTensorNodeOrdering(const std::vector &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 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 tmp; for (int k = 0, cnt = 0; k < n; ++k) @@ -344,9 +345,6 @@ std::vector tetTensorNodeOrdering(const std::vector &nodes, int n) std::vector prismTensorNodeOrdering(const std::vector &nodes, int n) { std::vector 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 prismTensorNodeOrdering(const std::vector &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 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; diff --git a/utilities/NekMesh/OutputModules/OutputGmsh.cpp b/utilities/NekMesh/OutputModules/OutputGmsh.cpp index b7d7f0b4d..c8050bce5 100644 --- a/utilities/NekMesh/OutputModules/OutputGmsh.cpp +++ b/utilities/NekMesh/OutputModules/OutputGmsh.cpp @@ -56,7 +56,6 @@ ModuleKey OutputGmsh::className = OutputGmsh::OutputGmsh(MeshSharedPtr m) : OutputModule(m) { map igelmap = InputGmsh::GenElmMap(); - map::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 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 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 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 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 hoQuad(faceIds, nodeList); hoQuad.Align(volFaceIds); -- GitLab From ecc4e7eb48f15320b5ce7fbcc2ebd3c3b264b627 Mon Sep 17 00:00:00 2001 From: David Moxey Date: Thu, 11 Aug 2016 17:33:40 +0100 Subject: [PATCH 14/25] Add some comments --- library/NekMeshUtils/MeshElements/Element.h | 2 +- utilities/NekMesh/InputModules/InputGmsh.cpp | 169 ++++++++++++++++++- 2 files changed, 161 insertions(+), 10 deletions(-) diff --git a/library/NekMeshUtils/MeshElements/Element.h b/library/NekMeshUtils/MeshElements/Element.h index ec32cada2..d1f6b4fb1 100644 --- a/library/NekMeshUtils/MeshElements/Element.h +++ b/library/NekMeshUtils/MeshElements/Element.h @@ -451,7 +451,7 @@ public: NEKMESHUTILS_EXPORT int GetMaxOrder(); - /// Complete this object. + /// Insert interior (i.e. volume) points into this element. NEKMESHUTILS_EXPORT virtual void MakeOrder( int order, SpatialDomains::GeometrySharedPtr geom, diff --git a/utilities/NekMesh/InputModules/InputGmsh.cpp b/utilities/NekMesh/InputModules/InputGmsh.cpp index 5f909662f..ceecce56a 100644 --- a/utilities/NekMesh/InputModules/InputGmsh.cpp +++ b/utilities/NekMesh/InputModules/InputGmsh.cpp @@ -64,8 +64,31 @@ ModuleKey InputGmsh::className = GetModuleFactory().RegisterCreatorFunction( std::map InputGmsh::elmMap = InputGmsh::GenElmMap(); /** - * @brief Reorder a quadrilateral to appear in Nektar++ ordering from - * Gmsh. + * @brief Reorder Gmsh nodes so that they appear in a tensor product format + * suitable for the interior of a Nektar++ quadrilateral. + * + * For an example, consider a second order quadrilateral. This routine will + * produce a permutation map that applies the following permutation: + * + * ``` + * 3---6---2 6---7---8 + * | | | | + * 7 8 5 ---> 3 4 5 + * | | | | + * 0---4---1 0---1---2 + * + * Gmsh tensor-product + * ``` + * + * We assume that Gmsh uses a recursive ordering system, so that interior nodes + * are reordered by calling this function recursively. + * + * @param nodes The integer IDs of the nodes to be reordered, in Gmsh format. + * @param n The number of nodes in one coordinate direction. If this is + * zero we assume no reordering needs to be done and return the + * identity permutation. + * + * @return The nodes vector in tensor-product ordering. */ std::vector quadTensorNodeOrdering(const std::vector &nodes, int n) { @@ -115,6 +138,35 @@ std::vector quadTensorNodeOrdering(const std::vector &nodes, int n) return nodeList; } +/** + * @brief Reorder Gmsh nodes so that they appear in a "tensor product" format + * suitable for the interior of a Nektar++ triangle. + * + * For an example, consider a third-order triangle. This routine will produce a + * permutation map that applies the following permutation: + * + * ``` + * 2 9 + * | \ | \ + * 7 6 7 7 + * | \ | \ + * 8 9 5 4 5 6 + * | \ | \ + * 0---3---4---1 0---1---2---3 + * + * Gmsh tensor-product + * ``` + * + * We assume that Gmsh uses a recursive ordering system, so that interior nodes + * are reordered by calling this function recursively. + * + * @param nodes The integer IDs of the nodes to be reordered, in Gmsh format. + * @param n The number of nodes in one coordinate direction. If this is + * zero we assume no reordering needs to be done and return the + * identity permutation. + * + * @return The nodes vector in tensor-product ordering. + */ std::vector triTensorNodeOrdering(const std::vector &nodes, int n) { std::vector nodeList; @@ -196,7 +248,40 @@ struct cmpop } }; - +/** + * @brief Reorder Gmsh nodes so that they appear in a "tensor product" format + * suitable for the interior of a Nektar++ tetrahedron. + * + * For an example, consider a second order tetrahedron. This routine will + * produce a permutation map that applies the following permutation: + * + * ``` + * 2 5 + * ,/|`\ ,/|`\ + * ,/ | `\ ,/ | `\ + * ,6 '. `5 ,3 '. `4 + * ,/ 8 `\ ,/ 8 `\ + * ,/ | `\ ,/ | `\ + * 0--------4--'.--------1 0--------1--'.--------2 + * `\. | ,/ `\. | ,/ + * `\. | ,9 `\. | ,7 + * `7. '. ,/ `6. '. ,/ + * `\. |/ `\. |/ + * `3 `9 + * + * Gmsh tensor-product + * ``` + * + * We assume that Gmsh uses a recursive ordering system, so that interior nodes + * are reordered by calling this function recursively. + * + * @param nodes The integer IDs of the nodes to be reordered, in Gmsh format. + * @param n The number of nodes in one coordinate direction. If this is + * zero we assume no reordering needs to be done and return the + * identity permutation. + * + * @return The nodes vector in tensor-product ordering. + */ std::vector tetTensorNodeOrdering(const std::vector &nodes, int n) { std::vector nodeList; @@ -342,6 +427,40 @@ std::vector tetTensorNodeOrdering(const std::vector &nodes, int n) return nodeList; } +/** + * @brief Reorder Gmsh nodes so that they appear in a "tensor product" format + * suitable for the interior of a Nektar++ prism. This routine is specifically + * designed for *interior* Gmsh nodes only. + * + * Prisms are a bit of a special case, in that interior nodes are heterogeneous + * in the number of points they use. As an example, for a second-order interior + * of a fourth-order prism, this routine calculates the mapping taking us + * + * ``` + * ,6 ,8 + * .-' | \ .-' | \ + * ,-3 | \ ,-5 | \ + * ,-' | \ | \ ,-' | \ | \ + * 2 | \7-------8 2 | \6-------7 + * | \ | ,-' \ ,-' | \ | ,-' \ ,-' + * | \,5------,0' | \,3------,4' + * |,-' \ ,-' |,-' \ ,-' + * 1-------4' 0-------1' + * + * Gmsh tensor-product + * ``` + * + * @todo Make this routine work for higher orders. The Gmsh ordering seems + * a little different than other elements, therefore the functionality is + * (for now) hard-coded to orders less than 5. + * + * @param nodes The integer IDs of the nodes to be reordered, in Gmsh format. + * @param n The number of nodes in one coordinate direction. If this is + * zero we assume no reordering needs to be done and return the + * identity permutation. + * + * @return The nodes vector in tensor-product ordering. + */ std::vector prismTensorNodeOrdering(const std::vector &nodes, int n) { std::vector nodeList; @@ -385,6 +504,39 @@ std::vector prismTensorNodeOrdering(const std::vector &nodes, int n) return nodeList; } +/** + * @brief Reorder Gmsh nodes so that they appear in a "tensor product" format + * suitable for the interior of a Nektar++ hexahedron. + * + * For an example, consider a second order hexahedron. This routine will produce + * a permutation map that applies the following permutation: + * + * ``` + * 3----13----2 6----7-----8 + * |\ |\ |\ |\ + * |15 24 | 14 |15 16 | 18 + * 9 \ 20 11 \ 3 \ 4 5 \ + * | 7----19+---6 | 24---25+---26 + * |22 | 26 | 23| |12 | 13 | 14| + * 0---+-8----1 | 0---+-1----2 | + * \ 17 25 \ 18 \ 21 22 \ 23 + * 10 | 21 12| 9 | 10 11| + * \| \| \| \| + * 4----16----5 18---19----20 + * + * Gmsh tensor-product + * ``` + * + * We assume that Gmsh uses a recursive ordering system, so that interior nodes + * are reordered by calling this function recursively. + * + * @param nodes The integer IDs of the nodes to be reordered, in Gmsh format. + * @param n The number of nodes in one coordinate direction. If this is + * zero we assume no reordering needs to be done and return the + * identity permutation. + * + * @return The nodes vector in tensor-product ordering. + */ std::vector hexTensorNodeOrdering(const std::vector &nodes, int n) { int i, j, k; @@ -555,12 +707,11 @@ InputGmsh::~InputGmsh() } /** - * Gmsh file contains a list of nodes and their coordinates, along with - * a list of elements and those nodes which define them. We read in and - * store the list of nodes in #m_node and store the list of elements in - * #m_element. Each new element is supplied with a list of entries from - * #m_node which defines the element. Finally some mesh statistics are - * printed. + * Gmsh file contains a list of nodes and their coordinates, along with a list + * of elements and those nodes which define them. We read in and store the list + * of nodes in #m_node and store the list of elements in #m_element. Each new + * element is supplied with a list of entries from m_node which defines the + * element. Finally some mesh statistics are printed. * * @param pFilename Filename of Gmsh file to read. */ -- GitLab From 416702587ca3d9e251eda14088ebc6eac25b78e7 Mon Sep 17 00:00:00 2001 From: David Moxey Date: Thu, 11 Aug 2016 17:42:46 +0100 Subject: [PATCH 15/25] Remove some debugging statements --- library/NekMeshUtils/MeshElements/Face.h | 3 --- library/NekMeshUtils/MeshElements/Prism.cpp | 4 ---- 2 files changed, 7 deletions(-) diff --git a/library/NekMeshUtils/MeshElements/Face.h b/library/NekMeshUtils/MeshElements/Face.h index 2e80743df..6cfa27bbf 100644 --- a/library/NekMeshUtils/MeshElements/Face.h +++ b/library/NekMeshUtils/MeshElements/Face.h @@ -265,9 +265,6 @@ public: m_faceNodes[cnt] = boost::shared_ptr( new Node(id++, x[0], x[1], x[2])); } - std::cout << "NFACEPTS = " << m_faceNodes.size() << std::endl; - - m_curveType = pType; } else if (m_vertexList.size() == 4) diff --git a/library/NekMeshUtils/MeshElements/Prism.cpp b/library/NekMeshUtils/MeshElements/Prism.cpp index 25943c094..21913e9f2 100644 --- a/library/NekMeshUtils/MeshElements/Prism.cpp +++ b/library/NekMeshUtils/MeshElements/Prism.cpp @@ -322,8 +322,6 @@ 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 xp(3); @@ -331,8 +329,6 @@ void Prism::MakeOrder(int order, xp[1] = py[i]; xp[2] = pz[i]; - cout << xp[0] << " " << xp[1] << " " << xp[2] << endl; - Array x(3, 0.0); for (int j = 0; j < coordDim; ++j) { -- GitLab From 453fc0414aea9ae2e5b2ffd45fc9f6a5585345a2 Mon Sep 17 00:00:00 2001 From: David Moxey Date: Mon, 15 Aug 2016 11:37:53 +0100 Subject: [PATCH 16/25] Start output of high-order boundary surfaces --- library/NekMeshUtils/MeshElements/Element.h | 3 +- library/NekMeshUtils/MeshElements/Face.h | 1 + .../NekMeshUtils/MeshElements/Hexahedron.cpp | 10 ++- .../NekMeshUtils/MeshElements/Hexahedron.h | 3 +- library/NekMeshUtils/MeshElements/Line.cpp | 61 +++++++++++++++++ library/NekMeshUtils/MeshElements/Line.h | 8 ++- library/NekMeshUtils/MeshElements/Mesh.cpp | 66 +++++++++++++++---- library/NekMeshUtils/MeshElements/Prism.cpp | 10 ++- library/NekMeshUtils/MeshElements/Prism.h | 4 +- .../MeshElements/Quadrilateral.cpp | 19 +++--- .../NekMeshUtils/MeshElements/Quadrilateral.h | 3 +- .../NekMeshUtils/MeshElements/Tetrahedron.cpp | 10 ++- .../NekMeshUtils/MeshElements/Tetrahedron.h | 3 +- .../NekMeshUtils/MeshElements/Triangle.cpp | 20 ++++-- library/NekMeshUtils/MeshElements/Triangle.h | 3 +- utilities/NekMesh/InputModules/InputGmsh.cpp | 16 ++--- .../NekMesh/OutputModules/OutputGmsh.cpp | 14 ++++ 17 files changed, 209 insertions(+), 45 deletions(-) diff --git a/library/NekMeshUtils/MeshElements/Element.h b/library/NekMeshUtils/MeshElements/Element.h index d1f6b4fb1..9c880a7f7 100644 --- a/library/NekMeshUtils/MeshElements/Element.h +++ b/library/NekMeshUtils/MeshElements/Element.h @@ -457,7 +457,8 @@ public: SpatialDomains::GeometrySharedPtr geom, LibUtilities::PointsType edgeType, int coordDim, - int &id) + int &id, + bool justConfig = false) { ASSERTL0(false, "This function should be implemented at a shape level."); diff --git a/library/NekMeshUtils/MeshElements/Face.h b/library/NekMeshUtils/MeshElements/Face.h index 6cfa27bbf..3e93e265d 100644 --- a/library/NekMeshUtils/MeshElements/Face.h +++ b/library/NekMeshUtils/MeshElements/Face.h @@ -249,6 +249,7 @@ public: const int nTriPts = nPoints * (nPoints + 1) / 2; const int nTriIntPts = (nPoints - 3) * (nPoints - 2) / 2; m_faceNodes.resize(nTriIntPts); + std::cout << "num tri int pts: " << nTriIntPts << " " << nTriPts << " " << nPoints << std::endl; for (int i = 3 + 3*(nPoints-2), cnt = 0; i < nTriPts; ++i, ++cnt) { diff --git a/library/NekMeshUtils/MeshElements/Hexahedron.cpp b/library/NekMeshUtils/MeshElements/Hexahedron.cpp index a307fa232..ea7af93c1 100644 --- a/library/NekMeshUtils/MeshElements/Hexahedron.cpp +++ b/library/NekMeshUtils/MeshElements/Hexahedron.cpp @@ -202,8 +202,10 @@ void Hexahedron::MakeOrder(int order, SpatialDomains::GeometrySharedPtr geom, LibUtilities::PointsType pType, int coordDim, - int &id) + int &id, + bool justConfig) { + m_curveType = pType; m_conf.m_order = order; m_volumeNodes.clear(); @@ -215,7 +217,11 @@ void Hexahedron::MakeOrder(int order, m_conf.m_faceNodes = true; m_conf.m_volumeNodes = true; - m_curveType = pType; + + if (justConfig) + { + return; + } int nPoints = order + 1; StdRegions::StdExpansionSharedPtr xmap = geom->GetXmap(); diff --git a/library/NekMeshUtils/MeshElements/Hexahedron.h b/library/NekMeshUtils/MeshElements/Hexahedron.h index 9b38e7588..fcd28ce80 100644 --- a/library/NekMeshUtils/MeshElements/Hexahedron.h +++ b/library/NekMeshUtils/MeshElements/Hexahedron.h @@ -84,7 +84,8 @@ public: SpatialDomains::GeometrySharedPtr geom, LibUtilities::PointsType pType, int coordDim, - int &id); + int &id, + bool justConfig = false); NEKMESHUTILS_EXPORT static unsigned int GetNumNodes(ElmtConfig pConf); NEKMESHUTILS_EXPORT virtual int GetFaceVertex(int i, int j) diff --git a/library/NekMeshUtils/MeshElements/Line.cpp b/library/NekMeshUtils/MeshElements/Line.cpp index 862838530..d8212761f 100644 --- a/library/NekMeshUtils/MeshElements/Line.cpp +++ b/library/NekMeshUtils/MeshElements/Line.cpp @@ -109,6 +109,67 @@ SpatialDomains::GeometrySharedPtr Line::GetGeom(int coordDim) return ret; } +void Line::MakeOrder(int order, + SpatialDomains::GeometrySharedPtr geom, + LibUtilities::PointsType pType, + int coordDim, + int &id, + bool justConfig) +{ + m_conf.m_order = order; + m_curveType = pType; + m_conf.m_volumeNodes = false; + m_volumeNodes.clear(); + + // Lines of order == 1 have no interior volume points. + if (order == 1) + { + m_conf.m_faceNodes = false; + return; + } + + m_conf.m_faceNodes = true; + + if (justConfig) + { + return; + } + + int nPoints = order + 1; + StdRegions::StdExpansionSharedPtr xmap = geom->GetXmap(); + + Array px; + LibUtilities::PointsKey pKey(nPoints, pType); + ASSERTL1(pKey.GetPointsDim() == 1, "Points distribution must be 1D"); + LibUtilities::PointsManager()[pKey]->GetPoints(px); + + Array > phys(coordDim); + + for (int i = 0; i < coordDim; ++i) + { + phys[i] = Array(xmap->GetTotPoints()); + xmap->BwdTrans(geom->GetCoeffs(i), phys[i]); + } + + int nQuadIntPts = (nPoints - 2) * (nPoints - 2); + m_volumeNodes.resize(nQuadIntPts); + + for (int i = 1, cnt = 0; i < nPoints-1; ++i) + { + Array xp(1); + xp[0] = px[i]; + + Array x(3, 0.0); + for (int k = 0; k < coordDim; ++k) + { + x[k] = xmap->PhysEvaluate(xp, phys[k]); + } + + m_volumeNodes[cnt] = boost::shared_ptr( + new Node(id++, x[0], x[1], x[2])); + } +} + /** * @brief Return the number of nodes defining a line. */ diff --git a/library/NekMeshUtils/MeshElements/Line.h b/library/NekMeshUtils/MeshElements/Line.h index 158475be6..b267d4b28 100644 --- a/library/NekMeshUtils/MeshElements/Line.h +++ b/library/NekMeshUtils/MeshElements/Line.h @@ -68,7 +68,13 @@ public: NEKMESHUTILS_EXPORT virtual SpatialDomains::GeometrySharedPtr GetGeom( int coordDim); - + NEKMESHUTILS_EXPORT virtual void MakeOrder( + int order, + SpatialDomains::GeometrySharedPtr geom, + LibUtilities::PointsType pType, + int coordDim, + int &id, + bool justConfig = false); NEKMESHUTILS_EXPORT static unsigned int GetNumNodes(ElmtConfig pConf); }; } diff --git a/library/NekMeshUtils/MeshElements/Mesh.cpp b/library/NekMeshUtils/MeshElements/Mesh.cpp index ca5e94fdb..dae08817c 100644 --- a/library/NekMeshUtils/MeshElements/Mesh.cpp +++ b/library/NekMeshUtils/MeshElements/Mesh.cpp @@ -96,27 +96,33 @@ void Mesh::MakeOrder(int order, boost::unordered_map faceGeoms; boost::unordered_map volGeoms; - // Decide on distribution of points to use for each shape type. + // Decide on distribution of points to use for each shape type based on the + // input we've been supplied. std::map pTypes; if (distType == LibUtilities::ePolyEvenlySpaced) { - pTypes[LibUtilities::eSegment] = LibUtilities::ePolyEvenlySpaced; - pTypes[LibUtilities::eTriangle] = LibUtilities::eNodalTriEvenlySpaced; + 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::eTetrahedron] = + LibUtilities::eNodalTetEvenlySpaced; + pTypes[LibUtilities::ePrism] = LibUtilities::eNodalPrismEvenlySpaced; + pTypes[LibUtilities::eHexahedron] = LibUtilities::ePolyEvenlySpaced; } else if (distType == LibUtilities::eGaussLobattoLegendre) { - 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; + pTypes[LibUtilities::eSegment] = LibUtilities::eGaussLobattoLegendre; + pTypes[LibUtilities::eTriangle] = LibUtilities::eNodalTriElec; + pTypes[LibUtilities::eQuadrilateral] = + LibUtilities::eGaussLobattoLegendre; + pTypes[LibUtilities::eTetrahedron] = LibUtilities::eNodalTetElec; + pTypes[LibUtilities::eHexahedron] = LibUtilities::eGaussLobattoLegendre; } + // Begin by generating Nektar++ geometry objects for edges, faces and + // elements so that we don't affect any neighbouring elements in the mesh as + // we process each element. for(eit = m_edgeSet.begin(); eit != m_edgeSet.end(); eit++) { SpatialDomains::Geometry1DSharedPtr geom = @@ -144,6 +150,8 @@ void Mesh::MakeOrder(int order, boost::unordered_set processedEdges, processedFaces, processedVolumes; + // Call MakeOrder with our generated geometries on each edge to fill in edge + // interior nodes. for(eit = m_edgeSet.begin(); eit != m_edgeSet.end(); eit++) { int edgeId = (*eit)->m_id; @@ -158,6 +166,8 @@ void Mesh::MakeOrder(int order, processedEdges.insert(edgeId); } + // Call MakeOrder with our generated geometries on each face to fill in face + // interior nodes. for(fit = m_faceSet.begin(); fit != m_faceSet.end(); fit++) { int faceId = (*fit)->m_id; @@ -174,6 +184,40 @@ void Mesh::MakeOrder(int order, processedFaces.insert(faceId); } + // Copy curvature into boundary conditions + for (int i = 0; i < m_element[1].size(); ++i) + { + ElementSharedPtr el = m_element[1][i]; + EdgeSharedPtr edge = el->GetEdgeLink(); + + if (!edge) + { + continue; + } + + // Copy face curvature + el->SetVolumeNodes(edge->m_edgeNodes); + el->MakeOrder(order, SpatialDomains::GeometrySharedPtr(), + pTypes[el->GetConf().m_e], m_spaceDim, id, true); + } + + for (int i = 0; i < m_element[2].size(); ++i) + { + ElementSharedPtr el = m_element[2][i]; + FaceSharedPtr face = el->GetFaceLink(); + + if (!face) + { + continue; + } + + // Copy face curvature + el->MakeOrder(order, SpatialDomains::GeometrySharedPtr(), + pTypes[el->GetConf().m_e], m_spaceDim, id, true); + el->SetVolumeNodes(face->m_faceNodes); + } + + // Finally, fill in volumes. const int nElmt = m_element[m_expDim].size(); for (int i = 0; i < nElmt; ++i) { diff --git a/library/NekMeshUtils/MeshElements/Prism.cpp b/library/NekMeshUtils/MeshElements/Prism.cpp index 21913e9f2..5d4e4097b 100644 --- a/library/NekMeshUtils/MeshElements/Prism.cpp +++ b/library/NekMeshUtils/MeshElements/Prism.cpp @@ -281,9 +281,11 @@ void Prism::MakeOrder(int order, SpatialDomains::GeometrySharedPtr geom, LibUtilities::PointsType pType, int coordDim, - int &id) + int &id, + bool justConfig) { m_conf.m_order = order; + m_curveType = pType; m_volumeNodes.clear(); if (order == 1) @@ -300,7 +302,11 @@ void Prism::MakeOrder(int order, m_conf.m_faceNodes = true; m_conf.m_volumeNodes = true; - m_curveType = pType; + + if (justConfig) + { + return; + } int nPoints = order + 1; StdRegions::StdExpansionSharedPtr xmap = geom->GetXmap(); diff --git a/library/NekMeshUtils/MeshElements/Prism.h b/library/NekMeshUtils/MeshElements/Prism.h index 65b4c506b..c0862f2bc 100644 --- a/library/NekMeshUtils/MeshElements/Prism.h +++ b/library/NekMeshUtils/MeshElements/Prism.h @@ -84,7 +84,9 @@ public: SpatialDomains::GeometrySharedPtr geom, LibUtilities::PointsType pType, int coordDim, - int &id); + int &id, + bool justConfig = false); + NEKMESHUTILS_EXPORT static unsigned int GetNumNodes(ElmtConfig pConf); NEKMESHUTILS_EXPORT virtual int GetFaceVertex(int i, int j) { diff --git a/library/NekMeshUtils/MeshElements/Quadrilateral.cpp b/library/NekMeshUtils/MeshElements/Quadrilateral.cpp index 84f74963e..6f6579378 100644 --- a/library/NekMeshUtils/MeshElements/Quadrilateral.cpp +++ b/library/NekMeshUtils/MeshElements/Quadrilateral.cpp @@ -138,24 +138,25 @@ void Quadrilateral::MakeOrder(int order, SpatialDomains::GeometrySharedPtr geom, LibUtilities::PointsType pType, int coordDim, - int &id) + int &id, + bool justConfig) { - m_conf.m_order = order; + m_conf.m_order = order; + m_curveType = pType; + m_conf.m_volumeNodes = false; + m_volumeNodes.clear(); + // Quadrilaterals of order == 1 have no interior volume points. if (order == 1) { - m_conf.m_volumeNodes = m_conf.m_faceNodes = false; + m_conf.m_faceNodes = false; return; } - m_conf.m_faceNodes = true; - m_conf.m_volumeNodes = false; - m_curveType = pType; + m_conf.m_faceNodes = true; - // Quadrilaterals of order < 2 have no interior volume points. - if (order < 2) + if (justConfig) { - m_volumeNodes.clear(); return; } diff --git a/library/NekMeshUtils/MeshElements/Quadrilateral.h b/library/NekMeshUtils/MeshElements/Quadrilateral.h index 3819ad323..620a228cc 100644 --- a/library/NekMeshUtils/MeshElements/Quadrilateral.h +++ b/library/NekMeshUtils/MeshElements/Quadrilateral.h @@ -84,7 +84,8 @@ public: SpatialDomains::GeometrySharedPtr geom, LibUtilities::PointsType pType, int coordDim, - int &id); + int &id, + bool justConfig = false); NEKMESHUTILS_EXPORT static unsigned int GetNumNodes(ElmtConfig pConf); }; diff --git a/library/NekMeshUtils/MeshElements/Tetrahedron.cpp b/library/NekMeshUtils/MeshElements/Tetrahedron.cpp index 895b4422b..ea7533038 100644 --- a/library/NekMeshUtils/MeshElements/Tetrahedron.cpp +++ b/library/NekMeshUtils/MeshElements/Tetrahedron.cpp @@ -257,9 +257,11 @@ void Tetrahedron::MakeOrder(int order, SpatialDomains::GeometrySharedPtr geom, LibUtilities::PointsType pType, int coordDim, - int &id) + int &id, + bool justConfig) { m_conf.m_order = order; + m_curveType = pType; m_volumeNodes.clear(); if (order == 1 || order == 2) @@ -279,10 +281,14 @@ void Tetrahedron::MakeOrder(int order, return; } - m_curveType = pType; m_conf.m_faceNodes = true; m_conf.m_volumeNodes = true; + if (justConfig) + { + return; + } + int nPoints = order + 1; StdRegions::StdExpansionSharedPtr xmap = geom->GetXmap(); diff --git a/library/NekMeshUtils/MeshElements/Tetrahedron.h b/library/NekMeshUtils/MeshElements/Tetrahedron.h index d8d1e4cf0..bf915f731 100644 --- a/library/NekMeshUtils/MeshElements/Tetrahedron.h +++ b/library/NekMeshUtils/MeshElements/Tetrahedron.h @@ -79,7 +79,8 @@ public: SpatialDomains::GeometrySharedPtr geom, LibUtilities::PointsType pType, int coordDim, - int &id); + int &id, + bool justConfig = false); NEKMESHUTILS_EXPORT static unsigned int GetNumNodes(ElmtConfig pConf); NEKMESHUTILS_EXPORT virtual int GetFaceVertex(int i, int j) diff --git a/library/NekMeshUtils/MeshElements/Triangle.cpp b/library/NekMeshUtils/MeshElements/Triangle.cpp index e465a45a4..7b2283cec 100644 --- a/library/NekMeshUtils/MeshElements/Triangle.cpp +++ b/library/NekMeshUtils/MeshElements/Triangle.cpp @@ -173,12 +173,25 @@ void Triangle::MakeOrder(int order, SpatialDomains::GeometrySharedPtr geom, LibUtilities::PointsType pType, int coordDim, - int &id) + int &id, + bool justConfig) { + m_conf.m_order = order; + m_curveType = pType; + m_conf.m_volumeNodes = false; + m_volumeNodes.clear(); + // Triangles of order < 3 have no interior volume points. - if (order < 3) + if (order == 1 || order == 2) + { + m_conf.m_faceNodes = false; + return; + } + + m_conf.m_faceNodes = true; + + if (justConfig) { - m_volumeNodes.clear(); return; } @@ -218,7 +231,6 @@ void Triangle::MakeOrder(int order, 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 = false; diff --git a/library/NekMeshUtils/MeshElements/Triangle.h b/library/NekMeshUtils/MeshElements/Triangle.h index b8f0009c7..33b32aab0 100644 --- a/library/NekMeshUtils/MeshElements/Triangle.h +++ b/library/NekMeshUtils/MeshElements/Triangle.h @@ -214,7 +214,8 @@ public: SpatialDomains::GeometrySharedPtr geom, LibUtilities::PointsType pType, int coordDim, - int &id); + int &id, + bool justConfig = false); NEKMESHUTILS_EXPORT static unsigned int GetNumNodes(ElmtConfig pConf); }; diff --git a/utilities/NekMesh/InputModules/InputGmsh.cpp b/utilities/NekMesh/InputModules/InputGmsh.cpp index ceecce56a..d0003f517 100644 --- a/utilities/NekMesh/InputModules/InputGmsh.cpp +++ b/utilities/NekMesh/InputModules/InputGmsh.cpp @@ -1679,15 +1679,15 @@ std::map InputGmsh::GenElmMap() tmp[ 4] = ElmtConfig(eTetrahedron, 1, false, false); tmp[ 5] = ElmtConfig(eHexahedron, 1, false, false); tmp[ 6] = ElmtConfig(ePrism, 1, false, false); - tmp[ 7] = ElmtConfig(ePyramid, 1, true, true); - tmp[ 8] = ElmtConfig(eSegment, 2, true, true); - tmp[ 9] = ElmtConfig(eTriangle, 2, true, true); - tmp[ 10] = ElmtConfig(eQuadrilateral, 2, true, true); - tmp[ 11] = ElmtConfig(eTetrahedron, 2, true, true); + tmp[ 7] = ElmtConfig(ePyramid, 1, false, false); + tmp[ 8] = ElmtConfig(eSegment, 2, true, false); + tmp[ 9] = ElmtConfig(eTriangle, 2, false, false); + tmp[ 10] = ElmtConfig(eQuadrilateral, 2, true, false); + tmp[ 11] = ElmtConfig(eTetrahedron, 2, false, false); tmp[ 12] = ElmtConfig(eHexahedron, 2, true, true); - tmp[ 13] = ElmtConfig(ePrism, 2, true, true); - tmp[ 14] = ElmtConfig(ePyramid, 2, true, true); - tmp[ 15] = ElmtConfig(ePoint, 1, true, false); + tmp[ 13] = ElmtConfig(ePrism, 2, true, false); + tmp[ 14] = ElmtConfig(ePyramid, 2, true, false); + tmp[ 15] = ElmtConfig(ePoint, 1, false, false); tmp[ 16] = ElmtConfig(eQuadrilateral, 2, false, false); tmp[ 17] = ElmtConfig(eHexahedron, 2, false, false); tmp[ 18] = ElmtConfig(ePrism, 2, false, false); diff --git a/utilities/NekMesh/OutputModules/OutputGmsh.cpp b/utilities/NekMesh/OutputModules/OutputGmsh.cpp index c8050bce5..bcf778303 100644 --- a/utilities/NekMesh/OutputModules/OutputGmsh.cpp +++ b/utilities/NekMesh/OutputModules/OutputGmsh.cpp @@ -212,11 +212,18 @@ void OutputGmsh::Process() tags.push_back(nodeList[j]->m_id); } + cout << "GOT " << nodeList.size() << " VERTS" << endl; + cout << "GOT " << edgeList.size() << " EDGES" << endl; + cout << "GOT " << faceList.size() << " FACES" << endl; + cout << "GOT " << volList.size() << " VOLS" << endl; + // Process edge-interior points for (int j = 0; j < edgeList.size(); ++j) { nodeList = edgeList[j]->m_edgeNodes; + cout << "edge " << j << ": " << edgeList[j]->m_edgeNodes.size() << endl; + if (e->GetEdgeOrient(j, edgeList[j]) == StdRegions::eForwards) { for (int k = 0; k < nodeList.size(); ++k) @@ -269,6 +276,7 @@ void OutputGmsh::Process() } // Process volume nodes + cout << volList.size() << endl; for (int j = 0; j < volList.size(); ++j) { tags.push_back(volList[j]->m_id); @@ -278,6 +286,12 @@ void OutputGmsh::Process() vector reordering = InputGmsh::CreateReordering(elmtType); vector inv(tags.size()); + cout << "TAGS LENGTH: " << tags.size() << endl; + cout << "REORDERING LENGTH: " << reordering.size() << endl; + + ASSERTL1(tags.size() == reordering.size(), + "Reordering map size not equal to element tags."); + for (int j = 0; j < tags.size(); ++j) { inv[reordering[j]] = j; -- GitLab From fb563319ed99ed4631ddbc7002c0204a624ceef1 Mon Sep 17 00:00:00 2001 From: David Moxey Date: Mon, 15 Aug 2016 14:38:11 +0100 Subject: [PATCH 17/25] Working high-order boundary information --- library/NekMeshUtils/MeshElements/Element.cpp | 16 ++++++- library/NekMeshUtils/MeshElements/Element.h | 4 +- library/NekMeshUtils/MeshElements/Face.h | 1 - utilities/NekMesh/Module.cpp | 47 ++++++++++++++----- .../NekMesh/OutputModules/OutputGmsh.cpp | 11 ----- 5 files changed, 51 insertions(+), 28 deletions(-) diff --git a/library/NekMeshUtils/MeshElements/Element.cpp b/library/NekMeshUtils/MeshElements/Element.cpp index a2e44d1d8..9cd7bb5a3 100644 --- a/library/NekMeshUtils/MeshElements/Element.cpp +++ b/library/NekMeshUtils/MeshElements/Element.cpp @@ -74,10 +74,16 @@ Element::Element(ElmtConfig pConf, unsigned int pNumNodes, * @param p Index of the vertex to replace. * @param pNew New vertex. */ -void Element::SetVertex(unsigned int p, NodeSharedPtr pNew) +void Element::SetVertex(unsigned int p, NodeSharedPtr pNew, bool descend) { NodeSharedPtr vOld = m_vertex[p]; m_vertex[p] = pNew; + + if (!descend) + { + return; + } + for (unsigned int i = 0; i < m_edge.size(); ++i) { if (m_edge[i]->m_n1 == vOld) @@ -122,10 +128,16 @@ void Element::SetVertex(unsigned int p, NodeSharedPtr pNew) * @param p Index of the edge to replace. * @param pNew New edge. */ -void Element::SetEdge(unsigned int p, EdgeSharedPtr pNew) +void Element::SetEdge(unsigned int p, EdgeSharedPtr pNew, bool descend) { EdgeSharedPtr vOld = m_edge[p]; m_edge[p] = pNew; + + if (!descend) + { + return; + } + for (unsigned int i = 0; i < m_face.size(); ++i) { for (unsigned int j = 0; j < m_face[i]->m_edgeList.size(); ++j) diff --git a/library/NekMeshUtils/MeshElements/Element.h b/library/NekMeshUtils/MeshElements/Element.h index 9c880a7f7..4cc60cb98 100644 --- a/library/NekMeshUtils/MeshElements/Element.h +++ b/library/NekMeshUtils/MeshElements/Element.h @@ -244,9 +244,9 @@ public: m_id = p; } /// Replace a vertex with another vertex object. - NEKMESHUTILS_EXPORT void SetVertex(unsigned int p, NodeSharedPtr pNew); + NEKMESHUTILS_EXPORT void SetVertex(unsigned int p, NodeSharedPtr pNew, bool descend = true); /// Replace an edge with another edge object. - NEKMESHUTILS_EXPORT void SetEdge(unsigned int p, EdgeSharedPtr pNew); + NEKMESHUTILS_EXPORT void SetEdge(unsigned int p, EdgeSharedPtr pNew, bool descend = true); /// Replace a face with another face object. NEKMESHUTILS_EXPORT void SetFace(unsigned int p, FaceSharedPtr pNew); /// Set a correspondence between this element and an edge diff --git a/library/NekMeshUtils/MeshElements/Face.h b/library/NekMeshUtils/MeshElements/Face.h index 3e93e265d..6cfa27bbf 100644 --- a/library/NekMeshUtils/MeshElements/Face.h +++ b/library/NekMeshUtils/MeshElements/Face.h @@ -249,7 +249,6 @@ public: const int nTriPts = nPoints * (nPoints + 1) / 2; const int nTriIntPts = (nPoints - 3) * (nPoints - 2) / 2; m_faceNodes.resize(nTriIntPts); - std::cout << "num tri int pts: " << nTriIntPts << " " << nTriPts << " " << nPoints << std::endl; for (int i = 3 + 3*(nPoints-2), cnt = 0; i < nTriPts; ++i, ++cnt) { diff --git a/utilities/NekMesh/Module.cpp b/utilities/NekMesh/Module.cpp index 35f6f7a73..191026c66 100644 --- a/utilities/NekMesh/Module.cpp +++ b/utilities/NekMesh/Module.cpp @@ -214,12 +214,13 @@ void Module::ProcessEdges(bool ReprocessEdges) // Create links for 1D elements for (int i = 0; i < m_mesh->m_element[1].size(); ++i) { - NodeSharedPtr v0 = m_mesh->m_element[1][i]->GetVertex(0); - NodeSharedPtr v1 = m_mesh->m_element[1][i]->GetVertex(1); + ElementSharedPtr elmt = m_mesh->m_element[1][i]; + NodeSharedPtr v0 = elmt->GetVertex(0); + NodeSharedPtr v1 = elmt->GetVertex(1); vector edgeNodes; EdgeSharedPtr E = boost::shared_ptr( - new Edge(v0, v1, edgeNodes, - m_mesh->m_element[1][i]->GetConf().m_edgeCurveType)); + new Edge(v0, v1, edgeNodes, elmt->GetConf().m_edgeCurveType)); + EdgeSet::iterator it = m_mesh->m_edgeSet.find(E); if (it == m_mesh->m_edgeSet.end()) { @@ -227,20 +228,25 @@ void Module::ProcessEdges(bool ReprocessEdges) << "1D element " << i << endl; abort(); } - m_mesh->m_element[1][i]->SetEdgeLink(*it); + elmt->SetEdgeLink(*it); // Update 2D element boundary map. pair eMap = (*it)->m_elLink.at(0); eMap.first->SetBoundaryLink(eMap.second, i); + // Update vertices + elmt->SetVertex(0, (*it)->m_n1, false); + elmt->SetVertex(1, (*it)->m_n2, false); + // Copy curvature to edge. if ((*it)->m_edgeNodes.size() > 0) { - ElementSharedPtr edge = m_mesh->m_element[1][i]; + ElementSharedPtr edge = elmt; if (edge->GetVertex(0) == (*it)->m_n1) { edge->SetVolumeNodes((*it)->m_edgeNodes); } + elmt->SetCurveType((*it)->m_curveType); } } } @@ -283,14 +289,14 @@ void Module::ProcessFaces(bool ReprocessFaces) { (*(testIns.first))->m_id = fid++; (*(testIns.first))->m_elLink.push_back( - pair(elmt[i],j)); + pair(elmt[i],j)); } else { elmt[i]->SetFace(j,*testIns.first); // Update face to element map. (*(testIns.first))->m_elLink.push_back( - pair(elmt[i],j)); + pair(elmt[i],j)); } } } @@ -299,12 +305,14 @@ void Module::ProcessFaces(bool ReprocessFaces) // Create links for 2D elements for (int i = 0; i < m_mesh->m_element[2].size(); ++i) { - vector vertices = m_mesh->m_element[2][i]->GetVertexList(); + ElementSharedPtr elmt = m_mesh->m_element[2][i]; + vector vertices = elmt->GetVertexList(); vector faceNodes; - vector edgeList = m_mesh->m_element[2][i]->GetEdgeList(); + vector edgeList = elmt->GetEdgeList(); FaceSharedPtr F = boost::shared_ptr( new Face(vertices, faceNodes, edgeList, - m_mesh->m_element[2][i]->GetConf().m_faceCurveType)); + elmt->GetConf().m_faceCurveType)); + FaceSet::iterator it = m_mesh->m_faceSet.find(F); if (it == m_mesh->m_faceSet.end()) { @@ -312,11 +320,26 @@ void Module::ProcessFaces(bool ReprocessFaces) << "element " << i << endl; abort(); } - m_mesh->m_element[2][i]->SetFaceLink(*it); + + elmt->SetFaceLink(*it); + + // Set edges/vertices + for (int j = 0; j < elmt->GetVertexCount(); ++j) + { + elmt->SetVertex(j, (*it)->m_vertexList[j], false); + elmt->SetEdge(j, (*it)->m_edgeList[j], false); + } // Update 3D element boundary map. pair eMap = (*it)->m_elLink.at(0); eMap.first->SetBoundaryLink(eMap.second, i); + + // Copy face curvature + if ((*it)->m_faceNodes.size() > 0) + { + elmt->SetVolumeNodes((*it)->m_faceNodes); + elmt->SetCurveType((*it)->m_curveType); + } } } diff --git a/utilities/NekMesh/OutputModules/OutputGmsh.cpp b/utilities/NekMesh/OutputModules/OutputGmsh.cpp index bcf778303..080bef898 100644 --- a/utilities/NekMesh/OutputModules/OutputGmsh.cpp +++ b/utilities/NekMesh/OutputModules/OutputGmsh.cpp @@ -212,18 +212,11 @@ void OutputGmsh::Process() tags.push_back(nodeList[j]->m_id); } - cout << "GOT " << nodeList.size() << " VERTS" << endl; - cout << "GOT " << edgeList.size() << " EDGES" << endl; - cout << "GOT " << faceList.size() << " FACES" << endl; - cout << "GOT " << volList.size() << " VOLS" << endl; - // Process edge-interior points for (int j = 0; j < edgeList.size(); ++j) { nodeList = edgeList[j]->m_edgeNodes; - cout << "edge " << j << ": " << edgeList[j]->m_edgeNodes.size() << endl; - if (e->GetEdgeOrient(j, edgeList[j]) == StdRegions::eForwards) { for (int k = 0; k < nodeList.size(); ++k) @@ -276,7 +269,6 @@ void OutputGmsh::Process() } // Process volume nodes - cout << volList.size() << endl; for (int j = 0; j < volList.size(); ++j) { tags.push_back(volList[j]->m_id); @@ -286,9 +278,6 @@ void OutputGmsh::Process() vector reordering = InputGmsh::CreateReordering(elmtType); vector inv(tags.size()); - cout << "TAGS LENGTH: " << tags.size() << endl; - cout << "REORDERING LENGTH: " << reordering.size() << endl; - ASSERTL1(tags.size() == reordering.size(), "Reordering map size not equal to element tags."); -- GitLab From 8232ec4692329790a592b62024993090a70e22d0 Mon Sep 17 00:00:00 2001 From: David Moxey Date: Mon, 15 Aug 2016 19:31:30 +0100 Subject: [PATCH 18/25] Start cleaning up files --- .../NekMesh/OutputModules/OutputGmsh.cpp | 38 ++++++++++++++----- .../NekMesh/OutputModules/OutputNekpp.cpp | 8 ++++ 2 files changed, 36 insertions(+), 10 deletions(-) diff --git a/utilities/NekMesh/OutputModules/OutputGmsh.cpp b/utilities/NekMesh/OutputModules/OutputGmsh.cpp index 080bef898..d347e6362 100644 --- a/utilities/NekMesh/OutputModules/OutputGmsh.cpp +++ b/utilities/NekMesh/OutputModules/OutputGmsh.cpp @@ -64,6 +64,8 @@ OutputGmsh::OutputGmsh(MeshSharedPtr m) : OutputModule(m) { elmMap[it->second] = it->first; } + + m_config["order"] = ConfigOption(false, "-1", "Enforce a polynomial order"); } OutputGmsh::~OutputGmsh() @@ -87,6 +89,9 @@ void OutputGmsh::Process() cout << "OutputGmsh: Writing file..." << endl; } + boost::unordered_map > orderingMap; + boost::unordered_map >::iterator oIt; + // Open the file stream. OpenStream(); @@ -98,23 +103,36 @@ void OutputGmsh::Process() int id = m_mesh->m_vertexSet.size(); vector toComplete; - int maxOrder = -1; + int order = m_config["order"].as(); - // Do first pass over elements of expansion dimension to determine - // which elements need completion. - for (int i = 0; i < m_mesh->m_element[m_mesh->m_expDim].size(); ++i) + if (order != -1) + { + if (m_mesh->m_verbose) + { + cout << "Making mesh of order " << order << endl; + } + } + else { - ElementSharedPtr e = m_mesh->m_element[m_mesh->m_expDim][i]; - if (e->GetMaxOrder() > maxOrder) + // Do first pass over elements of expansion dimension to determine + // which elements need completion. + for (int i = 0; i < m_mesh->m_element[m_mesh->m_expDim].size(); ++i) { - maxOrder = e->GetMaxOrder(); + ElementSharedPtr e = m_mesh->m_element[m_mesh->m_expDim][i]; + if (e->GetMaxOrder() > order) + { + order = e->GetMaxOrder(); + } } } // Convert this mesh into a high-order mesh of uniform order. - cout << "Mesh order of " << maxOrder << " detected" << endl; - maxOrder = 4; - m_mesh->MakeOrder(maxOrder, LibUtilities::ePolyEvenlySpaced); + if (m_mesh->m_verbose) + { + cout << "Mesh order of " << order << " detected" << endl; + } + + m_mesh->MakeOrder(order, LibUtilities::ePolyEvenlySpaced); // Add edge- and face-interior nodes to vertex set. EdgeSet::iterator eIt; diff --git a/utilities/NekMesh/OutputModules/OutputNekpp.cpp b/utilities/NekMesh/OutputModules/OutputNekpp.cpp index 43f9e89a8..6735cb9e9 100644 --- a/utilities/NekMesh/OutputModules/OutputNekpp.cpp +++ b/utilities/NekMesh/OutputModules/OutputNekpp.cpp @@ -68,6 +68,7 @@ OutputNekpp::OutputNekpp(MeshSharedPtr m) : OutputModule(m) m_config["test"] = ConfigOption( true, "0", "Attempt to load resulting mesh and create meshgraph."); m_config["uncompress"] = ConfigOption(true, "0", "Uncompress xml sections"); + m_config["order"] = ConfigOption(false, "-1", "Enforce a polynomial order"); } OutputNekpp::~OutputNekpp() @@ -98,6 +99,13 @@ void OutputNekpp::Process() cout << "OutputNekpp: Writing file..." << endl; } + int order = m_config["order"].as(); + + if (order != -1) + { + m_mesh->MakeOrder(order, LibUtilities::ePolyEvenlySpaced); + } + TiXmlDocument doc; TiXmlDeclaration *decl = new TiXmlDeclaration("1.0", "utf-8", ""); doc.LinkEndChild(decl); -- GitLab From b78ae33d152edd7a0492cbe66196115211c7c78f Mon Sep 17 00:00:00 2001 From: David Moxey Date: Mon, 15 Aug 2016 19:36:09 +0100 Subject: [PATCH 19/25] Add cache for element output --- .../NekMesh/OutputModules/OutputGmsh.cpp | 30 +++++++++++++------ 1 file changed, 21 insertions(+), 9 deletions(-) diff --git a/utilities/NekMesh/OutputModules/OutputGmsh.cpp b/utilities/NekMesh/OutputModules/OutputGmsh.cpp index d347e6362..77c587b05 100644 --- a/utilities/NekMesh/OutputModules/OutputGmsh.cpp +++ b/utilities/NekMesh/OutputModules/OutputGmsh.cpp @@ -292,22 +292,34 @@ void OutputGmsh::Process() tags.push_back(volList[j]->m_id); } - // Construct inverse of input reordering - vector reordering = InputGmsh::CreateReordering(elmtType); - vector inv(tags.size()); + // Construct inverse of input reordering. First try to find it in + // our cache. + oIt = orderingMap.find(elmtType); - ASSERTL1(tags.size() == reordering.size(), - "Reordering map size not equal to element tags."); - - for (int j = 0; j < tags.size(); ++j) + // If it's not created, then create it. + if (oIt == orderingMap.end()) { - inv[reordering[j]] = j; + vector reordering = InputGmsh::CreateReordering(elmtType); + vector inv(tags.size()); + + ASSERTL1(tags.size() == reordering.size(), + "Reordering map size not equal to element tags."); + + for (int j = 0; j < tags.size(); ++j) + { + inv[reordering[j]] = j; + } + + oIt = orderingMap.insert(make_pair(elmtType, inv)).first; } + ASSERTL1(tags.size() == oIt->second.size(), + "Reordering map size not equal to element tags."); + // Finally write element nodes. for (int j = 0; j < tags.size(); ++j) { - m_mshFile << tags[inv[j]] + 1 << " "; + m_mshFile << tags[oIt->second[j]] + 1 << " "; } m_mshFile << endl; -- GitLab From 393eda70b5c4ff3b7d950c3ace96bfd437104d4a Mon Sep 17 00:00:00 2001 From: David Moxey Date: Mon, 15 Aug 2016 22:20:04 +0100 Subject: [PATCH 20/25] Add some tests --- utilities/NekMesh/CMakeLists.txt | 5 + .../NekMesh/Tests/Gmsh/CubeHexLinear.msh | 68 ++++ .../NekMesh/Tests/Gmsh/CubeHexLinear.tst | 19 + .../NekMesh/Tests/Gmsh/CubePrismLinear.msh | 84 +++++ .../NekMesh/Tests/Gmsh/CubePrismLinear.tst | 19 + .../NekMesh/Tests/Gmsh/CubeTetLinear.msh | 343 ++++++++++++++++++ .../NekMesh/Tests/Gmsh/CubeTetLinear.tst | 19 + .../NekMesh/Tests/Gmsh/SquareQuadLinear.msh | 66 ++++ .../NekMesh/Tests/Gmsh/SquareQuadLinear.tst | 19 + .../NekMesh/Tests/Gmsh/SquareTriLinear.msh | 82 +++++ .../NekMesh/Tests/Gmsh/SquareTriLinear.tst | 19 + 11 files changed, 743 insertions(+) create mode 100644 utilities/NekMesh/Tests/Gmsh/CubeHexLinear.msh create mode 100644 utilities/NekMesh/Tests/Gmsh/CubeHexLinear.tst create mode 100644 utilities/NekMesh/Tests/Gmsh/CubePrismLinear.msh create mode 100644 utilities/NekMesh/Tests/Gmsh/CubePrismLinear.tst create mode 100644 utilities/NekMesh/Tests/Gmsh/CubeTetLinear.msh create mode 100644 utilities/NekMesh/Tests/Gmsh/CubeTetLinear.tst create mode 100644 utilities/NekMesh/Tests/Gmsh/SquareQuadLinear.msh create mode 100644 utilities/NekMesh/Tests/Gmsh/SquareQuadLinear.tst create mode 100644 utilities/NekMesh/Tests/Gmsh/SquareTriLinear.msh create mode 100644 utilities/NekMesh/Tests/Gmsh/SquareTriLinear.tst diff --git a/utilities/NekMesh/CMakeLists.txt b/utilities/NekMesh/CMakeLists.txt index a1ff2828c..704f207a2 100644 --- a/utilities/NekMesh/CMakeLists.txt +++ b/utilities/NekMesh/CMakeLists.txt @@ -105,15 +105,20 @@ ADD_NEKTAR_TEST (Nektar++/CylQuadBl) # Gmsh tests ADD_NEKTAR_TEST (Gmsh/CubeAllElements) ADD_NEKTAR_TEST (Gmsh/CubeHex) +ADD_NEKTAR_TEST (Gmsh/CubeHexLinear) ADD_NEKTAR_TEST (Gmsh/CubePrism) +ADD_NEKTAR_TEST (Gmsh/CubePrismLinear) ADD_NEKTAR_TEST (Gmsh/CubeTet) +ADD_NEKTAR_TEST (Gmsh/CubeTetLinear) IF (WIN32) ADD_NEKTAR_TEST (Gmsh/Scalar_Windows) ELSE () ADD_NEKTAR_TEST (Gmsh/Scalar) ENDIF () ADD_NEKTAR_TEST (Gmsh/SquareQuad) +ADD_NEKTAR_TEST (Gmsh/SquareQuadLinear) ADD_NEKTAR_TEST (Gmsh/SquareTri) +ADD_NEKTAR_TEST (Gmsh/SquareTriLinear) # Nektar tests ADD_NEKTAR_TEST (Nektar/HexLinear) ADD_NEKTAR_TEST (Nektar/Tube45) diff --git a/utilities/NekMesh/Tests/Gmsh/CubeHexLinear.msh b/utilities/NekMesh/Tests/Gmsh/CubeHexLinear.msh new file mode 100644 index 000000000..c8a119e38 --- /dev/null +++ b/utilities/NekMesh/Tests/Gmsh/CubeHexLinear.msh @@ -0,0 +1,68 @@ +$MeshFormat +2.2 0 8 +$EndMeshFormat +$Nodes +27 +1 0 0 0 +2 1 0 0 +3 0 1 0 +4 1 1 0 +5 0 0 1 +6 1 0 1 +7 1 1 1 +8 0 1 1 +9 0.5 0 0 +10 0.5 1 0 +11 0 0.5 0 +12 1 0.5 0 +13 0.5 0 1 +14 1 0.5 1 +15 0.5 1 1 +16 0 0.5 1 +17 0 0 0.5 +18 1 0 0.5 +19 1 1 0.5 +20 0 1 0.5 +21 0.5 0.5 0 +22 0.5 0 0.5 +23 1 0.5 0.5 +24 0.5 1 0.5 +25 0 0.5 0.5 +26 0.5 0.5 1 +27 0.5 0.5 0.5 +$EndNodes +$Elements +32 +1 3 2 1 5 1 9 21 11 +2 3 2 1 5 11 21 10 3 +3 3 2 1 5 9 2 12 21 +4 3 2 1 5 21 12 4 10 +5 3 2 2 14 1 9 22 17 +6 3 2 2 14 17 22 13 5 +7 3 2 2 14 9 2 18 22 +8 3 2 2 14 22 18 6 13 +9 3 2 3 18 2 12 23 18 +10 3 2 3 18 18 23 14 6 +11 3 2 3 18 12 4 19 23 +12 3 2 3 18 23 19 7 14 +13 3 2 4 22 3 20 24 10 +14 3 2 4 22 20 8 15 24 +15 3 2 4 22 10 24 19 4 +16 3 2 4 22 24 15 7 19 +17 3 2 5 26 1 17 25 11 +18 3 2 5 26 17 5 16 25 +19 3 2 5 26 11 25 20 3 +20 3 2 5 26 25 16 8 20 +21 3 2 6 27 5 13 26 16 +22 3 2 6 27 16 26 15 8 +23 3 2 6 27 13 6 14 26 +24 3 2 6 27 26 14 7 15 +25 5 2 0 1 1 9 21 11 17 22 27 25 +26 5 2 0 1 17 22 27 25 5 13 26 16 +27 5 2 0 1 11 21 10 3 25 27 24 20 +28 5 2 0 1 25 27 24 20 16 26 15 8 +29 5 2 0 1 9 2 12 21 22 18 23 27 +30 5 2 0 1 22 18 23 27 13 6 14 26 +31 5 2 0 1 21 12 4 10 27 23 19 24 +32 5 2 0 1 27 23 19 24 26 14 7 15 +$EndElements diff --git a/utilities/NekMesh/Tests/Gmsh/CubeHexLinear.tst b/utilities/NekMesh/Tests/Gmsh/CubeHexLinear.tst new file mode 100644 index 000000000..564d1dc6b --- /dev/null +++ b/utilities/NekMesh/Tests/Gmsh/CubeHexLinear.tst @@ -0,0 +1,19 @@ + + + Gmsh linear hex with order 7 output + NekMesh + -m jac:list CubeHexLinear.msh CubeHex.xml:xml:test:order=7 + + CubeHexLinear.msh + + + + ^Total negative Jacobians: (\d+) + + + 0 + + + + + diff --git a/utilities/NekMesh/Tests/Gmsh/CubePrismLinear.msh b/utilities/NekMesh/Tests/Gmsh/CubePrismLinear.msh new file mode 100644 index 000000000..89f0e150c --- /dev/null +++ b/utilities/NekMesh/Tests/Gmsh/CubePrismLinear.msh @@ -0,0 +1,84 @@ +$MeshFormat +2.2 0 8 +$EndMeshFormat +$Nodes +27 +1 0 0 0 +2 1 0 0 +3 0 1 0 +4 1 1 0 +5 0 0 1 +6 1 0 1 +7 1 1 1 +8 0 1 1 +9 0.5 0 0 +10 0.5 1 0 +11 0 0.5 0 +12 1 0.5 0 +13 0.5 0 1 +14 1 0.5 1 +15 0.5 1 1 +16 0 0.5 1 +17 0 0 0.5 +18 1 0 0.5 +19 1 1 0.5 +20 0 1 0.5 +21 0.5 0.5 0 +22 0.5 0 0.5 +23 1 0.5 0.5 +24 0.5 1 0.5 +25 0 0.5 0.5 +26 0.5 0.5 1 +27 0.5 0.5 0.5 +$EndNodes +$Elements +48 +1 2 2 1 5 1 9 21 +2 2 2 1 5 1 21 11 +3 2 2 1 5 11 21 10 +4 2 2 1 5 11 10 3 +5 2 2 1 5 9 2 12 +6 2 2 1 5 9 12 21 +7 2 2 1 5 21 12 4 +8 2 2 1 5 21 4 10 +9 2 2 6 27 5 13 26 +10 2 2 6 27 5 26 16 +11 2 2 6 27 16 26 15 +12 2 2 6 27 16 15 8 +13 2 2 6 27 13 6 14 +14 2 2 6 27 13 14 26 +15 2 2 6 27 26 14 7 +16 2 2 6 27 26 7 15 +17 3 2 2 14 1 9 22 17 +18 3 2 2 14 17 22 13 5 +19 3 2 2 14 9 2 18 22 +20 3 2 2 14 22 18 6 13 +21 3 2 3 18 2 12 23 18 +22 3 2 3 18 18 23 14 6 +23 3 2 3 18 12 4 19 23 +24 3 2 3 18 23 19 7 14 +25 3 2 4 22 3 20 24 10 +26 3 2 4 22 20 8 15 24 +27 3 2 4 22 10 24 19 4 +28 3 2 4 22 24 15 7 19 +29 3 2 5 26 1 17 25 11 +30 3 2 5 26 17 5 16 25 +31 3 2 5 26 11 25 20 3 +32 3 2 5 26 25 16 8 20 +33 6 2 0 1 1 9 21 17 22 27 +34 6 2 0 1 17 22 27 5 13 26 +35 6 2 0 1 1 21 11 17 27 25 +36 6 2 0 1 17 27 25 5 26 16 +37 6 2 0 1 11 21 10 25 27 24 +38 6 2 0 1 25 27 24 16 26 15 +39 6 2 0 1 11 10 3 25 24 20 +40 6 2 0 1 25 24 20 16 15 8 +41 6 2 0 1 9 2 12 22 18 23 +42 6 2 0 1 22 18 23 13 6 14 +43 6 2 0 1 9 12 21 22 23 27 +44 6 2 0 1 22 23 27 13 14 26 +45 6 2 0 1 21 12 4 27 23 19 +46 6 2 0 1 27 23 19 26 14 7 +47 6 2 0 1 21 4 10 27 19 24 +48 6 2 0 1 27 19 24 26 7 15 +$EndElements diff --git a/utilities/NekMesh/Tests/Gmsh/CubePrismLinear.tst b/utilities/NekMesh/Tests/Gmsh/CubePrismLinear.tst new file mode 100644 index 000000000..a94c71802 --- /dev/null +++ b/utilities/NekMesh/Tests/Gmsh/CubePrismLinear.tst @@ -0,0 +1,19 @@ + + + Gmsh linear prism with order 4 output + NekMesh + -m jac:list CubePrism.msh CubePrism.xml:xml:test:order=4 + + CubePrism.msh + + + + ^Total negative Jacobians: (\d+) + + + 0 + + + + + diff --git a/utilities/NekMesh/Tests/Gmsh/CubeTetLinear.msh b/utilities/NekMesh/Tests/Gmsh/CubeTetLinear.msh new file mode 100644 index 000000000..94c6ab6a7 --- /dev/null +++ b/utilities/NekMesh/Tests/Gmsh/CubeTetLinear.msh @@ -0,0 +1,343 @@ +$MeshFormat +2.2 0 8 +$EndMeshFormat +$Nodes +64 +1 0 0 0 +2 1 0 0 +3 0 1 0 +4 1 1 0 +5 0 0 1 +6 1 0 1 +7 1 1 1 +8 0 1 1 +9 0.3333333333333333 0 0 +10 0.6666666666666666 0 0 +11 0.3333333333333333 1 0 +12 0.6666666666666666 1 0 +13 0 0.3333333333333333 0 +14 0 0.6666666666666666 0 +15 1 0.3333333333333333 0 +16 1 0.6666666666666666 0 +17 0.3333333333333333 0 1 +18 0.6666666666666666 0 1 +19 1 0.3333333333333333 1 +20 1 0.6666666666666666 1 +21 0.6666666666666666 1 1 +22 0.3333333333333333 1 1 +23 0 0.6666666666666666 1 +24 0 0.3333333333333333 1 +25 0 0 0.3333333333333333 +26 0 0 0.6666666666666666 +27 1 0 0.3333333333333333 +28 1 0 0.6666666666666666 +29 1 1 0.3333333333333333 +30 1 1 0.6666666666666666 +31 0 1 0.3333333333333333 +32 0 1 0.6666666666666666 +33 0.3333333333333333 0.3333333333333333 0 +34 0.3333333333333333 0.6666666666666666 0 +35 0.6666666666666666 0.3333333333333333 0 +36 0.6666666666666666 0.6666666666666666 0 +37 0.3333333333333333 0 0.3333333333333333 +38 0.3333333333333333 0 0.6666666666666666 +39 0.6666666666666666 0 0.3333333333333333 +40 0.6666666666666666 0 0.6666666666666666 +41 1 0.3333333333333333 0.3333333333333333 +42 1 0.3333333333333333 0.6666666666666666 +43 1 0.6666666666666666 0.3333333333333333 +44 1 0.6666666666666666 0.6666666666666666 +45 0.3333333333333333 1 0.3333333333333333 +46 0.3333333333333333 1 0.6666666666666666 +47 0.6666666666666666 1 0.3333333333333333 +48 0.6666666666666666 1 0.6666666666666666 +49 0 0.3333333333333333 0.3333333333333333 +50 0 0.3333333333333333 0.6666666666666666 +51 0 0.6666666666666666 0.3333333333333333 +52 0 0.6666666666666666 0.6666666666666666 +53 0.3333333333333333 0.3333333333333333 1 +54 0.3333333333333333 0.6666666666666666 1 +55 0.6666666666666666 0.3333333333333333 1 +56 0.6666666666666666 0.6666666666666666 1 +57 0.3333333333333333 0.3333333333333333 0.3333333333333333 +58 0.3333333333333333 0.3333333333333333 0.6666666666666666 +59 0.3333333333333333 0.6666666666666666 0.3333333333333333 +60 0.3333333333333333 0.6666666666666666 0.6666666666666666 +61 0.6666666666666666 0.3333333333333333 0.3333333333333333 +62 0.6666666666666666 0.3333333333333333 0.6666666666666666 +63 0.6666666666666666 0.6666666666666666 0.3333333333333333 +64 0.6666666666666666 0.6666666666666666 0.6666666666666666 +$EndNodes +$Elements +270 +1 2 2 1 5 13 1 33 +2 2 2 1 5 1 9 33 +3 2 2 1 5 14 13 34 +4 2 2 1 5 13 33 34 +5 2 2 1 5 3 14 11 +6 2 2 1 5 14 34 11 +7 2 2 1 5 33 9 35 +8 2 2 1 5 9 10 35 +9 2 2 1 5 34 33 36 +10 2 2 1 5 33 35 36 +11 2 2 1 5 11 34 12 +12 2 2 1 5 34 36 12 +13 2 2 1 5 35 10 15 +14 2 2 1 5 10 2 15 +15 2 2 1 5 36 35 16 +16 2 2 1 5 35 15 16 +17 2 2 1 5 12 36 4 +18 2 2 1 5 36 16 4 +19 2 2 2 14 25 1 37 +20 2 2 2 14 1 9 37 +21 2 2 2 14 26 25 38 +22 2 2 2 14 25 37 38 +23 2 2 2 14 5 26 17 +24 2 2 2 14 26 38 17 +25 2 2 2 14 37 9 39 +26 2 2 2 14 9 10 39 +27 2 2 2 14 38 37 40 +28 2 2 2 14 37 39 40 +29 2 2 2 14 17 38 18 +30 2 2 2 14 38 40 18 +31 2 2 2 14 39 10 2 +32 2 2 2 14 39 2 27 +33 2 2 2 14 40 39 27 +34 2 2 2 14 40 27 28 +35 2 2 2 14 18 40 28 +36 2 2 2 14 18 28 6 +37 2 2 3 18 27 2 41 +38 2 2 3 18 2 15 41 +39 2 2 3 18 28 27 42 +40 2 2 3 18 27 41 42 +41 2 2 3 18 6 28 19 +42 2 2 3 18 28 42 19 +43 2 2 3 18 41 15 43 +44 2 2 3 18 15 16 43 +45 2 2 3 18 42 41 44 +46 2 2 3 18 41 43 44 +47 2 2 3 18 19 42 20 +48 2 2 3 18 42 44 20 +49 2 2 3 18 43 16 4 +50 2 2 3 18 43 4 29 +51 2 2 3 18 44 43 29 +52 2 2 3 18 44 29 30 +53 2 2 3 18 20 44 30 +54 2 2 3 18 20 30 7 +55 2 2 4 22 31 11 3 +56 2 2 4 22 31 45 11 +57 2 2 4 22 32 46 31 +58 2 2 4 22 31 46 45 +59 2 2 4 22 8 22 32 +60 2 2 4 22 32 22 46 +61 2 2 4 22 45 47 11 +62 2 2 4 22 11 47 12 +63 2 2 4 22 46 48 45 +64 2 2 4 22 45 48 47 +65 2 2 4 22 22 48 46 +66 2 2 4 22 22 21 48 +67 2 2 4 22 47 4 12 +68 2 2 4 22 47 29 4 +69 2 2 4 22 48 29 47 +70 2 2 4 22 48 30 29 +71 2 2 4 22 21 30 48 +72 2 2 4 22 21 7 30 +73 2 2 5 26 25 49 1 +74 2 2 5 26 1 49 13 +75 2 2 5 26 26 50 25 +76 2 2 5 26 25 50 49 +77 2 2 5 26 5 24 26 +78 2 2 5 26 26 24 50 +79 2 2 5 26 49 51 13 +80 2 2 5 26 13 51 14 +81 2 2 5 26 50 52 49 +82 2 2 5 26 49 52 51 +83 2 2 5 26 24 23 50 +84 2 2 5 26 50 23 52 +85 2 2 5 26 51 3 14 +86 2 2 5 26 51 31 3 +87 2 2 5 26 52 31 51 +88 2 2 5 26 52 32 31 +89 2 2 5 26 23 32 52 +90 2 2 5 26 23 8 32 +91 2 2 6 27 5 17 53 +92 2 2 6 27 5 53 24 +93 2 2 6 27 24 53 54 +94 2 2 6 27 24 54 23 +95 2 2 6 27 23 54 22 +96 2 2 6 27 23 22 8 +97 2 2 6 27 17 18 55 +98 2 2 6 27 17 55 53 +99 2 2 6 27 53 55 56 +100 2 2 6 27 53 56 54 +101 2 2 6 27 54 56 21 +102 2 2 6 27 54 21 22 +103 2 2 6 27 18 6 19 +104 2 2 6 27 18 19 55 +105 2 2 6 27 55 19 20 +106 2 2 6 27 55 20 56 +107 2 2 6 27 56 20 7 +108 2 2 6 27 56 7 21 +109 4 2 0 1 1 9 33 37 +110 4 2 0 1 37 25 57 33 +111 4 2 0 1 1 25 37 33 +112 4 2 0 1 25 37 57 38 +113 4 2 0 1 38 26 58 57 +114 4 2 0 1 25 26 38 57 +115 4 2 0 1 26 38 58 17 +116 4 2 0 1 17 5 53 58 +117 4 2 0 1 26 5 17 58 +118 4 2 0 1 1 33 13 49 +119 4 2 0 1 57 25 49 33 +120 4 2 0 1 33 25 49 1 +121 4 2 0 1 25 57 49 50 +122 4 2 0 1 58 26 50 57 +123 4 2 0 1 57 26 50 25 +124 4 2 0 1 26 58 50 24 +125 4 2 0 1 53 5 24 58 +126 4 2 0 1 58 5 24 26 +127 4 2 0 1 13 33 34 49 +128 4 2 0 1 57 49 59 33 +129 4 2 0 1 33 49 59 34 +130 4 2 0 1 49 57 59 50 +131 4 2 0 1 58 50 60 57 +132 4 2 0 1 57 50 60 59 +133 4 2 0 1 50 58 60 24 +134 4 2 0 1 53 24 54 60 +135 4 2 0 1 58 24 53 60 +136 4 2 0 1 13 34 14 51 +137 4 2 0 1 59 49 51 34 +138 4 2 0 1 34 49 51 13 +139 4 2 0 1 49 59 51 52 +140 4 2 0 1 60 50 52 59 +141 4 2 0 1 59 50 52 49 +142 4 2 0 1 50 60 52 23 +143 4 2 0 1 54 24 23 60 +144 4 2 0 1 60 24 23 50 +145 4 2 0 1 14 34 11 51 +146 4 2 0 1 59 51 45 34 +147 4 2 0 1 34 51 45 11 +148 4 2 0 1 51 59 45 52 +149 4 2 0 1 60 52 46 59 +150 4 2 0 1 59 52 46 45 +151 4 2 0 1 52 60 46 23 +152 4 2 0 1 54 23 22 60 +153 4 2 0 1 60 23 22 46 +154 4 2 0 1 14 11 3 51 +155 4 2 0 1 45 51 31 11 +156 4 2 0 1 11 51 31 3 +157 4 2 0 1 51 45 31 52 +158 4 2 0 1 46 52 32 31 +159 4 2 0 1 45 52 46 31 +160 4 2 0 1 52 46 32 23 +161 4 2 0 1 22 23 8 32 +162 4 2 0 1 46 23 22 32 +163 4 2 0 1 9 10 35 39 +164 4 2 0 1 39 37 61 35 +165 4 2 0 1 9 37 39 35 +166 4 2 0 1 37 39 61 40 +167 4 2 0 1 40 38 62 61 +168 4 2 0 1 37 38 40 61 +169 4 2 0 1 38 40 62 18 +170 4 2 0 1 18 17 55 62 +171 4 2 0 1 38 17 18 62 +172 4 2 0 1 9 35 33 37 +173 4 2 0 1 61 37 57 33 +174 4 2 0 1 35 37 61 33 +175 4 2 0 1 37 61 57 38 +176 4 2 0 1 62 38 58 57 +177 4 2 0 1 61 38 62 57 +178 4 2 0 1 38 62 58 17 +179 4 2 0 1 55 17 53 58 +180 4 2 0 1 62 17 55 58 +181 4 2 0 1 33 35 36 63 +182 4 2 0 1 61 57 63 33 +183 4 2 0 1 35 33 61 63 +184 4 2 0 1 57 61 63 64 +185 4 2 0 1 62 58 64 57 +186 4 2 0 1 61 57 62 64 +187 4 2 0 1 58 62 64 56 +188 4 2 0 1 55 53 56 58 +189 4 2 0 1 62 58 55 56 +190 4 2 0 1 33 36 34 63 +191 4 2 0 1 63 57 59 33 +192 4 2 0 1 34 33 63 59 +193 4 2 0 1 57 63 59 64 +194 4 2 0 1 64 58 60 57 +195 4 2 0 1 59 57 64 60 +196 4 2 0 1 58 64 60 56 +197 4 2 0 1 56 53 54 60 +198 4 2 0 1 58 53 56 60 +199 4 2 0 1 34 36 12 47 +200 4 2 0 1 63 59 47 34 +201 4 2 0 1 36 34 63 47 +202 4 2 0 1 59 63 47 48 +203 4 2 0 1 64 60 48 59 +204 4 2 0 1 63 59 64 48 +205 4 2 0 1 60 64 48 21 +206 4 2 0 1 56 54 21 60 +207 4 2 0 1 64 60 56 21 +208 4 2 0 1 34 12 11 47 +209 4 2 0 1 47 59 45 34 +210 4 2 0 1 11 34 47 45 +211 4 2 0 1 59 47 45 48 +212 4 2 0 1 48 60 46 59 +213 4 2 0 1 45 59 48 46 +214 4 2 0 1 60 48 46 22 +215 4 2 0 1 21 54 22 60 +216 4 2 0 1 48 60 21 22 +217 4 2 0 1 10 2 15 41 +218 4 2 0 1 27 39 41 2 +219 4 2 0 1 2 39 41 10 +220 4 2 0 1 39 27 41 42 +221 4 2 0 1 28 40 42 27 +222 4 2 0 1 27 40 42 39 +223 4 2 0 1 40 28 42 19 +224 4 2 0 1 6 18 19 28 +225 4 2 0 1 28 18 19 40 +226 4 2 0 1 10 15 35 41 +227 4 2 0 1 41 39 61 35 +228 4 2 0 1 10 39 41 35 +229 4 2 0 1 39 41 61 42 +230 4 2 0 1 42 40 62 61 +231 4 2 0 1 39 40 42 61 +232 4 2 0 1 40 42 62 19 +233 4 2 0 1 19 18 55 62 +234 4 2 0 1 40 18 19 62 +235 4 2 0 1 35 15 16 43 +236 4 2 0 1 41 61 43 35 +237 4 2 0 1 15 35 41 43 +238 4 2 0 1 61 41 43 44 +239 4 2 0 1 42 62 44 61 +240 4 2 0 1 41 61 42 44 +241 4 2 0 1 62 42 44 20 +242 4 2 0 1 19 55 20 62 +243 4 2 0 1 42 62 19 20 +244 4 2 0 1 35 16 36 43 +245 4 2 0 1 43 61 63 35 +246 4 2 0 1 36 35 43 63 +247 4 2 0 1 61 43 63 44 +248 4 2 0 1 44 62 64 61 +249 4 2 0 1 63 61 44 64 +250 4 2 0 1 62 44 64 20 +251 4 2 0 1 20 55 56 62 +252 4 2 0 1 64 62 20 56 +253 4 2 0 1 36 16 4 43 +254 4 2 0 1 43 63 29 36 +255 4 2 0 1 4 36 43 29 +256 4 2 0 1 63 43 29 44 +257 4 2 0 1 44 64 30 29 +258 4 2 0 1 63 64 44 29 +259 4 2 0 1 64 44 30 20 +260 4 2 0 1 20 56 7 30 +261 4 2 0 1 64 56 20 30 +262 4 2 0 1 36 4 12 47 +263 4 2 0 1 29 63 47 36 +264 4 2 0 1 4 36 29 47 +265 4 2 0 1 63 29 47 48 +266 4 2 0 1 30 64 48 29 +267 4 2 0 1 29 64 48 63 +268 4 2 0 1 64 30 48 21 +269 4 2 0 1 7 56 21 30 +270 4 2 0 1 30 56 21 64 +$EndElements diff --git a/utilities/NekMesh/Tests/Gmsh/CubeTetLinear.tst b/utilities/NekMesh/Tests/Gmsh/CubeTetLinear.tst new file mode 100644 index 000000000..85e201fb0 --- /dev/null +++ b/utilities/NekMesh/Tests/Gmsh/CubeTetLinear.tst @@ -0,0 +1,19 @@ + + + Gmsh tet cube, convert to order 9 + NekMesh + -m jac:list CubeTetLinear.msh CubeTet.xml:xml:order=9:test + + CubeTetLinear.msh + + + + ^Total negative Jacobians: (\d+) + + + 0 + + + + + diff --git a/utilities/NekMesh/Tests/Gmsh/SquareQuadLinear.msh b/utilities/NekMesh/Tests/Gmsh/SquareQuadLinear.msh new file mode 100644 index 000000000..6bfcf5868 --- /dev/null +++ b/utilities/NekMesh/Tests/Gmsh/SquareQuadLinear.msh @@ -0,0 +1,66 @@ +$MeshFormat +2.2 0 8 +$EndMeshFormat +$Nodes +25 +1 -0.5 -0.5 0 +2 0.5 -0.5 0 +3 -0.5 0.5 0 +4 0.5 0.5 0 +5 -0.25 -0.5 0 +6 0 -0.5 0 +7 0.25 -0.5 0 +8 -0.25 0.5 0 +9 0 0.5 0 +10 0.25 0.5 0 +11 -0.5 -0.25 0 +12 -0.5 0 0 +13 -0.5 0.25 0 +14 0.5 -0.25 0 +15 0.5 0 0 +16 0.5 0.25 0 +17 -0.25 -0.25 0 +18 -0.25 0 0 +19 -0.25 0.25 0 +20 0 -0.25 0 +21 0 0 0 +22 0 0.25 0 +23 0.25 -0.25 0 +24 0.25 0 0 +25 0.25 0.25 0 +$EndNodes +$Elements +32 +1 1 2 1 1 1 5 +2 1 2 1 1 5 6 +3 1 2 1 1 6 7 +4 1 2 1 1 7 2 +5 1 2 3 2 3 8 +6 1 2 3 2 8 9 +7 1 2 3 2 9 10 +8 1 2 3 2 10 4 +9 1 2 4 3 1 11 +10 1 2 4 3 11 12 +11 1 2 4 3 12 13 +12 1 2 4 3 13 3 +13 1 2 2 4 2 14 +14 1 2 2 4 14 15 +15 1 2 2 4 15 16 +16 1 2 2 4 16 4 +17 3 2 0 5 1 5 17 11 +18 3 2 0 5 11 17 18 12 +19 3 2 0 5 12 18 19 13 +20 3 2 0 5 13 19 8 3 +21 3 2 0 5 5 6 20 17 +22 3 2 0 5 17 20 21 18 +23 3 2 0 5 18 21 22 19 +24 3 2 0 5 19 22 9 8 +25 3 2 0 5 6 7 23 20 +26 3 2 0 5 20 23 24 21 +27 3 2 0 5 21 24 25 22 +28 3 2 0 5 22 25 10 9 +29 3 2 0 5 7 2 14 23 +30 3 2 0 5 23 14 15 24 +31 3 2 0 5 24 15 16 25 +32 3 2 0 5 25 16 4 10 +$EndElements diff --git a/utilities/NekMesh/Tests/Gmsh/SquareQuadLinear.tst b/utilities/NekMesh/Tests/Gmsh/SquareQuadLinear.tst new file mode 100644 index 000000000..94bfb7d02 --- /dev/null +++ b/utilities/NekMesh/Tests/Gmsh/SquareQuadLinear.tst @@ -0,0 +1,19 @@ + + + Gmsh linear mesh with order 7 output + NekMesh + -m jac:list SquareQuadLinear.msh SquareQuad.xml:xml:test:order=7 + + SquareQuadLinear.msh + + + + ^Total negative Jacobians: (\d+) + + + 0 + + + + + diff --git a/utilities/NekMesh/Tests/Gmsh/SquareTriLinear.msh b/utilities/NekMesh/Tests/Gmsh/SquareTriLinear.msh new file mode 100644 index 000000000..f713acfa7 --- /dev/null +++ b/utilities/NekMesh/Tests/Gmsh/SquareTriLinear.msh @@ -0,0 +1,82 @@ +$MeshFormat +2.2 0 8 +$EndMeshFormat +$Nodes +25 +1 -1 0 0 +2 1 0 0 +3 -1 2 0 +4 1 2 0 +5 -0.5 0 0 +6 0 0 0 +7 0.5 0 0 +8 -0.5 2 0 +9 0 2 0 +10 0.5 2 0 +11 -1 0.5 0 +12 -1 1 0 +13 -1 1.5 0 +14 1 0.5 0 +15 1 1 0 +16 1 1.5 0 +17 -0.5 0.5 0 +18 -0.5 1 0 +19 -0.5 1.5 0 +20 0 0.5 0 +21 0 1 0 +22 0 1.5 0 +23 0.5 0.5 0 +24 0.5 1 0 +25 0.5 1.5 0 +$EndNodes +$Elements +48 +1 1 2 1 1 1 5 +2 1 2 1 1 5 6 +3 1 2 1 1 6 7 +4 1 2 1 1 7 2 +5 1 2 3 2 3 8 +6 1 2 3 2 8 9 +7 1 2 3 2 9 10 +8 1 2 3 2 10 4 +9 1 2 4 3 1 11 +10 1 2 4 3 11 12 +11 1 2 4 3 12 13 +12 1 2 4 3 13 3 +13 1 2 2 4 2 14 +14 1 2 2 4 14 15 +15 1 2 2 4 15 16 +16 1 2 2 4 16 4 +17 2 2 0 5 1 5 17 +18 2 2 0 5 1 17 11 +19 2 2 0 5 11 17 18 +20 2 2 0 5 11 18 12 +21 2 2 0 5 12 18 19 +22 2 2 0 5 12 19 13 +23 2 2 0 5 13 19 8 +24 2 2 0 5 13 8 3 +25 2 2 0 5 5 6 20 +26 2 2 0 5 5 20 17 +27 2 2 0 5 17 20 21 +28 2 2 0 5 17 21 18 +29 2 2 0 5 18 21 22 +30 2 2 0 5 18 22 19 +31 2 2 0 5 19 22 9 +32 2 2 0 5 19 9 8 +33 2 2 0 5 6 7 23 +34 2 2 0 5 6 23 20 +35 2 2 0 5 20 23 24 +36 2 2 0 5 20 24 21 +37 2 2 0 5 21 24 25 +38 2 2 0 5 21 25 22 +39 2 2 0 5 22 25 10 +40 2 2 0 5 22 10 9 +41 2 2 0 5 7 2 14 +42 2 2 0 5 7 14 23 +43 2 2 0 5 23 14 15 +44 2 2 0 5 23 15 24 +45 2 2 0 5 24 15 16 +46 2 2 0 5 24 16 25 +47 2 2 0 5 25 16 4 +48 2 2 0 5 25 4 10 +$EndElements diff --git a/utilities/NekMesh/Tests/Gmsh/SquareTriLinear.tst b/utilities/NekMesh/Tests/Gmsh/SquareTriLinear.tst new file mode 100644 index 000000000..4292e8577 --- /dev/null +++ b/utilities/NekMesh/Tests/Gmsh/SquareTriLinear.tst @@ -0,0 +1,19 @@ + + + Gmsh linear tri square, order 9 output + NekMesh + -m jac:list SquareTriLinear.msh SquareTri.xml:xml:test:order=9 + + SquareTriLinear.msh + + + + ^Total negative Jacobians: (\d+) + + + 0 + + + + + -- GitLab From db7642d8b8b2467ba87c222ae99a6204384ef839 Mon Sep 17 00:00:00 2001 From: David Moxey Date: Tue, 16 Aug 2016 10:19:46 +0100 Subject: [PATCH 21/25] Add user guide documentation, fix Windows compilation error. --- docs/user-guide/utilities/nekmesh.tex | 32 ++++++++++++++++-------- library/NekMeshUtils/MeshElements/Mesh.h | 4 +-- 2 files changed, 23 insertions(+), 13 deletions(-) diff --git a/docs/user-guide/utilities/nekmesh.tex b/docs/user-guide/utilities/nekmesh.tex index 6b4c927a3..3b0a9e4c8 100644 --- a/docs/user-guide/utilities/nekmesh.tex +++ b/docs/user-guide/utilities/nekmesh.tex @@ -26,13 +26,13 @@ through a basic example to show how a mesh can be converted from the widely-used mesh-generator \gmsh to the XML file format. \begin{notebox} - The default since Jan 2016 is to output the .xml files in a + The default since January 2016 is to output the \inltt{.xml} files in a compressed form where the VERTEX, EDGES, FACES, ELEMENTS and CURVED - information is compressed into binary format which is then - converted into base64. This is identified for each section by the attribute - \inltt{COMPRESSED="B64Z-LittleEndian''}. To output - in ascii format add the module option ``:xml:uncompress'' to the .xml file, - i.e. \\ \inltt{ \mc file.msh newfile.xml:xml:uncompress} + information is compressed into binary format which is then converted into + base64. This is identified for each section by the attribute + \inltt{COMPRESSED="B64Z-LittleEndian''}. To output in ascii format add the + module option ``:xml:uncompress'' to the \inltt{.xml} file, i.e. \\ \inltt{ + \mc file.msh newfile.xml:xml:uncompress} \end{notebox} \section{Exporting a mesh from \gmsh} @@ -274,7 +274,9 @@ The following output formats are supported: \toprule \textbf{Format} & \textbf{Extension} & \textbf{High-order} & \textbf{Notes}\\ \midrule - Gmsh & \texttt{msh} & \cmark & Curvature output is highly experimental.\\ + Gmsh & \texttt{msh} & \cmark & High-order hexes, quads, tetrahedra and + triangles are supported up to arbitrary order. Prisms supported up to order + 4, pyramids up to order 1.\\ Nektar++ & \texttt{xml} & \cmark & Most functionality supported. \\ VTK & \texttt{vtk} & \xmark & Experimental. Only ASCII triangular data is supported. \\ \bottomrule @@ -288,13 +290,21 @@ meshes since robustness is not guaranteed. The default for \texttt{xml} is into binary data which has been converted into base64. If you wish to see an ascii output you need to specify the output module option \inltt{uncompress} by executing: - +% \begin{lstlisting}[style=BashInputStyle] NekMesh Mesh.msh output.xml:xml:uncompress \end{lstlisting} - -In the rest of these subsections, we discuss the various processing modules -available within \mc. +% +Finally, both the Gmsh and Nektar++ output modules support an \inltt{order} +parameter, which allows you to generate a mesh of a uniform polynomial +order. This is used in the same manner as the above, so that the command +% +\begin{lstlisting}[style=BashInputStyle] +NekMesh Mesh.msh output.msh:msh:order=7 +\end{lstlisting} +% +will generate an order 7 Gmsh mesh. In the rest of these subsections, we discuss +the various processing modules available within \mc. \subsection{Extract surfaces from a mesh} diff --git a/library/NekMeshUtils/MeshElements/Mesh.h b/library/NekMeshUtils/MeshElements/Mesh.h index 6cd633ce6..19d826b13 100644 --- a/library/NekMeshUtils/MeshElements/Mesh.h +++ b/library/NekMeshUtils/MeshElements/Mesh.h @@ -133,8 +133,8 @@ public: /// Returns the total number of entities in the mesh. NEKMESHUTILS_EXPORT unsigned int GetNumEntities(); - void MakeOrder(int order, - LibUtilities::PointsType distType); + NEKMESHUTILS_EXPORT void MakeOrder(int order, + LibUtilities::PointsType distType); }; /// Shared pointer to a mesh. typedef boost::shared_ptr MeshSharedPtr; -- GitLab From 76c3d3648af1db1f414bc887cd16e6fb19d2ddf6 Mon Sep 17 00:00:00 2001 From: David Moxey Date: Tue, 16 Aug 2016 10:24:02 +0100 Subject: [PATCH 22/25] Add a little doxygen --- library/NekMeshUtils/MeshElements/Mesh.cpp | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/library/NekMeshUtils/MeshElements/Mesh.cpp b/library/NekMeshUtils/MeshElements/Mesh.cpp index dae08817c..900052306 100644 --- a/library/NekMeshUtils/MeshElements/Mesh.cpp +++ b/library/NekMeshUtils/MeshElements/Mesh.cpp @@ -82,7 +82,18 @@ unsigned int Mesh::GetNumEntities() } /** - * @brief Convert this mesh into a high-order mesh. + * @brief Convert this mesh into a mesh of uniform polynomial order @p order + * with a curve point distribution @p distType. + * + * This routine adds curvature points into a mesh so that the resulting elements + * are all of a uniform order @p order and all high-order vertices are + * consistently ordered. It proceeds in a bottom-up fashion: + * + * - First construct all edge, face and elemental geometry mappings. + * - Then call the local MakeOrder functions on each edge, face and element of + * dimension Mesh::m_expDim. + * - Finally, any boundary elements are updated so that they have the same + * interior degrees of freedom as their corresponding edge or face links. */ void Mesh::MakeOrder(int order, LibUtilities::PointsType distType) @@ -119,6 +130,10 @@ void Mesh::MakeOrder(int order, pTypes[LibUtilities::eTetrahedron] = LibUtilities::eNodalTetElec; pTypes[LibUtilities::eHexahedron] = LibUtilities::eGaussLobattoLegendre; } + else + { + ASSERTL1(false, "Mesh::MakeOrder does not support this points type."); + } // Begin by generating Nektar++ geometry objects for edges, faces and // elements so that we don't affect any neighbouring elements in the mesh as -- GitLab From 57a87b1d2db0238ddd2c93833b335cd8c04938b5 Mon Sep 17 00:00:00 2001 From: David Moxey Date: Tue, 16 Aug 2016 10:27:04 +0100 Subject: [PATCH 23/25] Update CHANGELOG --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index a93590d4f..77315a775 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -43,6 +43,7 @@ v4.4.0 - New module for inserting an alternate high-order surface into the working mesh (!669) - Improvements to mesh linearisation module (!659) +- Add support for Gmsh high-order output (!679) **FieldConvert:** - Move all modules to a new library, FieldUtils, to support post-processing -- GitLab From 3790ce249abb253702ac586431fba0ef17517fbe Mon Sep 17 00:00:00 2001 From: David Moxey Date: Tue, 16 Aug 2016 20:38:55 +0100 Subject: [PATCH 24/25] Fix another Windows compile problem, add some more comments --- library/NekMeshUtils/MeshElements/Edge.h | 1 + library/NekMeshUtils/MeshElements/Element.cpp | 12 ++++--- library/NekMeshUtils/MeshElements/Element.h | 31 +++++++++++++++++-- library/NekMeshUtils/MeshElements/Face.h | 4 ++- .../NekMeshUtils/MeshElements/Quadrilateral.h | 2 +- 5 files changed, 41 insertions(+), 9 deletions(-) diff --git a/library/NekMeshUtils/MeshElements/Edge.h b/library/NekMeshUtils/MeshElements/Edge.h index 2c2d6c519..bd5b5c941 100644 --- a/library/NekMeshUtils/MeshElements/Edge.h +++ b/library/NekMeshUtils/MeshElements/Edge.h @@ -170,6 +170,7 @@ public: return ret; } + /// Make this edge an order @p order edge. @see Element::MakeOrder. void MakeOrder(int order, SpatialDomains::GeometrySharedPtr geom, LibUtilities::PointsType edgeType, diff --git a/library/NekMeshUtils/MeshElements/Element.cpp b/library/NekMeshUtils/MeshElements/Element.cpp index 9cd7bb5a3..f7402f634 100644 --- a/library/NekMeshUtils/MeshElements/Element.cpp +++ b/library/NekMeshUtils/MeshElements/Element.cpp @@ -71,8 +71,10 @@ Element::Element(ElmtConfig pConf, unsigned int pNumNodes, * searched and the corresponding edge/face nodes are updated to * maintain consistency. * - * @param p Index of the vertex to replace. - * @param pNew New vertex. + * @param p Index of the vertex to replace. + * @param pNew New vertex. + * @param descend If true, we loop over edges and faces and replace the + * corresponding vertices with @p pNew. */ void Element::SetVertex(unsigned int p, NodeSharedPtr pNew, bool descend) { @@ -125,8 +127,10 @@ void Element::SetVertex(unsigned int p, NodeSharedPtr pNew, bool descend) * When an edge is replaced, the element faces are also searched and * the corresponding face edges are updated to maintain consistency. * - * @param p Index of the edge to replace. - * @param pNew New edge. + * @param p Index of the edge to replace. + * @param pNew New edge. + * @param descend If true, we loop over faces and replace the corresponding + * face edge with @p pNew. */ void Element::SetEdge(unsigned int p, EdgeSharedPtr pNew, bool descend) { diff --git a/library/NekMeshUtils/MeshElements/Element.h b/library/NekMeshUtils/MeshElements/Element.h index 4cc60cb98..03e678473 100644 --- a/library/NekMeshUtils/MeshElements/Element.h +++ b/library/NekMeshUtils/MeshElements/Element.h @@ -244,9 +244,11 @@ public: m_id = p; } /// Replace a vertex with another vertex object. - NEKMESHUTILS_EXPORT void SetVertex(unsigned int p, NodeSharedPtr pNew, bool descend = true); + NEKMESHUTILS_EXPORT void SetVertex( + unsigned int p, NodeSharedPtr pNew, bool descend = true); /// Replace an edge with another edge object. - NEKMESHUTILS_EXPORT void SetEdge(unsigned int p, EdgeSharedPtr pNew, bool descend = true); + NEKMESHUTILS_EXPORT void SetEdge( + unsigned int p, EdgeSharedPtr pNew, bool descend = true); /// Replace a face with another face object. NEKMESHUTILS_EXPORT void SetFace(unsigned int p, FaceSharedPtr pNew); /// Set a correspondence between this element and an edge @@ -451,7 +453,23 @@ public: NEKMESHUTILS_EXPORT int GetMaxOrder(); - /// Insert interior (i.e. volume) points into this element. + /** + * @brief Insert interior (i.e. volume) points into this element to make the + * geometry an order @p order representation. + * + * @param order The desired polynomial order. + * @param geom The geometry object used to describe the curvature + * mapping. + * @param edgeType The points distribution to use on the volume. + * @param coordDim The coordinate (i.e. space) dimension. + * @param id Counter which should be incremented to supply + * consistent vertex IDs. + * @param justConfig If true, then the configuration Element::m_conf + * will be updated but no nodes will be + * generated. This is used when considering boundary + * elements, which just require copying of face or + * edge interior nodes. + */ NEKMESHUTILS_EXPORT virtual void MakeOrder( int order, SpatialDomains::GeometrySharedPtr geom, @@ -464,6 +482,10 @@ public: "This function should be implemented at a shape level."); } + /** + * @brief Get the edge orientation of @p edge with respect to the local + * element, which lies at edge index @p edgeId. + */ NEKMESHUTILS_EXPORT virtual StdRegions::Orientation GetEdgeOrient( int edgeId, EdgeSharedPtr edge) { @@ -472,6 +494,9 @@ public: return StdRegions::eNoOrientation; } + /** + * @brief Returns the local index of vertex @p j of face @p i. + */ NEKMESHUTILS_EXPORT virtual int GetFaceVertex(int i, int j) { ASSERTL0(false, diff --git a/library/NekMeshUtils/MeshElements/Face.h b/library/NekMeshUtils/MeshElements/Face.h index 6cfa27bbf..944cd3967 100644 --- a/library/NekMeshUtils/MeshElements/Face.h +++ b/library/NekMeshUtils/MeshElements/Face.h @@ -111,7 +111,8 @@ public: } /// Assemble a list of nodes on curved face - NEKMESHUTILS_EXPORT void GetCurvedNodes(std::vector &nodeList) const + NEKMESHUTILS_EXPORT void GetCurvedNodes( + std::vector &nodeList) const { // Treat 2D point distributions differently to 3D. if (m_curveType == LibUtilities::eNodalTriFekete || @@ -215,6 +216,7 @@ public: return s.str(); } + /// Make this face an order @p order face. @see Element::MakeOrder. void MakeOrder(int order, SpatialDomains::GeometrySharedPtr geom, LibUtilities::PointsType pType, diff --git a/library/NekMeshUtils/MeshElements/Quadrilateral.h b/library/NekMeshUtils/MeshElements/Quadrilateral.h index 620a228cc..ebdbbce33 100644 --- a/library/NekMeshUtils/MeshElements/Quadrilateral.h +++ b/library/NekMeshUtils/MeshElements/Quadrilateral.h @@ -114,7 +114,7 @@ template struct HOQuadrilateral void ReverseX() { - int np = (int)(sqrt(surfVerts.size()) + 0.5); + int np = (int)(sqrt((NekDouble)surfVerts.size()) + 0.5); for (int i = 0; i < np; ++i) { for (int j = 0; j < np/2; ++j) -- GitLab From c755032ea539857240d351b0861449b55f957a93 Mon Sep 17 00:00:00 2001 From: David Moxey Date: Wed, 17 Aug 2016 09:54:42 +0100 Subject: [PATCH 25/25] Fixing more Windows bugs --- library/NekMeshUtils/MeshElements/Quadrilateral.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/library/NekMeshUtils/MeshElements/Quadrilateral.h b/library/NekMeshUtils/MeshElements/Quadrilateral.h index ebdbbce33..560a07d76 100644 --- a/library/NekMeshUtils/MeshElements/Quadrilateral.h +++ b/library/NekMeshUtils/MeshElements/Quadrilateral.h @@ -126,7 +126,7 @@ template struct HOQuadrilateral void ReverseY() { - int np = (int)(sqrt(surfVerts.size()) + 0.5); + int np = (int)(sqrt((NekDouble)surfVerts.size()) + 0.5); // Reverse y direction for (int j = 0; j < np; ++j) { @@ -139,7 +139,7 @@ template struct HOQuadrilateral void Transpose() { - int np = (int)(sqrt(surfVerts.size()) + 0.5); + int np = (int)(sqrt((NekDouble)surfVerts.size()) + 0.5); std::vector tmp(surfVerts.size()); for (int i = 0; i < np; ++i) -- GitLab