\chapter{NekMesh}
\label{s:utilities:nekmesh}

\newcommand{\nm}{\texttt{NekMesh}\xspace}
\newcommand{\gmsh}{\texttt{Gmsh}\xspace}

\nm is a utility bundled with \nekpp which has two purposes:
\begin{itemize}
  \item allow foreign mesh file formats to be converted into \nekpp's XML
  format;
  \item aide in the generation of high-order meshes through a series of supplied
  processing modules.
\end{itemize}

\begin{notebox}
  \nm replaces a previous utility called MeshConvert. This change is to reflect
  the fact that the program no longer only converts and manipulates meshes but
  can now also generate them from a CAD definition. This mesh generator is in
  an early stage of development and as such is disabled by default. For the
  time being those not using the mesh generator can use \nm as they would
  have used MeshConvert, none of the functionality or methodology has changed.
\end{notebox}

There is also some limited support for other output formats. We begin by running
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 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 \inltt{.xml} file, i.e. \\ \inltt{
    \nm file.msh newfile.xml:xml:uncompress}
\end{notebox}


\section{Exporting a mesh from \gmsh}

To demonstrate how \nm works, we will define a simple channel-like 3D geometry.
First, we must define the \gmsh geometry to be used. The \gmsh definition is
given below, and is visualised in figure~\ref{fig:util:mc:gmsh-example}.

\begin{lstlisting}[style=XmlStyle]
Point(1) = {-1, 0, 0, 1.0};
Point(2) = {-0.3, 0, 0, 1.0};
Line(3) = {1, 2};
s[] = Extrude {0, 0, 7} {
  Line{3}; Layers{5}; Recombine;
};
v[] = Extrude {{0, 0, 1}, {0, 0, 0}, Pi} {
  Surface{s[1]}; Layers{10}; Recombine;
};
\end{lstlisting}

Whilst a full tutorial on \gmsh is far beyond the scope of this document, note
the use of the \texttt{Recombine} argument. This allows us to generate a
structured hexahedral mesh; remove the first \texttt{Recombine} to generate a
prismatic mesh and both occurances to generate a tetrahedral mesh. Increasing
the \texttt{Layers} numbers refines the mesh in the radial and azimuthal
direction respectively.

\section{Defining physical surfaces and volumes}

\begin{figure}
  \begin{center}
    \includegraphics[width=0.4\textwidth]{img/mc-example-gmsh}
    \includegraphics[width=0.4\textwidth]{img/mc-example-mesh}
  \end{center}
  \caption{Geometry definition in Gmsh (left) and resulting high-order mesh
    visualised in ParaView (right).}
  \label{fig:util:mc:gmsh-example}
\end{figure}

In order for us to use the mesh, we need to define the physical surfaces which
correspond to the inflow, outflow and walls so that we can set appropriate
boundary conditions. The numbering resulting from the extrusions in this case is
not straightforward. In the graphical interface, select \inlsh{Geometry >
  Physical Groups > Add > Surface}, and then hover over each of the surfaces
which are shown by the dashed gray lines. The numbering will be revealed in the
toolbar underneath the geometry as a ruled surface. In this case:
%
\begin{itemize}
\item \textbf{Walls:} surfaces 7, 8, 28, 29.
\item \textbf{Inflow:} surface 16.
\item \textbf{Outflow:} surface 24.
\end{itemize}
%
We also need to define the physical volumes, which can be done in a similar
fashion. For this example, there is only one volume having ID 1. Adding these
groups to the end of the \texttt{.geo} file is very straightforward:

\begin{lstlisting}[style=XmlStyle]
Physical Volume(0) = {1};
Physical Surface(1)= {7,8,28,29};
Physical Surface(2) = {16};
Physical Surface(3) = {24};
\end{lstlisting}
Either choose the option \inlsh{File->Save Mesh} or, assuming this is saved in
a file named \inlsh{test.geo}, run the command
\begin{lstlisting}[style=BashInputStyle]
gmsh -3 test.geo
\end{lstlisting}
which will produce the resulting MSH file \inlsh{test.msh}. One can generate a
high-order mesh by specifying the order on the command line, for example
\begin{lstlisting}[style=BashInputStyle]
gmsh -3 -order 6 test.geo
\end{lstlisting}
will generate a sixth-order mesh. Note that you will need to use a current
version of \gmsh in order to do this, most likely from subversion.

\section{Converting the MSH to Nektar++ format}
Assuming that you have compiled \nekpp according to the compilation
instructions, run the command
%
\begin{lstlisting}[style=BashInputStyle]
NekMesh test.msh test.xml
\end{lstlisting}
%
to generate the XML file.
%
\begin{notebox}
  This file contains only the geometry definition (and a default
  \inltt{EXPANSIONS} definition). In order to use this mesh, a
  \inltt{CONDITIONS} section must be supplied detailing the solver and
  parameters to use.
\end{notebox}
%
To validate the mesh visually, we can use a utility such as Paraview or
VisIt. To do this, we can use the \inltt{FieldConvert} command using:
%
\begin{lstlisting}[style=BashInputStyle]
FieldConvert test.xml test.vtu
\end{lstlisting}
%
which generates an unstructured VTK file \inlsh{test.vtu}.

It is possible that, when the high-order information was inserted into the mesh
by \gmsh, invalid elements are generated which self intersect. In this case, the
Jacobian of the mapping defining the curvature will have negative regions, which
will generate warnings such as:
\begin{lstlisting}[style=BashInputStyle]
Warning: Level 0 assertion violation
3D deformed Jacobian not positive (element ID = 48) (first vertex ID = 105)
\end{lstlisting}
This tells you the element ID that is invalid, and the ID of the first vertex of
the element. Whilst a resulting simulation may run, the results may not be valid
because of this problem, or excessively large amounts of time may be needed to
solve the resulting linear system.


\section{NekMesh in NekPy}

The Python interface allows the user to instantiate input,
output, and process modules by calling the static \verb+Create+ method of the
\verb+InputModule+, \verb+ProcessModule+, and \verb+OutputModule+,
register configuration options, and process them.
Consider the following example:

\begin{lstlisting}[style=C++Style, language=Python]
import sys
from NekPy.NekMesh import Mesh, ProcessModule, OutputModule

mesh          = Mesh()
mesh.expDim   = 3
mesh.spaceDim = 3
mesh.nummode  = 5
mesh.verbose  = True

# Load the CAD file
ProcessModule.Create("loadcad", mesh, \
                     filename="input.stp", verbose=True).Process()
# Load the octree
ProcessModule.Create("loadoctree", mesh, mindel=0.04,\
                      maxdel=0.2, eps=0.02).Process()
# Create a surface mesh
ProcessModule.Create("surfacemesh", mesh).Process()
# Output a 2D manifold mesh
mesh.expDim = 2
# Create a high-order surface
ProcessModule.Create("hosurface", mesh).Process()
# Dump out elemental Jacobians
ProcessModule.Create("jac", mesh, list=True).Process()
# Dump out the surface mesh.
OutputModule.Create("xml", mesh, test=True,\
                     outfile="output.xml").Process()
\end{lstlisting}

After importing the \verb+Mesh+, \verb+ProcessModule+, and \verb+OutputModule+
classes, first we create a \verb+Mesh+ object by calling the constructor.
This object will be shared by the modules.
Then we manipulate the \verb+Mesh+ object by creating different
\verb+ProcessModule+s, at the end we write out the result into a xml file
using an \verb+OutputModule+.
The configuration options for a given module are passed to the static
\verb+Create+ method of the \verb+InputModule+, \verb+ProcessModule+, and
\verb+OutputModule+. This creates the corresponding module and the
modules can be processed immediately after instantiation. Note that the first
parameter of the \verb+Create+ method has to be the key for a given module,
the second is the previously created \verb+Mesh+ object. The remaining keyword
arguments can specify additional parameters for a module.

