Commit b86b5337 authored by Chris Cantwell's avatar Chris Cantwell
Browse files

Merge branch 'feature/mc-hex-bl' into 'master'

Various MeshConvert additions

This MR adds various features and corrects quite a few bugs present in `MeshConvert`:

*Features*:
- Ability to split hexahedra in the boundary layer refinement module.
- New module to add curvature defined by a scalar function z=f(x,y).
- Improved Gmsh support for high-order elements.
- A small test suite of relevant examples which we can build on later.

*Fixes*:
- Many fixes relating to face-interior curvature.
- All `MeshConvert` elements are now orientated exactly in the same manner as their Nektar++ counterparts. (This will allow for some speed improvements in a future branch.)
- Lots of small bugs in `InputNek` relating to high-order curvature.

See merge request !401
parents 57e8235b 5d3ac263
......@@ -133,7 +133,7 @@ solve the resulting linear system.
\mc is designed to provide a pipeline approach to mesh generation. To do this,
we break up tasks into three different types. Each task is called a
\emph{module} and a chain of modules specifies the pipeline. There are
\emph{module} and a chain of modules specifies the pipeline.
%
\begin{itemize}
\item \textbf{Input} modules read meshes in a variety of formats;
......@@ -200,14 +200,15 @@ line argument:
\subsubsection{Input modules}
Input and output modules use file extension names to determine the correct
module to use. The table below indicates support currently implemented.
module to use. Not every module is capable of reading high-order information,
where it exists. The table below indicates support currently implemented.
\begin{center}
\begin{tabularx}{\linewidth}{llcX}
\toprule
\textbf{Format} & \textbf{Extension} & \textbf{High-order} & \textbf{Notes}\\
\midrule
Gmsh & \texttt{msh} & \cmark & Only reads elements and physical groups. Only edge curvature currently imported.\\
Gmsh & \texttt{msh} & \cmark & Only reads nodes, elements and physical groups (which are mapped to composites).\\
Nektar & \texttt{rea} & \cmark & Reads elements, fluid boundary conditions. Most curve types are unsupported: high-order information must be defined in an accompanying .hsf file. \\
Nektar++ & \texttt{xml} & \cmark & Fully supported. \\
PLY & \texttt{ply} & \xmark & Reads only the ASCII format.. \\
......@@ -226,6 +227,24 @@ Note that you can override the module used on the command line. For example,
MeshConvert pipe-3d:sem pipe-3d.xml
\end{lstlisting}
Typically, mesh generators allow physical surfaces and volumes to contain many
element types; for example a cube could be constructed from a mixture of hexes
and prisms. In \nekpp, a composite can only contain a single element
type. Whilst the converter will attempt to preserve the numbering of composites
from the original mesh type, sometimes a renumbering will occur when a domain
contains many element types. For example, for a domain with the tag \inltt{150}
containing quadrilaterals and triangles, the Gmsh reader will print a
notification along the lines of:
\begin{lstlisting}[style=BashInputStyle]
Multiple elements in composite detected; remapped:
- Tag 150 => 150 (Triangle), 151 (Quadrilateral)
\end{lstlisting}
The resulting file therefore has two composites of IDs \inltt{150} and
\inltt{151} respectively, containing the triangular and quadrilateral elements
of the original mesh.
\subsubsection{Output modules}
The following output formats are supported:
......@@ -251,17 +270,28 @@ available within \mc.
\subsubsection{Negative Jacobian detection}
To detect elements with negative Jacobians use the \inltt{jac} module:
To detect elements with negative Jacobian determinant, use the \inltt{jac}
module:
%
\begin{lstlisting}[style=BashInputStyle]
MeshConvert -m jac Mesh.xml output.xml
\end{lstlisting}
%
To get a detailed list of elements which have negative Jacobians, one may use
the \inltt{list} option:
%
\begin{lstlisting}[style=BashInputStyle]
MeshConvert -m jac:list Mesh.xml output.xml
\end{lstlisting}
%
and to extract the elements for the purposes of visualisation within the domain,
use the \inltt{extract} boolean parameter:
%
\begin{lstlisting}[style=BashInputStyle]
MeshConvert -m jac:extract Mesh.xml MeshWithNegativeElements.xml
\end{lstlisting}
\subsubsection{Spherigon Patches}
\subsubsection{Spherigon patches}
Where high-order information is not available (e.g. when using meshes from
imaging software), various techniques can be used to apply a smoothing to the
......@@ -356,13 +386,13 @@ MeshConvert -m peralign:surf1=11:surf2=12:dir=y:orient input.dat output.xml
\inltt{peralign} module before the \inltt{spherigon} module.
\end{notebox}
\subsubsection{Prism splitting}
\subsubsection{Boundary layer splitting}
Often it is the case that one can generate a coarse boundary layer grid of a
mesh. \mc has a method for splitting prismatic elements into finer elements
based on the work presented in~\cite{MoHaJoSh14} and~\cite{MoHaPeSh14b}. You
must have a prismatic mesh that is $O$-type -- that is, you can modify the
prismatic layer without modifying the rest of the mesh.
mesh. \mc has a method for splitting prismatic and hexahedral elements into
finer elements based on the work presented in~\cite{MoHaJoSh14}
and~\cite{MoHaPeSh14b}. You must have a prismatic mesh that is $O$-type -- that
is, you can modify the boundary layer without modifying the rest of the mesh.
Given $n$ layers, and a ratio $r$ which defines the relative heights of elements
in different layers, the method works by defining a geometric progression of
......@@ -370,9 +400,9 @@ points
\[
x_k = x_{k-1} + ar^k, \quad a = \frac{2(1-r)}{1 - r^{n+1}}
\]
in the standard segment $[-1,1]$. These are then projected into the coarse prism
to construct a sequence of increasingly refined elements, as depicted in
figure~\ref{fig:util:mc:split}.
in the standard segment $[-1,1]$. These are then projected into the coarse
elements to construct a sequence of increasingly refined elements, as depicted
in figure~\ref{fig:util:mc:split}.
\begin{figure}
\begin{center}
......@@ -409,7 +439,9 @@ of 2 and 7 integration points per element use the following command:
\begin{notebox}
You can also use an expression in terms of coordinates $(x,y,z)$ for $r$ to
make the ratio spatially varying; e.g. \inltt{r=x\^2+y\^2}.
make the ratio spatially varying; e.g. \inltt{r=sin(x)}. In this case the
function should be sufficiently smooth to prevent the elements
self-intersecting.
\end{notebox}
\subsubsection{High-order cylinder generation}
......@@ -466,6 +498,24 @@ module:
MeshConvert -m detect volume.xml volumeWithBoundaryComposite.xml
\end{lstlisting}
\subsubsection{Scalar function curvature}
This module imposes curvature on a surface given a scalar function
$z=f(x,y)$. For example, if on surface 1 we wish to apply a surface defined by a
Gaussian $z = \exp[-(x^2+y^2)]$ using 7 quadrature points in each direction, we
may issue the command
\begin{lstlisting}[style=BashInputStyle]
MeshConvert -m scalar:surf=1:nq=7:scalar=exp\(x*x+y*y\) mesh.xml deformed.xml
\end{lstlisting}
\begin{notebox}
This module makes no attempt to apply the curvature to the interior of the
domain. Elements must therefore be coarse in order to prevent
self-intersection. If a boundary layer is required, one option is to use this
module in combination with the splitting module described earlier.
\end{notebox}
%%% Local Variables:
%%% mode: latex
%%% TeX-master: "../user-guide"
......
......@@ -16,6 +16,7 @@ SET(MeshConvertHeaders
ProcessExtractSurf.h
ProcessJac.h
ProcessPerAlign.h
ProcessScalar.h
ProcessSpherigon.h
ProcessTetSplit.h
)
......@@ -39,6 +40,7 @@ SET(MeshConvertSources
ProcessExtractSurf.cpp
ProcessJac.cpp
ProcessPerAlign.cpp
ProcessScalar.cpp
ProcessSpherigon.cpp
ProcessTetSplit.cpp
)
......@@ -57,3 +59,22 @@ IF (NEKTAR_USE_VTK)
TARGET_LINK_LIBRARIES(MeshConvert vtkCommonCore vtkIOLegacy)
ENDIF ()
ENDIF (NEKTAR_USE_VTK)
# Gmsh tests
ADD_NEKTAR_TEST (Gmsh/CubeAllElements)
ADD_NEKTAR_TEST (Gmsh/CubeHex)
ADD_NEKTAR_TEST (Gmsh/CubePrism)
ADD_NEKTAR_TEST (Gmsh/CubeTet)
IF (WIN32)
ADD_NEKTAR_TEST (Gmsh/Scalar_Windows)
ELSE ()
ADD_NEKTAR_TEST (Gmsh/Scalar)
ENDIF ()
ADD_NEKTAR_TEST (Gmsh/SquareQuad)
ADD_NEKTAR_TEST (Gmsh/SquareTri)
# Nektar tests
ADD_NEKTAR_TEST (Nektar/HexLinear)
ADD_NEKTAR_TEST (Nektar/Tube45)
# StarCCM tests
ADD_NEKTAR_TEST (StarTec/CubePer)
ADD_NEKTAR_TEST_LENGTHY(StarTec/StraightRW)
......@@ -64,9 +64,15 @@ namespace Nektar
* Element map; takes a msh id to an %ElmtConfig object.
*/
static std::map<unsigned int, ElmtConfig> elmMap;
private:
int GetNnodes(unsigned int InputGmshEntity);
int GetNnodes(unsigned int InputGmshEntity);
vector<int> CreateReordering(unsigned int InputGmshEntity);
vector<int> TriReordering(ElmtConfig conf);
vector<int> QuadReordering(ElmtConfig conf);
vector<int> HexReordering(ElmtConfig conf);
vector<int> PrismReordering(ElmtConfig conf);
vector<int> TetReordering(ElmtConfig conf);
};
}
}
......
......@@ -33,7 +33,6 @@
//
////////////////////////////////////////////////////////////////////////////////
#include "MeshElements.h"
#include "InputNek.h"
......@@ -85,7 +84,7 @@ namespace Nektar
int i, j, k, nodeCounter = 0;
int nComposite = 0;
LibUtilities::ShapeType elType;
double vertex[3][6];
double vertex[3][8];
map<LibUtilities::ShapeType,int> domainComposite;
map<LibUtilities::ShapeType,vector< vector<NodeSharedPtr> > > elNodes;
map<LibUtilities::ShapeType,vector<int> > elIds;
......@@ -221,8 +220,7 @@ namespace Nektar
else if (line.find("Pyr") != string::npos ||
line.find("pyr") != string::npos)
{
cerr << "Pyramid elements not yet supported." << endl;
abort();
elType = LibUtilities::ePyramid;
}
else if (line.find("Qua") != string::npos ||
line.find("qua") != string::npos)
......@@ -269,7 +267,7 @@ namespace Nektar
vertex[1][k], vertex[2][k]));
nodeList.push_back(n);
}
elNodes[elType].push_back(nodeList);
elIds [elType].push_back(i);
}
......@@ -281,7 +279,7 @@ namespace Nektar
{
LibUtilities::ShapeType elType = elmOrder[i];
vector<vector<NodeSharedPtr> > &tmp = elNodes[elType];
for (j = 0; j < tmp.size(); ++j)
{
vector<int> tags;
......@@ -383,7 +381,7 @@ namespace Nektar
cerr << "Unable to read number of curved sides" << endl;
abort();
}
int nCurvedSides;
int faceId, elId;
map<string,pair<NekCurve, string> >::iterator it;
......@@ -403,6 +401,8 @@ namespace Nektar
elId = elMap[elId-1];
ElementSharedPtr el = m_mesh->m_element[m_mesh->m_expDim][elId];
int origFaceId = faceId;
if (el->GetConf().m_e == LibUtilities::ePrism && faceId % 2 == 0)
{
boost::shared_ptr<Prism> p =
......@@ -485,109 +485,76 @@ namespace Nektar
}
else if (it->second.first == eFile)
{
vector<unsigned int> vertId(3);
FaceSharedPtr f = el->GetFace(faceId);
static int tetFaceVerts[4][3] = {
{0,1,2},{0,1,3},{1,2,3},{0,2,3}};
vector<int> vertId(3);
s >> vertId[0] >> vertId[1] >> vertId[2];
// Find vertex combination in hoData.
hoIt = hoData[word].find(HOSurfSharedPtr(
new HOSurf(vertId)));
if (hoIt == hoData[word].end())
{
cerr << "Unable to find high-order surface data "
<< "for element id " << elId+1 << endl;
abort();
}
// Depending on order of vertices in rea file, surface
// information may need to be rotated or
// reflected. These procedures are taken from
// nektar/Hlib/src/HOSurf.C
HOSurfSharedPtr surf = *hoIt;
if (vertId[0] == surf->vertId[0])
{
if (vertId[1] == surf->vertId[1] ||
vertId[1] == surf->vertId[2])
{
if (vertId[1] == surf->vertId[2])
{
surf->Rotate(1);
surf->Reflect();
}
}
}
else if (vertId[0] == surf->vertId[1])
{
if (vertId[1] == surf->vertId[0] ||
vertId[1] == surf->vertId[2])
{
if (vertId[1] == surf->vertId[0])
{
surf->Reflect();
}
else
{
surf->Rotate(2);
}
}
}
else if (vertId[0] == surf->vertId[2])
// Prisms and tets may have been rotated by OrientPrism
// routine which reorders vertices. This block rotates
// vertex ids accordingly.
if (el->GetConf().m_e == LibUtilities::eTetrahedron)
{
if (vertId[1] == surf->vertId[0] ||
vertId[1] == surf->vertId[1])
boost::shared_ptr<Tetrahedron> tet =
boost::static_pointer_cast<Tetrahedron>(el);
vector<int> tmpVertId = vertId;
for (j = 0; j < 3; ++j)
{
if (vertId[1] == surf->vertId[1])
{
surf->Rotate(2);
surf->Reflect();
}
else
int v = tet->GetVertex(
tet->m_origVertMap[
tetFaceVerts[origFaceId][j]])->m_id;
for (k = 0; k < 3; ++k)
{
surf->Rotate(1);
int w = f->m_vertexList[k]->m_id;
if (v == w)
{
vertId[k] = tmpVertId[j];
break;
}
}
}
}
// If the element is a prism, check to see if
// orientation has changed and update order of surface
// vertices.
int reverseSide = 2;
// Prisms may have been rotated by OrientPrism routine
// and break curved faces. This block rotates faces
// accordingly. TODO: Add similar routine for tets.
if (el->GetConf().m_e == LibUtilities::ePrism)
else if (el->GetConf().m_e == LibUtilities::ePrism)
{
boost::shared_ptr<Prism> pr =
boost::shared_ptr<Prism> pr =
boost::static_pointer_cast<Prism>(el);
if (pr->m_orientation == 1)
{
// Prism has been rotated counter-clockwise;
// rotate face, reverse what was the last edge
// (now located at edge 0).
(*hoIt)->Rotate(1);
reverseSide = 0;
swap(vertId[2], vertId[1]);
swap(vertId[0], vertId[1]);
}
else if (pr->m_orientation == 2)
{
// Prism has been rotated clockwise; rotate
// face, reverse what was the last edge (now
// located at edge 1).
(*hoIt)->Rotate(2);
reverseSide = 1;
swap(vertId[0], vertId[1]);
swap(vertId[2], vertId[1]);
}
}
// Find vertex combination in hoData.
hoIt = hoData[word].find(HOSurfSharedPtr(
new HOSurf(vertId)));
// Finally, add high order data to appropriate
// face. NOTE: this is a bit of a hack since the
// elements are technically linear, but should work just
// fine.
FaceSharedPtr f = el->GetFace(faceId);
if (hoIt == hoData[word].end())
{
cerr << "Unable to find high-order surface data "
<< "for element id " << elId+1 << endl;
abort();
}
// Depending on order of vertices in rea file, surface
// information may need to be rotated or reflected.
HOSurfSharedPtr surf = *hoIt;
surf->Align(vertId);
// Finally, add high order data to appropriate face.
int Ntot = (*hoIt)->surfVerts.size();
int N = ((int)sqrt(8.0*Ntot+1.0)-1)/2;
EdgeSharedPtr edge;
// Apply high-order map to convert face data to Nektar++
// ordering (vertices->edges->internal).
vector<NodeSharedPtr> tmpVerts = (*hoIt)->surfVerts;
......@@ -600,36 +567,44 @@ namespace Nektar
{
NodeSharedPtr a = (*hoIt)->surfVerts[j];
}
vector<int> faceVertIds(3);
faceVertIds[0] = f->m_vertexList[0]->m_id;
faceVertIds[1] = f->m_vertexList[1]->m_id;
faceVertIds[2] = f->m_vertexList[2]->m_id;
for (j = 0; j < f->m_edgeList.size(); ++j)
{
edge = f->m_edgeList[j];
// Skip over edges which have already been
// populated, apart from those which need to be
// reoriented.
if (edge->m_edgeNodes.size() > 0 && reverseSide == 2)
// populated.
if (edge->m_edgeNodes.size() > 0)
{
continue;
}
edge->m_edgeNodes.clear();
edge->m_curveType = LibUtilities::eGaussLobattoLegendre;
//edge->m_edgeNodes.clear();
edge->m_curveType
= LibUtilities::eGaussLobattoLegendre;
for (int k = 0; k < N-2; ++k)
{
edge->m_edgeNodes.push_back(
(*hoIt)->surfVerts[3+j*(N-2)+k]);
}
// Reverse order of modes along correct side.
if (j == reverseSide)
// Nodal triangle data is always
// counter-clockwise. Add this check to reorder
// where necessary.
if (edge->m_n1->m_id != faceVertIds[j])
{
reverse(edge->m_edgeNodes.begin(),
edge->m_edgeNodes.end());
}
}
// Insert interior face curvature.
f->m_curveType = LibUtilities::eNodalTriElec;
for (int j = 3+3*(N-2); j < Ntot; ++j)
{
......@@ -1026,7 +1001,7 @@ namespace Nektar
getline(hsf, line);
// Read in nodal points for each face.
map<unsigned int, vector<NodeSharedPtr> > faceMap;
map<int, vector<NodeSharedPtr> > faceMap;
for (int i = 0; i < Nface; ++i)
{
getline(hsf, line);
......@@ -1052,9 +1027,9 @@ namespace Nektar
getline(hsf, line);
for (int i = 0; i < Nface; ++i)
{
string tmp;
int fid;
vector<unsigned int> nodeIds(3);
string tmp;
int fid;
vector<int> nodeIds(3);
getline(hsf, line);
ss.clear(); ss.str(line);
......@@ -1090,6 +1065,7 @@ namespace Nektar
case LibUtilities::eTriangle: nNodes = 3; break;
case LibUtilities::eQuadrilateral: nNodes = 4; break;
case LibUtilities::eTetrahedron: nNodes = 4; break;
case LibUtilities::ePyramid: nNodes = 5; break;
case LibUtilities::ePrism: nNodes = 6; break;
case LibUtilities::eHexahedron: nNodes = 8; break;
default:
......@@ -1111,9 +1087,9 @@ namespace Nektar
{
return false;
}
vector<unsigned int> ids1 = p1->vertId;
vector<unsigned int> ids2 = p2->vertId;
vector<int> ids1 = p1->vertId;
vector<int> ids2 = p2->vertId;
sort(ids1.begin(), ids1.end());
sort(ids2.begin(), ids2.end());
......@@ -1125,64 +1101,5 @@ namespace Nektar
return true;
}
/**
* Rotates the triangle of data points inside #surfVerts
* counter-clockwise nrot times.
*
* @param nrot Number of times to rotate triangle.
*/
void HOSurf::Rotate(int nrot)
{
int n, i, j, cnt;
int np = ((int)sqrt(8.0*surfVerts.size()+1.0)-1)/2;
NodeSharedPtr* tmp = new NodeSharedPtr[np*np];
//NodeSharedPtr tmp[np][np];
for (n = 0; n < nrot; ++n)
{
for (cnt = i = 0; i < np; ++i)
{
for (j = 0; j < np-i; ++j, cnt++)
{
tmp[i*np+j] = surfVerts[cnt];
}
}
for (cnt = i = 0; i < np; ++i)
{
for (j = 0; j < np-i; ++j,cnt++)
{
surfVerts[cnt] = tmp[(np-1-i-j)*np+i];
}
}
}
delete[] tmp;
}
void HOSurf::Reflect()
{
int i, j, cnt;
int np = ((int)sqrt(8.0*surfVerts.size()+1.0)-1)/2;
NodeSharedPtr* tmp = new NodeSharedPtr[np*np];
for (cnt = i = 0; i < np; ++i)
{
for (j = 0; j < np-i; ++j,cnt++)
{
tmp[i*np+np-i-1-j] = surfVerts[cnt];
}
}
for(cnt = i = 0; i < np; ++i)
{
for(j = 0; j < np-i; ++j,cnt++)
{
surfVerts[cnt] = tmp[i*np+j];
}
}
delete[] tmp;
}
}
}
......@@ -42,42 +42,15 @@ namespace Nektar
{
namespace Utilities
{
/**
* Wrapper structure for high-order mesh information stored in a hsf
* file.
*/
struct HOSurf