Commit d199ea62 authored by Douglas Serson's avatar Douglas Serson

Merge branch 'fix/output-gmsh' into 'master'

Add support for Gmsh high-order output

This MR primarily adds support for Gmsh output at arbitrarily high-order for quads, hexes, triangles and tetrahedra, and up to order 4 output for prisms.

See merge request !679
parents bff54bcf c755032e
......@@ -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
......
......@@ -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}
\label{s:utilities:nekmesh:extract}
......
......@@ -38,6 +38,7 @@
#include <iomanip>
#include <LibUtilities/Foundations/ManagerAccess.h>
#include <SpatialDomains/SegGeom.h>
#include <NekMeshUtils/NekMeshUtilsDeclspec.h>
......@@ -169,6 +170,45 @@ public:
return ret;
}
/// Make this edge an order @p order edge. @see Element::MakeOrder.
void MakeOrder(int order,
SpatialDomains::GeometrySharedPtr geom,
LibUtilities::PointsType edgeType,
int coordDim,
int &id)
{
int nPoints = order + 1;
StdRegions::StdExpansionSharedPtr xmap = geom->GetXmap();
Array<OneD, NekDouble> edgePoints;
LibUtilities::PointsKey edgeKey(nPoints, edgeType);
LibUtilities::PointsManager()[edgeKey]->GetPoints(edgePoints);
Array<OneD, Array<OneD, NekDouble> > phys(coordDim);
for (int i = 0; i < coordDim; ++i)
{
phys[i] = Array<OneD, NekDouble>(xmap->GetTotPoints());
xmap->BwdTrans(geom->GetCoeffs(i), phys[i]);
}
m_edgeNodes.resize(nPoints - 2);
for (int i = 1; i < nPoints - 1; ++i)
{
Array<OneD, NekDouble> x(3, 0.0);
for (int j = 0; j < coordDim; ++j)
{
x[j] = xmap->PhysEvaluate(edgePoints + i, phys[j]);
}
m_edgeNodes[i-1] = boost::shared_ptr<Node>(
new Node(id++, x[0], x[1], x[2]));
}
m_curveType = edgeType;
}
/// ID of edge.
unsigned int m_id;
/// First vertex node.
......
......@@ -71,13 +71,21 @@ 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)
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)
......@@ -119,13 +127,21 @@ void Element::SetVertex(unsigned int p, NodeSharedPtr pNew)
* 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)
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)
......
......@@ -244,9 +244,11 @@ 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
......@@ -448,13 +450,60 @@ public:
"This function should be implemented at a shape level.");
return boost::shared_ptr<SpatialDomains::Geometry>();
}
NEKMESHUTILS_EXPORT int GetMaxOrder();
/// Complete this object.
NEKMESHUTILS_EXPORT virtual void Complete(int order)
/**
* @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,
LibUtilities::PointsType edgeType,
int coordDim,
int &id,
bool justConfig = false)
{
ASSERTL0(false,
"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)
{
ASSERTL0(false,
"This function should be implemented at a shape level.");
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,
"This function should be implemented at a shape level.");
return 0;
}
NEKMESHUTILS_EXPORT void Print()
{
int i, j;
......
......@@ -111,7 +111,8 @@ public:
}
/// Assemble a list of nodes on curved face
NEKMESHUTILS_EXPORT void GetCurvedNodes(std::vector<NodeSharedPtr> &nodeList) const
NEKMESHUTILS_EXPORT void GetCurvedNodes(
std::vector<NodeSharedPtr> &nodeList) const
{
// Treat 2D point distributions differently to 3D.
if (m_curveType == LibUtilities::eNodalTriFekete ||
......@@ -215,6 +216,114 @@ 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,
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<OneD, NekDouble> px, py;
LibUtilities::PointsKey pKey(nPoints, pType);
ASSERTL1(pKey.GetPointsDim() == 2, "Points distribution must be 2D");
LibUtilities::PointsManager()[pKey]->GetPoints(px, py);
Array<OneD, Array<OneD, NekDouble> > phys(coordDim);
for (int i = 0; i < coordDim; ++i)
{
phys[i] = Array<OneD, NekDouble>(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<OneD, NekDouble> xp(2);
xp[0] = px[i];
xp[1] = py[i];
Array<OneD, NekDouble> x(3, 0.0);
for (int j = 0; j < coordDim; ++j)
{
x[j] = xmap->PhysEvaluate(xp, phys[j]);
}
m_faceNodes[cnt] = boost::shared_ptr<Node>(
new Node(id++, x[0], x[1], x[2]));
}
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<OneD, NekDouble> px;
LibUtilities::PointsKey pKey(nPoints, pType);
ASSERTL1(pKey.GetPointsDim() == 1, "Points distribution must be 1D");
LibUtilities::PointsManager()[pKey]->GetPoints(px);
Array<OneD, Array<OneD, NekDouble> > phys(coordDim);
for (int i = 0; i < coordDim; ++i)
{
phys[i] = Array<OneD, NekDouble>(xmap->GetTotPoints());
xmap->BwdTrans(geom->GetCoeffs(i), phys[i]);
}
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<OneD, NekDouble> xp(2);
xp[0] = px[j];
xp[1] = px[i];
Array<OneD, NekDouble> x(3, 0.0);
for (int k = 0; k < coordDim; ++k)
{
x[k] = xmap->PhysEvaluate(xp, phys[k]);
}
m_faceNodes[cnt] = boost::shared_ptr<Node>(
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.
NEKMESHUTILS_EXPORT SpatialDomains::Geometry2DSharedPtr GetGeom(int coordDim)
......
......@@ -47,6 +47,12 @@ LibUtilities::ShapeType Hexahedron::m_type =
GetElementFactory().RegisterCreatorFunction(
LibUtilities::eHexahedron, Hexahedron::create, "Hexahedron");
/// Vertex IDs that make up hexahedron faces.
int Hexahedron::m_faceIds[6][4] = {
{0, 1, 2, 3}, {0, 1, 5, 4}, {1, 2, 6, 5},
{3, 2, 6, 7}, {0, 3, 7, 4}, {4, 5, 6, 7}
};
/**
* @brief Create a hexahedral element.
*/
......@@ -102,12 +108,6 @@ Hexahedron::Hexahedron(ElmtConfig pConf,
// Create faces
int face_edges[6][4];
int face_ids[6][4] = {{0, 1, 2, 3},
{0, 1, 5, 4},
{1, 2, 6, 5},
{3, 2, 6, 7},
{0, 3, 7, 4},
{4, 5, 6, 7}};
for (int j = 0; j < 6; ++j)
{
vector<NodeSharedPtr> faceVertices;
......@@ -115,9 +115,9 @@ Hexahedron::Hexahedron(ElmtConfig pConf,
vector<NodeSharedPtr> faceNodes;
for (int k = 0; k < 4; ++k)
{
faceVertices.push_back(m_vertex[face_ids[j][k]]);
NodeSharedPtr a = m_vertex[face_ids[j][k]];
NodeSharedPtr b = m_vertex[face_ids[j][(k + 1) % 4]];
faceVertices.push_back(m_vertex[m_faceIds[j][k]]);
NodeSharedPtr a = m_vertex[m_faceIds[j][k]];
NodeSharedPtr b = m_vertex[m_faceIds[j][(k + 1) % 4]];
for (unsigned int i = 0; i < m_edge.size(); ++i)
{
if (((*(m_edge[i]->m_n1) == *a) &&
......@@ -176,6 +176,96 @@ 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,
bool justConfig)
{
m_curveType = pType;
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;
if (justConfig)
{
return;
}
int nPoints = order + 1;
StdRegions::StdExpansionSharedPtr xmap = geom->GetXmap();
Array<OneD, NekDouble> px;
LibUtilities::PointsKey pKey(nPoints, pType);
ASSERTL1(pKey.GetPointsDim() == 1, "Points distribution must be 1D");
LibUtilities::PointsManager()[pKey]->GetPoints(px);
Array<OneD, Array<OneD, NekDouble> > phys(coordDim);
for (int i = 0; i < coordDim; ++i)
{
phys[i] = Array<OneD, NekDouble>(xmap->GetTotPoints());
xmap->BwdTrans(geom->GetCoeffs(i), phys[i]);
}
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<OneD, NekDouble> xp(3);
xp[0] = px[k];
xp[1] = px[j];
xp[2] = px[i];
Array<OneD, NekDouble> x(3, 0.0);
for (int k = 0; k < coordDim; ++k)
{
x[k] = xmap->PhysEvaluate(xp, phys[k]);
}
m_volumeNodes[cnt] = boost::shared_ptr<Node>(
new Node(id++, x[0], x[1], x[2]));
}
}
}
}
/**
* @brief Return the number of nodes defining a hexahedron.
*/
......
......@@ -77,8 +77,24 @@ 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,
bool justConfig = false);
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];
};
}
}
......
......@@ -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<OneD, NekDouble> px;
LibUtilities::PointsKey pKey(nPoints, pType);
ASSERTL1(pKey.GetPointsDim() == 1, "Points distribution must be 1D");
LibUtilities::PointsManager()[pKey]->GetPoints(px);
Array<OneD, Array<OneD, NekDouble> > phys(coordDim);
for (int i = 0; i < coordDim; ++i)
{
phys[i] = Array<OneD, NekDouble>(xmap->GetTotPoints());
xmap->BwdTrans(geom->GetCoeffs(i), phys[i]);
}
int nQuadIntPts = (nPoints - 2) * (nPoints - 2);
m_volumeNodes.resize(nQuadIntPts);
for (int i = 1, cnt = 0; i < nPoints-1; ++i)
{
Array<OneD, NekDouble> xp(1);
xp[0] = px[i];
Array<OneD, NekDouble> x(3, 0.0);
for (int k = 0; k < coordDim; ++k)
{
x[k] = xmap->PhysEvaluate(xp, phys[k]);
}
m_volumeNodes[cnt] = boost::shared_ptr<Node>(
new Node(id++, x[0], x[1], x[2]));
}
}
/**
* @brief Return the number of nodes defining a line.
*/
......
......@@ -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);
};
}
......
......@@ -34,6 +34,7 @@
////////////////////////////////////////////////////////////////////////////////
#include <NekMeshUtils/MeshElements/Mesh.h>
#include <LibUtilities/Foundations/ManagerAccess.h>
using namespace std;
......@@ -79,5 +80,167 @@ unsigned int Mesh::GetNumEntities()
return nEnt;
}
/**
* @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)
{
int id = m_vertexSet.size();
EdgeSet::iterator eit;
FaceSet::iterator fit;
boost::unordered_map<int, SpatialDomains::Geometry1DSharedPtr> edgeGeoms;
boost::unordered_map<int, SpatialDomains::Geometry2DSharedPtr> faceGeoms;
boost::unordered_map<int, SpatialDomains::GeometrySharedPtr> volGeoms;
// Decide on distribution of points to use for each shape type based on the
// input we've been supplied.
std::map<LibUtilities::ShapeType, LibUtilities::PointsType> pTypes;
if (distType == LibUtilities::ePolyEvenlySpaced)
{
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;
}
else if (distType == LibUtilities::eGaussLobattoLegendre)
{
// Prism still to do.
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;
}
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
// we process each element.
for(eit = m_edgeSet.begin(); eit != m_edgeSet.end(); eit++)
{
SpatialDomains::Geometry1DSharedPtr geom =