The Python interface allows the user to create new modules by inheriting from
one of the possible base classes
(\verb+InputModule+, \verb+ProcessModule+, \verb+OutputModule+).

The following is a simple example when we inherit from the \verb+InputModule+
and override the \verb+Process+ method:
\begin{lstlisting}[style=C++Style, language=Python]
import sys
from NekPy.LibUtilities import ShapeType
import NekPy.NekMesh as NekMesh
import numpy as np

# StructuredGrid creates a 2D structured grid of triangles.
class StructuredGrid(NekMesh.InputModule):
    def __init__(self, mesh):
        super(StructuredGrid, self).__init__(mesh)
        self.mesh.spaceDim = 2
        self.mesh.expDim = 2
        # Define some configuration options for this module.
        self.AddConfigOption("nx", "2", "Number of points in x direction")
        self.AddConfigOption("ny", "2", "Number of points in y direction")
        self.AddConfigOption("lx", "0", "Lower-left x-coordinate")
        self.AddConfigOption("rx", "0", "Upper-right x-coordinate")
        self.AddConfigOption("ly", "0", "Lower-left y-coordinate")
        self.AddConfigOption("ry", "0", "Upper-right y-coordinate")
        self.AddConfigOption("compid", "0", "Composite ID")

    def Process(self):
        # Get the input variables from our configuration options.
        coord_1x   = self.GetFloatConfig("lx")
        coord_1y   = self.GetFloatConfig("ly")
        coord_2x   = self.GetFloatConfig("rx")
        coord_2y   = self.GetFloatConfig("ry")
        nx         = self.GetIntConfig("nx")
        ny         = self.GetIntConfig("ny")
        compID     = self.GetIntConfig("compid")
        x_points   = np.linspace(coord_1x, coord_2x, nx)
        y_points   = np.linspace(coord_1y, coord_2y, ny)

        nodes = []
        id_cnt = 0

        for y in range(ny):
            tmp = []
            for x in range(nx):
                tmp.append(NekMesh.Node(id_cnt, x_points[x], y_points[y], 0.0))
                id_cnt += 1
            nodes.append(tmp)
        self._create_triangles(nodes, nx, ny, compID)
        # Call the Module functions to create all of the edges, faces and
        # composites.
        self.ProcessVertices()
        self.ProcessEdges()
        self.ProcessFaces()
        self.ProcessElements()
        self.ProcessComposites()

    def _create_triangles(self, nodes, nx, ny, compID):
        ...

# Register our TestInput module with the factory.
NekMesh.Module.Register(
    NekMesh.ModuleType.Input, "StructuredGrid", StructuredGrid)

if __name__ == '__main__':
    # Create a 'pipeline' of the input and output modules.
    mesh = NekMesh.Mesh()

    # First, call our input module's create function from the NekMesh factory.
    NekMesh.InputModule.Create(
        "StructuredGrid", mesh,
        nx = sys.argv[1], ny = sys.argv[2], lx = sys.argv[3],
        ly = sys.argv[4], rx = sys.argv[5], ry = sys.argv[6],
        compid = sys.argv[7], shape = sys.argv[8]).Process()

    # Then ensure there's no negative Jacobians.
    NekMesh.ProcessModule.Create("jac", mesh, list=True).Process()

    # Finally, output the resulting file
    NekMesh.OutputModule.Create(
        "xml", mesh, test=True, outfile=sys.argv[9]).Process()
\end{lstlisting}


\section{NekMesh modules}

\nm 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.
%
\begin{itemize}
  \item \textbf{Input} modules read meshes in a variety of formats;
  \item \textbf{Processing} modules modify meshes to aide in generation processes;
  \item \textbf{Output} modules write meshes in a variety of formats.
\end{itemize}
%
The figure below depicts how these might be coupled together to form a pipeline:
%
\begin{figure}
  \begin{center}
    \begin{tikzpicture}[node distance=90pt]
      \tikzstyle{rect}=[line width=1.6pt,rounded corners=4pt,inner sep=6pt]
      \node[rect,top color=blue!40,bottom color=blue!20,draw=blue!50!black] (A) {Input};
      \node[rect,top color=red!40,bottom color=red!20,draw=red!50!black,right of=A] (B) {Process 1};
      \node[rect,top color=red!40,bottom color=red!20,draw=red!50!black,right of=B] (C) {Process 2};
      \node[rect,top color=green!40,bottom color=green!20,draw=green!50!black,right of=C] (D) {Output};

      \draw[-latex,thick] (A) -- (B);
      \draw[-latex,thick] (B) -- (C);
      \draw[-latex,thick] (C) -- (D);
    \end{tikzpicture}
  \end{center}
  \caption{Illustrative pipeline of the \nm process.}
  \label{fig:util:mc:pipeline}
\end{figure}
%
On the command line, we would define this as:
%
\begin{lstlisting}[style=BashInputStyle]
  NekMesh -m process1 -m process2 input.msh output.xml
\end{lstlisting}
%
Process modules can also have parameters passed to them, that can take
arguments, or not.
%
\begin{lstlisting}[style=BashInputStyle]
  NekMesh -m process1:p1=123:booleanparam input.msh output.xml
\end{lstlisting}
%
To list all available modules use the \inltt{-l} command line argument:
%
\begin{lstlisting}[style=BashInputStyle]
  Available classes:
    Input: dat:
      Reads Tecplot polyhedron ascii format converted from Star CCM (.dat).
  ...
\end{lstlisting}
%
and then to see the options for a particular module, use the \inltt{-p} command
line argument:
%
\begin{lstlisting}[style=BashInputStyle]
  Options for module detect:
          vol: Tag identifying surface to process.
\end{lstlisting}
%
\begin{notebox}
  Module names change when you use the \inltt{-p} option. Input modules should
  be preceded by \inltt{in:}, processing modules by \inltt{proc:} and output
  modules by \inltt{out:}.
\end{notebox}

\subsection{Input modules}

Input and output modules use file extension names to determine the correct
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
    CGNS & \texttt{cgns} & \cmark & Based on the output currently provide by ANSA. Have support for high order elements supported. Requires NEKTAR\_USE\_CGNS option to be activated in cmake which then requires cgns library 4.4 to be  available or compiled \\
    Gmsh & \texttt{msh} & \cmark & Only reads nodes, elements and physical groups (which are mapped to composites). File format versions 2.x and 4.x currently supported.\\
    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.. \\
    Semtex & \texttt{sem} & \cmark & Reads elements and boundary conditions. In order to read high-order information, run \inltt{meshpr session.sem > session.msh} and place in the same directory as the session file.\\
    Star-CCM+ & \texttt{dat} & \xmark & Star outputs plt file which currently needs to be coverted to ascii using Tecplot. Reads mesh only, only support for quads and triangles (2D) and hexes, prisms, tetrahedra (3D).\\
    Star-CCM+ & \texttt{ccm} & \xmark & Reads start ccm format. Reads mesh only, only support for quads and triangles (2D) and hexes, prisms, tetrahedra (3D). Requires NEKTAR\_USE\_CCM option to be activated in cmake and then requires ccmio library to be compiled by user.  \\
    VTK & \texttt{vtk} & \xmark & Experimental support. Only ASCII triangular data is supported. \\
    \bottomrule
  \end{tabularx}
\end{center}

Note that you can override the module used on the command line. For example,
\texttt{Semtex} session files rarely have extensions. So for a session called
\inltt{pipe-3d} we can convert this using the syntax
%
\begin{lstlisting}[style=BashInputStyle]
NekMesh pipe-3d:sem pipe-3d.xml
\end{lstlisting}

The NekMesh input module also has an option to re-process all
composites. By default only composites are reprocessed. However when
extracting a smaller mesh from a larger mesh definition by redefining
the composite of volumetric elements you can force the edges and/or
faces to be reprocessed removing the definition of any edges and/or
faces not required in the smaller mesh. This option is called
\inltt{processall} and has the syntax

\begin{lstlisting}[style=BashInputStyle]
NekMesh input.xml:xml:processall output.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. We note there is one exception to this
convention in three-dimensional meshes where a face composite can
contain both triangular and quadrilateral elements.

Typically a NekMesh call requires both an input and output module to be
called, however, by specifying the output file name or file extension
as \inltt{stdout} no output file will be created. This option is typically
used for mesh statistics processing or inspecting composite values etc.
An example call would be:

\begin{lstlisting}[style=BashInputStyle]
  NekMesh input.xml stdout
\end{lstlisting}

or 

\begin{lstlisting}[style=BashInputStyle]
  NekMesh input.xml name.stdout
\end{lstlisting}


\subsection{Output modules}

The following output formats are supported:

\begin{center}
  \begin{tabularx}{\linewidth}{llcX}
    \toprule
    \textbf{Format} & \textbf{Extension} & \textbf{High-order} & \textbf{Notes}\\
    \midrule
    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. \\
    HDF5 & \texttt{nekg} & \cmark & Most functionality supported. \\
    VTK & \texttt{vtk} & \xmark & Outputs mesh only, supports line segments in 1D, 
    triangles and quadrilaterals in 2D and tetrahedra, hexahedra, prisms and 
    pyramids in 3D. The VTK legacy format and XML format, both compressed and 
    uncompressed, are supported. Requires NEKTAR\_USE\_VTK option to be activated 
    in cmake. \\
    \bottomrule
  \end{tabularx}
\end{center}

Note that for \gmsh, it is highly likely that you will
need to experiment with the source code in order to successfully generate
meshes since robustness is not guaranteed.

The default for \texttt{xml} and \texttt{vtk} 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}. For the uncompressed 
\texttt{xml} the user can execute:
%
\begin{lstlisting}[style=BashInputStyle]
NekMesh Mesh.msh output.xml:xml:uncompress
\end{lstlisting}
%
If the user wishes to obtain the \texttt{vtk} output in the legacy format,
the output module option \inltt{legacy} should be specified by executing:
%
\begin{lstlisting}[style=BashInputStyle]
NekMesh Mesh.xml output.vtk:vtk:legacy
\end{lstlisting}
%
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 \nm.


It is possible to use FieldConvert to extract a smaller region of a mesh
from a larger mesh using the ``-r xmin,xmax,ymin,ymax,zmin'' range
option, e.g.
\begin{lstlisting}[style=BashInputStyle]
FieldConvert -r xmin,xmax,ymin,ymax,zmin,zmax bigMesh.xml smallMesh.xml
\end{lstlisting}
However this will not provide a composite of he faces on the
boundary if they were not part of the original boundary composites of
the original \inltt{bigMesh.xml}. These can however be recovered by an
output module option \inltt{chkbndcomp} in \nm. To do this one should
specify
\begin{lstlisting}[style=BashInputStyle]
NekMesh smallMesh.xml newSmallMesh.xml:xml:chkbndcomp 
\end{lstlisting}
this will then add a composite for any face that is on the boundary
and not part of an existing boundary condition and put it into a
composite with \inltt{id=9999}. This capability is currently only set up for 3D meshes

\subsubsection{HDF5 format}

NekMesh and all solvers within Nektar++ - along with subsequent FieldConvert modules - also support the HDF5 format. This allows for faster loading of geometries and meshes within each solver - and is a significant improvement over the XML format. HDF5 is recommended input format for any larger cases.

Converting from XML to HDF5 is a simple task that only requires the one NekMesh command:
\begin{lstlisting}[style=BashInputStyle]
    NekMesh XMLMesh.xml HDF5Mesh.nekg
\end{lstlisting}

This will create two files \inltt{HDF5Mesh.xml} and \inltt{HDF5Mesh.nekg} which are both needed in the same directory to run the simulation. An additional flag in the session file is required, ensuring it is placed before the expansion list being:

\begin{lstlisting}[style=XmlStyle]
    <GEOMETRY DIM="3" SPACE="3" HDF5FILE="HDF5Mesh.nekg" />
\end{lstlisting}

HDF5 also has the additional advantage of ensuring the mesh and session file are split - which allows for easy ammending of the session file - whilst allowing for use of FieldCovnert modules that require only 1 XML input file - rather than having to concatenate the session and mesh XML files. Solvers and any FieldConvert modules can be run by referencing only the session file after the \inltt{GEOMETRY} tag is included. 

\subsection{Extract surfaces from a mesh}
\label{s:utilities:nekmesh:extract}

Often one wants to visualise surfaces of a 3D mesh, or extract the values of
variables on the surface and visualise them. To support this, \nm can extract
two-dimensional surfaces which can be visualised using \inltt{FieldConvert} in
order to extract the value of a 3D field on a given surface.

As an example, we can extract composite surfaces 2 and 3-5 from a mesh using the
\inltt{extract} module:
%
\begin{lstlisting}[style=BashInputStyle]
NekMesh -m extract:surf=2,3-5 Mesh.xml output.xml
\end{lstlisting}
%
If you also wish to have the boundaries of the extracted surface detected add
the \inltt{detectbnd} option
%
\begin{lstlisting}[style=BashInputStyle]
NekMesh -m extract:surf=2,3-5:detectbnd Mesh.xml output.xml
\end{lstlisting}
which will produce new composites for the extracted boundary.

\subsection{Negative Jacobian detection and histogram visualisation}

To detect elements with negative Jacobian determinant, use the \inltt{jac}
module:
%
\begin{lstlisting}[style=BashInputStyle]
NekMesh -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]
NekMesh -m jac:list Mesh.xml output.xml
\end{lstlisting}
%
To extract the elements whose Jacobian is under a specific value for the purposes of visualisation within the domain,
use the \inltt{extract} boolean parameter. The \inltt{value} should be a number below 1.0 and is 0.0 by default. This means if \inltt{extract} is set by user with no value, the output mesh file will only have invalid mesh elements:
%
\begin{lstlisting}[style=BashInputStyle]
NekMesh -m jac:extract=value Mesh.xml MeshWithNegativeElements.xml
\end{lstlisting}

To show the histogram of Jacobian on the screen, use the \inltt{histo} parameter.

\begin{lstlisting}[style=BashInputStyle]
NekMesh -m jac:histo=value1,value2,value3:quality:detail:histofile=filename Mesh.xml output.xml
\end{lstlisting}

The \inltt{histo} takes 3 input values: The maximum value shown on the histogram, the number of positive bins on the histogram and the number of negative bins. Each value is separated with a comma. By default the values are: 1.0,10,1. This means by default the histogram will have 10 positive bins between (0.0, 1.0) with an interval 0.1, and will have 1 negative bins which shows all the invalid elememnts. An example can be found below in figure \ref{fig:util:mesh:histo} (a).

This command \inltt{histo} also generates a text file called "jacs.txt" which contains the information of each element in the histogram. The name of the file can be changed by using the command \inltt{histofile}. The information in this text file are: Element ID, Jacobian value, element type, if this element is a boundary element, the boundary edge/face ID, vertex coordinates and neighbor element ID will also be included. Note that only elements whose Jacobian is smaller than "value1" will be counted, and they will be outputted in the histogram and text file. If command \inltt{detail} is added, the composite ID and name of each element will also be outputted into the text file. An example is shown in figure \ref{fig:util:mesh:histo} (c). 

The command \inltt{quality} is an option which show the percentage of each bin of the histogram. The results of \inltt{quality} always match the histogram. An example is shown in figure \ref{fig:util:mesh:histo} (b). 

An example with the following command is shown in figure \ref {fig:util:mesh:histo}. The results of \inltt{quality} and text file always correspond to the histogram.

\begin{lstlisting}[style=BashInputStyle]
NekMesh -m jac:histo=0.9,3,2:quality:detail Mesh.xml output.xml
\end{lstlisting}

\begin{figure}[!htbp]
  \begin{center}
    \includegraphics[width = 0.47 \textwidth]{img/ProcessJac_histo.png}
    \includegraphics[width = 0.47 \textwidth]{img/ProcessJac_quality.png}
    \includegraphics[width = 1.0 \textwidth]{img/ProcessJac_output.png}
    \caption{(a) Results of option "histo": histogram shown on the screen, (b) Results of option "quality", (c) Part of the text file "jacs.txt"}
    \label{fig:util:mesh:histo}
  \end{center}
\end{figure}

To turn off curvature associated with negative jacobians one can try to use the
\inltt{linearise} module:
\begin{lstlisting}[style=BashInputStyle]
NekMesh -m linerise:invalid  Mesh.xml output.xml
\end{lstlisting}
This option will remove all high order curvature on all element types
with singular jacobians.

\subsection{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
high-order element. In \nm we use \emph{spherigons}, a kind of patch used in the
computer graphics community used for efficiently smoothing polygon surfaces.

Spherigons work through the use of surface normals, where in this sense
`surface' refers to the underlying geometry. If we have either the exact or
approximate surface normal at each given vertex, spherigon patches approximate
the edges connecting two vertices by arcs of a circle. In \nm we can either
approximate the surface normals from the linear elements which connect to each
vertex (this is done by default), or supply a file which gives the surface
normals.

To apply spherigon patches on two connected surfaces 11 and 12 use the following
command:
%
\begin{lstlisting}[style=BashInputStyle]
NekMesh -m spherigon:surf=11,12 \
    MeshWithStraighEdges.xml MeshWithSpherigons.xml
\end{lstlisting}
%
If the two surfaces "11" and "12" are not connected, or connect at a sharp edge
which is $C^0$ continuous but not $C^1$ smooth, use two separate instances of
the spherigon module.
%
\begin{lstlisting}[style=BashInputStyle]
NekMesh -m spherigon:surf=11 -m spherigon:surf=12 \
    MeshWithStraighEdges.xml MeshWithSpherigons.xml
\end{lstlisting}
%
This is to avoid the approximated surface normals being incorrect at the edge.

If you have a high-resolution mesh of the surfaces 11 and 12 in \inltt{ply}
format it can be used to improve the normal definition of the spherigons. Run:
\begin{lstlisting}[style=BashInputStyle]
NekMesh -m spherigon:surf=11,12:usenormalfile=Surf_11-12_Mesh.ply \
    MeshWithStraighEdges.xml MeshWithSpherigons.xml
\end{lstlisting}

This can be useful, for example, when meshing the Leading edge of an
airfoil. Starting from a linear mesh (left figure) the spherigon patches curve
the surface elements producing leading edge closer to the underlying geometry:

\begin{figure}[!htbp]
  \begin{center}
    \includegraphics[width = 0.47 \textwidth]{img/noSphnoBL.jpg}
    \includegraphics[width = 0.47 \textwidth]{img/SphnoBL.jpg}
    \caption{(a) Leading edge without spherigons, (b) Leading edge with
      spherigons}
  \end{center}
\end{figure}

\subsection{Periodic boundary condition alignment}

When using periodic boundary conditions, the order of the elements within the
boundary composite determines which element edges are periodic with the
corresponding boundary composite.

To facilitate this alignment, \nm has a periodic alignment module
which attempts to identify pairs of mutually periodic edges. Given two
surfaces \inltt{surf1} and \inltt{surf2}, which for example correspond
to the physical surface IDs specified in \gmsh, and an axis which
defines the periodicity direction, the following command attempts to
reorder the composites:
%
\begin{lstlisting}[style=BashInputStyle]
NekMesh -m peralign:surf1=11:surf2=12:dir=y \
    -m peralign:surf1=13:surf2=14:dir=z Mesh.xml Mesh_aligned.xml
\end{lstlisting}
%
Here the surfaces with IDs 11 and 12 will be aligned normal to the $y$-axis and
the surfaces 13 and 14 will be aligned normal to the $z$-axis.

Note that this command cannot perform magic -- it assumes that any given edge or
face lying on the surface is periodic with another face on the opposing surface,
that there are the same number of elements on both surfaces, and the
corresponding edge or face is the same size and shape but translated along the
appropriate axis.


When using periodic boundary conditions that are rotationally aligned the
following rotational options should be applied:
%
\begin{lstlisting}[style=BashInputStyle]
NekMesh -m peralign:surf1=11:surf2=12:dir=x:rot=PI/6 \
              Mesh.xml Mesh_aligned.xml
\end{lstlisting}
%
where \inltt{rot} specifies the rotation angle in radians from
surf1 to surf2 about the axis specified by \inltt{dir} (i.e. the ``x'' axis in
this example).

The rotation/translation is assumed to be exact within a relative tolerance.
An optional factor, which is used to scale the tolerance, \inltt{tolfact}
can also be specified. The default tolerance factor is 4, and it needs to be
$tolfact \ge 1$. For example:
%
\begin{lstlisting}[style=BashInputStyle]
NekMesh -m peralign:surf1=11:surf2=12:dir=x:rot=PI/6:tolfact=100 \
              Mesh.xml Mesh_aligned.xml
\end{lstlisting}
%

In 3D, where prismatic or tetrahedral elements are connected to one or both of
the surfaces, additional logic is needed to guarantee connectivity in the XML
file. In this case we append the \inltt{orient} parameter:
%
\begin{lstlisting}[style=BashInputStyle]
NekMesh -m peralign:surf1=11:surf2=12:dir=y:orient input.dat output.xml
\end{lstlisting}

\begin{notebox}
  One of the present shortcomings of \inltt{orient} is that it throws away all
  high-order information and works only on the linear element. This can be
  gotten around if you are just doing e.g. spherigon patches by running this
  \inltt{peralign} module before the \inltt{spherigon} module.
\end{notebox}

\subsection{Boundary layer splitting}

Often it is the case that one can generate a coarse boundary layer grid of a
mesh. \nm has a method for splitting prismatic and hexahedral elements into
finer elements based on the work presented in~\cite{MoHaPeSh14}
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
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
elements to construct a sequence of increasingly refined elements, as depicted
in figure~\ref{fig:util:mc:split}.

\begin{figure}
  \begin{center}
    \begin{tikzpicture}
      \node[above right] at (0,0){%
        \includegraphics{img/stdprism_split}};
      \node[above right] at (6,0){%
        \includegraphics[width=5cm]{img/prism_split}};
      \draw[-latex,thick,bend left] (4.5,3.75) to[bend left]
      node[midway,above] {$\chi^e(\mathbf{\xi})$} (8.5,3.75);
    \end{tikzpicture}
  \end{center}
  \caption{Splitting $\Omega_{\text{st}}$ and applying the mapping $\chi^e$ to
    obtain a high-order layer of prisms from the macro-element.}
  \label{fig:util:mc:split}
\end{figure}

To split a prism boundary layer on surface 11 into 3 layers with a growth rate
of 2 and 7 integration points per element use the following command:
\begin{lstlisting}[style=BashInputStyle]
  NekMesh -m bl:surf=11:layers=3:r=2:nq=7 MeshWithOnePrismLayer.xml \
        MeshWith3PrismsLayers.xml
\end{lstlisting}
%
\begin{figure}[!htbp]
  \begin{center}
    \includegraphics[width = 0.47 \textwidth]{img/SphnoBL.jpg}
    \includegraphics[width = 0.47 \textwidth]{img/SphBL.jpg}
    \caption{(a) LE with Spherigons but only one prism layer for resolving the
      boundary layer, (b) LE with Spherigons with 3 growing layers of prisms for
      better resolving the boundary layer.}
  \end{center}
\end{figure}

\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=sin(x)}. In this case the
  function should be sufficiently smooth to prevent the elements
  self-intersecting.
\end{notebox}

\subsection{High-order cylinder generation}

Generating accurate high-order curved geometries in \gmsh is quite challenging.
This module processes an existing linear cylindrical mesh, with axis aligned
with the $z$-coordinate axis, to generate accurate high-order curvature
information along the edges.

\begin{lstlisting}[style=BashInputStyle]
NekMesh -m cyl:surf=2:r=1.0:N=5 LinearCylinder.xml HighOrderCylinder.xml
\end{lstlisting}

The module parameters are:

\begin{itemize}
  \item \inlsh{surf}: Surface on which to apply curvature. This should be the
  outer surface of the cylinder.
  \item \inlsh{r}: Radius of the cylinder.
  \item \inlsh{N}: Number of high-order points along each element edge.
\end{itemize}

\begin{notebox}
  The module could also be used to apply curvature along the interior of a
  hollow cylinder. However, there are no checks to ensure the resulting elements
  are not self-intersecting.
\end{notebox}

\subsection{Linearisation}

The ability to remove all the high-order information in a mesh can be useful at
times when mesh generation is tricky or impossible in the presence of curvature

To do this in NekMesh use the command:

\begin{lstlisting}[style=BashInputStyle]
  NekMesh -m linearise:all high-order-mesh.xml linear-mesh.xml
\end{lstlisting}

The output will contain only the linear mesh information, all curved information
is removed. Alternatively

\begin{lstlisting}[style=BashInputStyle]
  NekMesh -m linearise:invalid high-order-mesh.xml linear-mesh.xml
\end{lstlisting}

attempts to remove curvature from elements only where necessary. This is a
simple algorithm that removes curvature from invalid elements and repeats until
all elements are valid. Either \inlsh{all} or \inlsh{invalid} must be specified.

\begin{itemize}
  \item \inlsh{all}: remove curvature from all elements.
  \item \inlsh{invalid}: remove curvature from invalid elements.
  \item \inlsh{prismonly}: consider only prisms when removing curvature. This is
  useful in the presence of a prismatic boundary layer.
\end{itemize}

\subsection{Extracting interface between tetrahedra and prismatic elements}

When the mesh is three-dimensional and comprised of a prismatic boundary layer
with tetrahedra in the interior of the domain, this module extracts the
prismatic elements only, and constructs a boundary region for the interface
between the tetrahedra and prisms. This is useful in, for example, the study of
aortic flows, where the prismatic boundary layer can be extracted and refined to
study unsteady advection-diffusion problems on a more refined grid inside the
boundary layer.

To use this module you therefore use the command:

\begin{lstlisting}[style=BashInputStyle]
  NekMesh -m extracttetprisminterface input.xml output.xml
\end{lstlisting}

There are no configuration options for this module, as it is highly specific to
a certain class of meshes.

\subsection{Boundary identification}

Some mesh formats lack the ability to identify boundaries of the domain they
discretise. \nm has a rudimentary boundary identification routine for conformal
meshes, which will create a composite of edges (2D) or faces (3D) which are
connected to precisely one element. This can be done using the \inltt{detect}
module:

\begin{lstlisting}[style=BashInputStyle]
  NekMesh -m detect volume.xml volumeWithBoundaryComposite.xml
\end{lstlisting}

\subsection{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]
  NekMesh -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}

\subsection{Link Checking}

It is quite possible that a mesh contains some sort of hanging entity or
element connectivity error. The check link module is a fast check that, a)
elements are correctly connected and b) the boundary entities (composites)
match the interior domain:

\begin{lstlisting}[style=BashInputStyle]
  NekMesh -m linkcheck mesh.xml mesh2.xml
\end{lstlisting}

This module should be added to the module chain if the user suspects there may
be a mesh issue. The module will print a warning if there is a connectivity
error.

\subsection{2D mesh extrusion}

This module allows a 2D mesh, quads, triangles or both, to be extruded in the
$z$ direction to make a simple 3D mesh made of prisms and hexahedra. It is also
capable of extruding the high-order curvature within the 2D mesh. The module
requires two parameters:

\begin{lstlisting}[style=BashInputStyle]
  NekMesh -m extrude:layers=n:length=l 2D.xml 3D.xml
\end{lstlisting}

length which determines how long the $z$ extrusion will be and layers, the
number of elements in the $z$ direction.

\subsection{2D mesh revolution}

This module allows a 2D mesh, quads, triangles or both, to be revolved around
the $y$ axis to make a simple 3D mesh made of prisms and hexahedra. It is also
capable of revolving the high-order curvature within the 2D mesh. The module
requires two parameters:

\begin{lstlisting}[style=BashInputStyle]
  NekMesh -m revolve:layers=n:angle=a 2D.xml 3D.xml
\end{lstlisting}
angle which determines the angle of revolution ($\theta$) around the $y$ axis and 
layers, the number of elements in the $\theta$ direction.

The angle can also be set to "full" for a full revolution of $2\pi$, although
this is the default setting:

\begin{lstlisting}[style=BashInputStyle]
NekMesh -m revolve:layers=n:angle=full 2D.xml 3D.xml
\end{lstlisting}

The angle must be between $0$ and $2\pi$.
The 2D mesh must not cross or touch the $y$ axis i.e. $x>0$ for all points.
By default the mesh is revolved through the full $2\pi$.

\subsection{Variational Optimisation}
\label{sec:varopti}
This module can correct invalid and improve the quality of elements in
high-order meshes by applying curvilinear deformation to the interiors of
domains. It achieves this by solving a solid mechanics system which, using
variational calculus has been cast is a non-linear energy optimsation problem.
It is basis of the work in \cite{TuPeMo16}.

It works by considering the boundary (curved) mesh entities to be fixed and
moving the interior nodes to a lower  energy configuration. This new
configuration in most scenarios is a higher quality mesh. The energy is
evaluated depending on which functional is chosen. We find hyperleastic to be
the most reliable but it can also model the mesh and a linearelastic solid as
well as functionals based on the Winslow equation and the distortion method
proposed by Roca et al. \cite{GaRoPeSa15}.

There are a large number of options which can be viewed using the help function
but the basic usage is:
\begin{lstlisting}[style=BashInputStyle]
  NekMesh -m varopti:type inital.xml optimised.xml
\end{lstlisting}

where type can be \texttt{hyperelastic}, \texttt{linearelastic},
\texttt{winslow} or \texttt{roca}.

\subsection{\emph{r}-adaptation: interpolation file}
This module can deform an existing mesh by using the variational optimiser
presented above. A file must be provided that contains a list of points
and a scaling value for each of them. This scaling factor is then used to
target an element size based on the initial size of the element. Scaling
values are interpolated throughout the domain based on the interpolation
method of the main library. The file should look like
\begin{lstlisting}
0 0 0 2.0
0 1 0 2.0
1 0 0 0.5
1 1 0 0.5
\end{lstlisting}
where the first three columns are \texttt{x}, \texttt{y}, \texttt{z} and
the last column is the scaling factor.

The call is identical to the variational optimisation module above:
\begin{lstlisting}[style=BashInputStyle]
  NekMesh -m varopti:type:scalingfile=file.txt:subiter=x inital.xml adapted.xml
\end{lstlisting}
where \texttt{subiter} is an additional parameter to the variational
optimiser that defines the frequency at which individual elements
update their target scaling based on their latest location in the domain.
\texttt{subiter} should be a scalar and is the number of steps between updates.
It is often recommended to run \emph{r}-adaptation on a linear mesh for
stability and performance reasons. Note also that the mesh must have
CAD information in order for nodes to slide on curves and surfaces.

\subsection{\emph{r}-adaptation: CAD curves}
Another option for \emph{r}-adaptation can be defined where instead of providing
an interpolation file, the user can provide a list of CAD curves towards which
the scaling takes place. With this option, the \emph{r}-adaption uses a fixed
scaling factor for all elements which are on or within a specified distance from
the curve.

The call is similar to the \emph{r}-adaption with interpolation file, with the difference of replacing the file with a list of CAD curve IDs, the scaling factor and an optional distance from the curve:
\begin{lstlisting}[style=BashInputStyle]
  NekMesh -m varopti:type:radaptcurves=CurveIDs:radaptscale=scale:radaptrad=rad:subiter=x inital.xml adapted.xml
\end{lstlisting}
where \texttt{radaptcurves}, \texttt{radaptscale} and \texttt{radaptscale} are
additional containing the list of CAD curve IDs (integers separated by commas,
or with a dash for ranges, i.e. 1-3,6 would select curves with ID 1, 2, 3, and
6), the scaling factor as a float and the distance from the CAD curves as a
float, respectively. As with the \emph{r}-adaption with interpolation file,
\texttt{subiter} is to define the frequency at which individual elements update
their target scaling based on their current proximity to the selected CAD
curves.

\subsection{Mesh Curving}
This module can take any linear mesh or high-order mesh, given that it is a close
representation of the CAD and project the boundary of the mesh onto the CAD.
This will reconstruct the connectivity between the CAD and the surface mesh.
Then it will curve and optimize the surface of the mesh. The method has a number
of fail-safes ensuring that even bad CAD or poor linear meshes should be able 
to be curved. If the method encounters an issue, such as the linear mesh being a
large distance from the CAD, it will simply leave that element straight sided. A
well made CAD model and accurate linear mesh should be curved with little issue.
This module has been tested with hybrid prism-tet-pyramid, fully hexahedral and
tetrahedral meshes. As a best practice, it is recommended to use a single thick
O-type prismatic or hexahedral layer that is with sufficient thickness to account
for the mesh curving. 

The module needs to be informed of the CAD file to curve the mesh on the CAD and the
order at which to curve the surface:
\begin{lstlisting}[style=BashInputStyle]
  NekMesh -m projectcad:file=cadfile.step:order=x:cLength=y: inital.xml optimised.xml
\end{lstlisting}

The new version of the module reuses all state-of-the-art mesh curving modules with additional functionalities. 

\begin{itemize}
  \item \inlsh{surfopti}: True by default.Activates/Deactivates mesh curving and surface optimization. 
  \item \inlsh{ho}: Needs to be passed if the input mesh is already high-order.
  \item \inlsh{tolv1}: (optional) Tolerance for max distance of a node to CADSurface.
  \item \inlsh{tolv2}: (optional) Tolerance for max distance  (< tolv1) of a node to CADSurface.
  \item \inlsh{cLength}: (optional) Scales the tolerances by the factor, expands the
  bounding boxes and speeds up the bounding box generation. This is recommended
  to be tunned by the users need if the default values do not produce a good result.
  \item \inlsh{extract}: (optional) Exports the CAD-mesh connectivity to CAD.txt.
  \item \inlsh{cadtxt}:  (optional) Imports the CAD-mesh connectivity from CAD.txt directly.
\end{itemize}

Note: Please ensure that the STEP file is using the same dimension as the linear mesh.
A common error is to have a STP file in MiliMeters, whereas the mesh is in Meters or 
vice-versa. A way to check is to see the bounding box of printed in the beginning of
the module is in the correct order. 

\subsection{Mesh combine}
This module can combine two mesh files into one output file.
This will read all information from two mesh files and generate the combine meish 
file with all these information. The process will duplicate composite IDs 
from mesh 1 detected in mesh 2 and will remap into the output file.

\begin{lstlisting}[style=BashInputStyle]
  NekMesh -m combine:file=mesh1.xml mesh2.xml combine.xml
\end{lstlisting}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Mesh generation}

In addition to the functionality described previously, \nm is capable of
generating high-order meshes directly from a CAD definition. By default this
functionality is not activated, a user wishing to utilise the mesh generation
capability of \nm must compile \nekpp with the \texttt{NEKTAR\_USE\_MESHGEN}
option on. As well as compiling the relevant routines into \nm it will also
download a number of other packages which are required.

The most critical dependancy of the mesh generation routines is
\texttt{OpenCascade} which powers the CAD engine. \nm is capable of finding and
using existing installations of \texttt{OpenCascade 6.8} or \texttt{OCE 0.17}.
If either are not present on the installation machine \nm will install
\texttt{OCE 0.17} from source. This is a very big installation and will take
some time so it is advised that the user ensures \texttt{OpenCascade} is
availble on the machine.

As with all tasks within \nm the mesh generation capability exists as its
own separate module which is of type Input. Due to the vast amount of code
associated with the generation of high-order meshes and the comparatively small
nature of modules in the \nm program a new library has been created for \nekpp
called {\em NekMeshUtils}, which contains all the core routines and classes for
the \nm mesh format as well as a series of classes for the generation of meshes.
This library also contains the CAD API for \nekpp which is used to generate the
meshes.

\subsection{Methodology}

This section outlines the approach taken by \nm to generate high-order meshes.
%
To simplify the sometimes very complicated high-order mesh generation processes
in other programs, \nm executes all the stages required to produce a high-order
mesh in one single pipeline which once started requires no interaction from the
user. In broad terms these stages are:
\begin{itemize}
	\item Specification of the element sizes in the mesh,
	\item Coarse linear mesh generation of the domain,
	\item Generation of optimised high-order surface on the geometric boundary,
\end{itemize}
%
and are outlined in more detail in the following sections.

\subsubsection{CAD Interaction}

At the core of all the ideas in the \nm generator is that the final mesh is a
high quality representation of the underlying geometry. As such all of the
entities in the mesh must know where they are located with respect to the CAD
and the system to be able to query any geometric information at any point in the
domain easily and with accuracy. To handle this \nm has been interfaced with the
third-party suite of CAD libraries called OpenCascade. In its normal state
OpenCascade is a very large collection of libraries with tens of thousands of
functions which are simply not needed for our purposes, because of this its
installation is a very arduous and long process. Combine this with the fact that
there are dozens of versions and types of OpenCascade, such as OpenCascade
Community Edition, it is simply impossible for \nm to use already existing
OpenCascade installations on a given machine. To solve these issues, when
installing \nekpp with the mesh generator it will download pre-compiled binaries
for the relevant OS and link against those, any previously installed versions of
OpenCascade will not be searched for and therefore ignored.
%
To reduce the massively complex libraries in OpenCascade down to a manageable
set of functions to be used in \nm a set of interface classes have been created
which act as buffer between it and \nekpp. These CAD classes mean that
development of mesh generation routines is significantly easier and in the
future \nekpp developers will be able to utilise CAD information in all aspects
of the framework without having to learn OpenCascade. Another advantage with
this approach is that adding support for other CAD engines, as well as
OpenCascade, in the future should be relativity simple and will not require the
rewriting of any of the \nm code.

\subsubsection{Automatic specification of the mesh}

One of the key challenges of generating a high-order mesh is the creation of a
suitable coarse linear mesh. It is
quite difficult for a user to define a full set spacings over a whole domain
which will produce a good quality especially when aiming for coarseness. This
is tackled in \nm with a system for automatically defining a
set of smooth and coarse mesh spacings throughout the whole domain. This is
achieved using an octree description of the domain. The domain is recursively
subdivided into octants which each describe a small portion of the domain. The
level to which the domain subdivides is based on the curvature of the geometric
boundary. Higher curvature regions will subdivide to a finer level allowing for
increased control on the mesh specification and smoothness. The geometric
curvature is then related to a mesh sizing parameter and propagated throughout
the domain ensuring a smooth mesh. For those unfamiliar with octrees, it is
best to think of it as a non-conforming hexahedral mesh

\subsubsection{Linear Mesh Generation}

The first challenge mentioned in the previous section is addressed with the
\nm approach to linear mesh generation. Primarily because of the
difficulties in interfacing existing linear mesh generators for high-order
applications the decision was made to include a bespoke linear mesh generator
within the program. Compared with the mesh generators included in commercial
packages this linear mesh generator takes the quite unconventional and more
historic approach in building the mesh in a bottom up fashion from 0D to 3D.
Using this approach means it is possible to guarantee a level of boundary
conformity which direct to 3D approaches cannot at the desired level of
coarseness. In this approach, first mesh nodes are placed on the vertices of
the CAD model (0D), then the curves in the CAD are meshed in 1D using the
vertex nodes as boundaries, then the surfaces are meshed in their 2D parameter
plane using the curve meshes as boundaries and finally the 3D volume is meshed
using the surface mesh as the boundary to complete the linear mesh. In \nm, to
achieve greater robustness, the 2D mesh generation library Triangle is used and
the TetGen library for the 3D. Both of which are highly developed Delaunay based
mesh generators. As with all additional libraies in \nekpp these are
automatically downloaded and installed if needed.

\subsubsection{High-order Surface Generation}

Addition of the high-order nodes to and the curving of the mesh is very open
problem, no high-order mesh generator has solved this and while the methods
used in \nm are not 100\% full-proof, the system currently in place can create
good quality high-order curved meshes with a reasonable robustness. This area
will receive the greatest level of development in the future. The most critical
part of defining the high-order mesh is the addition of high-order nodes on the
geometric surface. The mesh generator must achieve the greatest level of
geometric accuracy as it can otherwise it will greatly affect the final flow
solutions. If the linear surface triangulation is taken to be fixed during this
process, the problem can be addressed in a element by element fashion. If the
high-order nodes are placed by simply using an affine mapping to the CAD surface
and back the resulting high-order triangle will inherit the same distortions as
the CAD surface. To solve this \nm uses a system node location optimisation in
the parameter plane of the CAD surface to ensure the high-order triangles have
as little distortion as possible while remaining exactly on the geometric
surface. To do this the system models the high-order edges and triangles as a
network of springs with an associated spring energy which is minimised using a
multidimensional Newton type optimisation procedure with a Gauss-Seidel matrix
solver.

\subsubsection{Mesh Correction}

Due to the fact that, for the time being, no consideration is given to the
curving of mesh interior entities explicitly in the mesh generation process, the
curving the geometric surface can produce meshes with invalid elements,
especially in the case of Euler type (Tetrahedra only) meshes. Three strategies
exist within \nekpp to correct these elements. Firstly removing the curvature,
by removing the curvature of invalid elements they become valid. However this
has the massive downside of compromising the geometric accuracy of the mesh but
is quick and effective, this can be enacted using the command:
\begin{lstlisting}[style=BashInputStyle]
    NekMesh -m linearise:invalid invalidMesh.xml validMesh.xml
\end{lstlisting}
%
An alternative to this is to use the linear elastic solver within \nekpp to
deform the mesh interior entities. Its use is very computationally expensive,
as with all PDE solvers, and is also not particularly robust. It can be used
with the set of commands outlined in the FieldConvert deform and displacement
modules and the section on the Linear Elastic Solver.

The final and possibly most useful approach is to use the Variational
Optimsation module to curve the interior of the domain. This is explained in
\ref{sec:varopti}.

\subsection{Mesh generation manual}

The mesh generation is executed with the command:
\begin{lstlisting}[style=BashInputStyle]
    NekMesh session.mcf mesh.xml
\end{lstlisting}
where session.mcf is a mesh configuration file which contains all the options
and parameters needed for mesh generation. Below is an example of a simple
example which generates a 2D NACA wing.
\begin{lstlisting}[style=XmlStyle]
    <NEKTAR>
        <MESHING>

            <INFORMATION>
                <I PROPERTY="CADFile"   VALUE="6412"  />
                <I PROPERTY="MeshType"  VALUE="2D"    />
            </INFORMATION>

            <PARAMETERS>
                <P PARAM="MinDelta"     VALUE="0.01"  />
                <P PARAM="MaxDelta"     VALUE="1.0"   />
                <P PARAM="EPS"          VALUE="0.1"   />
                <P PARAM="Order"        VALUE="4"     />

                <!-- 2D Domain !-->
                <P PARAM="Xmin"         VALUE="-1.0"  />
                <P PARAM="Ymin"         VALUE="-2.0"  />
                <P PARAM="Xmax"         VALUE="3.0"   />
                <P PARAM="Ymax"         VALUE="2.0"   />
                <P PARAM="AOA"          VALUE="15.0"  />
            </PARAMETERS>

        </MESHING>
    </NEKTAR>
\end{lstlisting}
In all cases the mesh generator needs two pieces of information and four
parameters. It firstly needs to know the CAD file with which to work. In the
example above this is listed as a 4 digit number, this is because the mesh
generator is equiped with a NACA wing generator. In all other cases this
parameter would be the name of a CAD file (in either STEP or GEO format).
Secondly, what type of mesh to make, the options
are \inltt{EULER} and \inltt{BndLayer} for 3D meshes and \inltt{2D} and \inltt{2DBndLayer}
for 2D meshes. In the case of \inltt{EULER} the mesh will be made with only
tetrahedra. For \inltt{BndLayer} the mesh generator will attempt to insert a single
macro prism layer onto the geometry surface. This option requires additional
parameters. This is similar for the 2D scenarios. The automatic mesh
specification system requires three parameters to build the specification of a
smooth, curvature refined mesh. Firstly \inltt{MinDelta} which is the size of
the smallest element to be found in the final mesh. Secondly \inltt{MaxDelta}
which is the maximum size of an element in the mesh and lastly \inltt{EPS} which
is a sensitivity to curvature parameter with a range $1 \geq \varepsilon > 0$
which heuristically controls the size of the elements for a given degree of
curvature on the geometric surface. \inltt{Order} is the polynomial order of the
mesh to be generated.
%
When generating a boundary layer mesh a few additional parameters must be given.
An example is shown.
\begin{lstlisting}[style=XmlStyle]
    <NEKTAR>
        <MESHING>

            <INFORMATION>
                <I PROPERTY="CADFile"           VALUE="6412"        />
                <I PROPERTY="MeshType"          VALUE="2DBndLayer"  />
            </INFORMATION>

            <PARAMETERS>
                <P PARAM="MinDelta"             VALUE="0.01"        />
                <P PARAM="MaxDelta"             VALUE="1.0"         />
                <P PARAM="EPS"                  VALUE="0.1"         />
                <P PARAM="Order"                VALUE="4"           />

                <!-- Boundary layer !-->
                <P PARAM="BndLayerSurfaces"     VALUE="5,6"         />
                <P PARAM="BndLayerThickness"    VALUE="0.03"        />
                <P PARAM="BndLayerLayers"       VALUE="4"           />
                <P PARAM="BndLayerProgression"  VALUE="2.0"         />

                <!-- 2D Domain !-->
                <P PARAM="Xmin"                 VALUE="-1.0"        />
                <P PARAM="Ymin"                 VALUE="-2.0"        />
                <P PARAM="Xmax"                 VALUE="3.0"         />
                <P PARAM="Ymax"                 VALUE="2.0"         />
                <P PARAM="AOA"                  VALUE="15.0"        />
            </PARAMETERS>

        </MESHING>
    </NEKTAR>
\end{lstlisting}
A list of the CAD surfaces which will have a prism generated on is described by
\inltt{BndLayerSurfaces} and the thickness of the boundary to aim for is \inltt{BndLayerThickness}.
%
The mesh generator has been created with a range of error messages to aid in
debugging. If you encounter an error and the mesh generator fails, run \nm with
the verbose \inltt{-v} flag and send the stdout with the .mcf and .step files to
\inltt{d.moxey@exeter.ac.uk}. Without the feedback this functionality cannot
improve.

\subsubsection{Handling 3D geometries with voids}

Although NekMesh supports the definition of 3D geometries that contain voids --
for example, a sphere contained within a cube -- at present it does require the
definition of a point per-void that lies strictly on the interior of the
void. This is so that tetrahedra on the interior of the void can be removed
before the mesh is generated.
%
For example, if one defines a geometry where two spheres of radius 1, centred at
$(0,0,0)$ and $(2,0,0)$, were contained within a larger domain, the void points
can be specified through the \inltt{VOIDPOINTS} tag in the MCF as follows:

\begin{lstlisting}[style=XmlStyle]
    <NEKTAR>
        <MESHING>
            <!-- other parameters here... -->
            <VOIDPOINTS>
                <V> 0 0 0 </V>
                <V> 2 0 0 </V>
            </VOIDPOINTS>
       </MESHING>
    </NEKTAR>
\end{lstlisting}

\subsubsection{GEO format}

Recent developments have been made to facilitate the generation of meshes from
simple 2D and 3D geometries. The GEO file format, used by Gmsh, is a popular
option that allows the user to script geometrical and meshing operations without
the need of a GUI. A simplified reader has been implemented in NekMesh for 2D
and 3D geometries.  Although very basic this reader may be extended in the
future to cover a wider range of geometrical features.

For a full description of the GEO format the user should refer to Gmsh's
documentation. The following commands are currently supported:
\begin{itemize}
  \item \inltt{//} and \inltt{/* */} (i.e. comments)
  \item \inltt{Point}
  \item \inltt{Line}
  \item \inltt{Spline} (through points)
  \item \inltt{BSpline} (i.e. a B\'{e}zier curve)
  \item \inltt{Ellipse} (arc): as defined in Gmsh's OpenCASCADE kernel, the first
    point defines the start of the arc, the second point the centre and the fourth
    point the end. The third point is not used. The start point along with the centre
    point form the major axis and the minor axis is then computed so that the end
    point falls onto the arc. The major axis must always be greater or equal to the
    minor axis.
  \item \inltt{Circle} (arc): the circle is a special case of the ellipse where the
    third point is skipped. The distances between the start and end points and the
    centre point must be equal or an error will be thrown.
  \item \inltt{Line Loop}
  \item \inltt{Plane Surface}
  \item \inltt{Ruled Surface} or, in newer versions of Gmsh, \inltt{Surface}
  \item \inltt{Surface Loop}
  \item \inltt{Volume}
\end{itemize}

At the present time, NekMesh does not support the full scripting capabilities of
the GEO format, but the evaluation of simple variables is supported. The used
GEO files should be a straightforward succession of entity creations (see list
above). This should however allow for the creation of quite a wide range of 2D
and 3D geometries by transformation of arbitrary curves into generic splines and
arcs.

\subsubsection{CCM format}
For very complex 3D geometries with manifold bodies, \nekpp can take advantage of the meshing tool embedded in the Star-CCM package. \nm projects linear meshes generated in Star-CCM on the CAD definition to generate high-order meshes. \nekpp must be built \texttt{NEKTAR\_USE\_CCM} in conjunction with the \texttt{ccmio} package to use this functionality. Projection is performed by specifying the option \inltt{projectcad}.

The supported element types for the linear mesh are :
\begin{itemize}
  \item \inltt{Tetrahedron}
  \item \inltt{Prism}
  \item \inltt{Pyramid}
\end{itemize}

Therefore it is recommended to use the \texttt{Surface Remesher}, the \texttt{Tetrahedral Mesher} and the \texttt{Prism Layer Mesher} in Star-CCM. The \texttt{CAD Projection} and \texttt{Generate Only Standard Prismatic Cells} options should be activated. 
The prismatic layer for the linear mesh in Star-CCM should not be splitted e.g. there should be only one layer with the total thickness of the boundary layer (see figure \ref{fig:util:mc:Car_mesh}. Splitting will then be perfomed directly by \nm on surfaces specified by the argument \inltt{surf} and the number of layers (\inltt{layers} argument) and growth rate (\inltt{r} argument) are specified during the mesh conversion process with the \inltt{bl} option. Surfaces correspond to the regions created in Star-CCM and the reference index is displayed when running \nm with verbosity (\inltt{-v} option).

\begin{lstlisting}[style=BashInputStyle]
  NekMesh -m projectcad:file=cadFile.STEP:order=meshOrder -m bl:surf=surf1,surf2,surf3:r=growthRate:layers=nLayers linearMesh.ccm outputMesh.xml
\end{lstlisting}

The command \inltt{linearise} can be added to correct bad elements.

\begin{figure}
  \begin{center}
    \includegraphics[width=1.0\textwidth]{img/Car_mesh.png}
  \end{center}
  \caption{Linear mesh generated in Star-CCM (visualised with ParaView) showing the one-layer prismatic cells layer}
  \label{fig:util:mc:Car_mesh}
\end{figure}
%%% Local Variables:
%%% mode: latex
%%% TeX-master: "../user-guide"
%%% End:
