Commit 2e242353 authored by Chris Cantwell's avatar Chris Cantwell

Merge branch 'master' into fix/for-visit

parents a15ad44d 508bc973
......@@ -25,7 +25,13 @@ ENABLE_TESTING()
# Use the CMake 2.4 policy for CMake 2.6
IF(COMMAND cmake_policy)
cmake_policy(SET CMP0003 OLD)
# Silence warnings when using generator expressions in
ENDIF(COMMAND cmake_policy)
# Find the modules included with Nektar
......@@ -2,7 +2,7 @@ FIND_LIBRARY(WIN32_BLAS NAMES libblas PATHS ${TPSRC})
SET(WIN_ZIP_MD5_VALUE "b5ad4f7335ca964bbdafab129e44a046")
SET(WIN_ZIP_SHA1_VALUE "adb331fa195db264e95e46b902887f12971dbf48")
MESSAGE(STATUS "Found Win32 Lapack: ${WIN32_LAPACK}")
\ No newline at end of file
MESSAGE(STATUS "Found Win32 Lapack: ${WIN32_LAPACK}")
......@@ -62,6 +62,7 @@ IF( NEKTAR_USE_MPI )
......@@ -20,13 +20,14 @@ EXTERNALPROJECT_ADD(
TMP_DIR ${TPBUILD}/modmetis-5.1.0-tmp
......@@ -37,6 +37,7 @@ IF (THIRDPARTY_BUILD_TINYXML)
TMP_DIR ${TPBUILD}/tinyxml-2.6.2-tmp
......@@ -21,7 +21,12 @@ IF( NEKTAR_USE_VTK )
BINARY_DIR ${TPBUILD}/vtk-5.10.1
TMP_DIR ${TPBUILD}/vtk-5.10.1-tmp
SET(VTK_DIR ${TPDIST}/lib/vtk-5.10)
......@@ -43,6 +43,7 @@ IF (THIRDPARTY_BUILD_ZLIB)
TMP_DIR ${TPBUILD}/zlib-1.2.7-tmp
......@@ -52,6 +53,10 @@ IF (THIRDPARTY_BUILD_ZLIB)
IF (WIN32)
"Zlib library" FORCE)
"Zlib library" FORCE)
"Zlib library" FORCE)
\chapter{Data Structures and Algorithms}
\ No newline at end of file
\ No newline at end of file
......@@ -38,3 +38,24 @@ When using a fully-Fourier expansion, specifies the number of processes to use i
\lstinline[style=BashInputStyle]{--npz [int]}\\
When using Fourier expansions, specifies the number of processes to use in the z-coordinate direction.
Prints detailed information about the generated partitioning, such as number of
elements, number of local degrees of freedom and the number of boundary degrees
of freedom.
\lstinline[style=BashInputStyle]{--part-only [int]}\\
Partition the mesh only into the specified number of partitions, write to file
and exit. This can be used to pre-partition a very large mesh on a single
high-memory node, prior to being executed on a multi-node cluster.
Forces the use of METIS for mesh partitioning. If \nekpp{} is compiled with
Scotch support, the default is to use Scotch.
Forces the use of Scotch for mesh partitioning.
......@@ -29,7 +29,7 @@ There are two ways to obtain the source code for \nekpp:
\item Using authenticated access. This will allow you to directly contribute
back into the code.
git clone nektar++
git clone nektar++
You can easily switch to using the authenticated access from anonymous
......@@ -71,15 +71,21 @@ CardiacEPSolver session.xml
\section{Session file configuration}
\subsection{Solver Info}
\item \inltt{Eqtype} Specifies the PDE system to solve. This should take
the value \inltt{Monodomain}.
\item \inltt{Eqtype} Specifies the PDE system to solve. The following values
are supported:
\item \inltt{Monodomain}: solve the monodomain equation.
\item \inltt{BidomainRoth}: solve the bidomain equations using the Roth
\item \inltt{CellModel} Specifies the cell model to use. Available cell
models are
Value & Description & No. of Var. & Ref. \\
\inltt{AlievPanfilov} & Phenomological & 1 & \cite{AlPa96} \\
\inltt{CourtemancheRamirezNattel98} & Human atrial & 20 & \cite{CoRaNa98} \\
\inltt{FitzHughNagumo} & & & \\
......@@ -88,7 +94,7 @@ CardiacEPSolver session.xml
\inltt{PanditGilesDemir03} & & & \\
\inltt{TenTusscher06} & Human ventricular & 18 & \cite{TuPa06} \\
\inltt{Winslow99} & & & \\
......@@ -154,7 +160,7 @@ $\mathbf{\sigma}$ of the tissue.
\item \inltt{CheckpointCellModel} checkpoints the cell model. Can be
used along with the \inltt{Checkpoint Filter} to record complete simulation
used along with the \inltt{Checkpoint} filter to record complete simulation
state and regular intervals.
<FILTER TYPE="CheckpointCellModel">
......@@ -190,6 +196,30 @@ $\mathbf{\sigma}$ of the tissue.
\item \inltt{Points} specifies a list of coordinates at which electrograms
are desired. \emph{They must not lie within the domain.}
\item \inltt{Benchmark} Records spatially distributed event times for
activation and repolarisation (recovert) during a simulation, for
undertaking benchmark test problems.
<FILTER TYPE="Benchmark">
<PARAM NAME="ThresholdValue"> -40.0 </PARAM>
<PARAM NAME="InitialValue"> 0.0 </PARAM>
<PARAM NAME="OutputFile"> benchmark </PARAM>
<PARAM NAME="StartTime"> 0.0 </PARAM>
\item \inltt{ThresholdValue} specifies the value above which tissue is
considered to be depolarised and below which is considered
\item \inltt{InitialValue} specifies the initial value of the
activation or repolarisation time map.
\item \inltt{OutputFile} specifies the base filename of activation and
repolarisation maps output from the filter. This name is appended
with the index of the event and the suffix `.fld`.
\item \inltt{StartTime} (optional) specifies the simulation time at
which to start detecting events.
......@@ -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.
\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.
\textbf{Format} & \textbf{Extension} & \textbf{High-order} & \textbf{Notes}\\
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
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:
Multiple elements in composite detected; remapped:
- Tag 150 => 150 (Triangle), 151 (Quadrilateral)
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}
MeshConvert -m jac Mesh.xml output.xml
To get a detailed list of elements which have negative Jacobians, one may use
the \inltt{list} option:
MeshConvert -m jac:list Mesh.xml output.xml
and to extract the elements for the purposes of visualisation within the domain,
use the \inltt{extract} boolean parameter:
MeshConvert -m jac:extract Mesh.xml MeshWithNegativeElements.xml
\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.
\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
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}.
......@@ -409,7 +439,9 @@ of 2 and 7 integration points per element use the following command:
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
\subsubsection{High-order cylinder generation}
......@@ -466,6 +498,24 @@ module:
MeshConvert -m detect volume.xml volumeWithBoundaryComposite.xml
\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
MeshConvert -m scalar:surf=1:nq=7:scalar=exp\(x*x+y*y\) mesh.xml deformed.xml
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.
%%% Local Variables:
%%% mode: latex
%%% TeX-master: "../user-guide"
......@@ -155,16 +155,13 @@ every 10 timesteps, we use the syntax:
\subsection {Threshold value}
\subsection {ThresholdMax}
The threshold value filter writes a field output containing a variable $m$,
defined by the time at which the solution first exceeds a specified threshold
This filter assumes that the solution field to be examined is the first
specified in the session file.
defined by the time at which the selected variable first exceeds a specified
threshold value. The default name of the output file is the name of the session
with the suffix \inlsh{\_max.fld}. Thresholding is applied based on the first
variable listed in the session by default.
The following parameters are supported:
......@@ -174,12 +171,16 @@ The following parameters are supported:
\textbf{Option name} & \textbf{Required} & \textbf{Default} &
\textbf{Description} \\
\inltt{OutputFile} & \cmark & - &
Prefix of the output filename to which the checkpoints are written.\\
\inltt{OutputFile} & \xmark & \emph{session}\_max.fld &
Output filename to which the threshold times are written.\\
\inltt{ThresholdVar} & \xmark & \emph{first variable name} &
Specifies the variable on which the threshold will be applied.\\
\inltt{ThresholdValue} & \cmark & - &
Specifies the threshold value.\\
\inltt{InitialValue} & \cmark & - &
Specifies the initial time.\\
\inltt{StartTime} & \xmark & 0.0 &
Specifies the time at which to start recording.\\
......@@ -188,13 +189,19 @@ An example is given below:
<FILTER TYPE="ThresholdMax">
<PARAM NAME="OutputFile"> threshold_max.fld </PARAM>
<PARAM NAME="ThresholdVar"> u </PARAM>
<PARAM NAME="ThresholdValue"> 0.1 </PARAM>
<PARAM NAME="InitialValue"> 0.4 </PARAM>
<PARAM NAME="OutputFile"> ThresholdFile </PARAM>
which produces a field file \inltt{ThresholdFile.fld}.
which produces a field file \inlsh{threshold\_max.fld}.
\subsection{ThresholdMin value}
Performs the same function as the \inltt{ThresholdMax} filter but records the
time at which the threshold variable drops below a prescribed value.
\subsection{Modal energy}
......@@ -86,17 +86,15 @@ namespace Nektar
void MeshPartition::PartitionMesh(bool shared)
void MeshPartition::PartitionMesh(int nParts, bool shared)
ASSERTL0(m_comm->GetRowComm()->GetSize() > 1,
"Partitioning only necessary in parallel case.");
ASSERTL0(m_meshElements.size() >= m_comm->GetRowComm()->GetSize(),
ASSERTL0(m_meshElements.size() >= nParts,
"Too few elements for this many processes.");
m_shared = shared;
if (m_weightingRequired) WeightElements();
PartitionGraph(m_mesh, m_localPartition);
PartitionGraph(m_mesh, nParts, m_localPartition);
void MeshPartition::WriteLocalPartition(LibUtilities::SessionReaderSharedPtr& pSession)
......@@ -131,7 +129,7 @@ namespace Nektar
void MeshPartition::WriteAllPartitions(LibUtilities::SessionReaderSharedPtr& pSession)
for (int i = 0; i < m_comm->GetRowComm()->GetSize(); ++i)
for (int i = 0; i < m_localPartition.size(); ++i)
TiXmlDocument vNew;
TiXmlDeclaration * decl = new TiXmlDeclaration("1.0", "utf-8", "");
......@@ -482,11 +480,12 @@ namespace Nektar
void MeshPartition::PrintPartInfo(std::ostream &out)
int nElmt = boost::num_vertices(m_mesh);
int nPart = m_comm->GetRowComm()->GetSize();
int nPart = m_localPartition.size();
out << "# Partition information:" << std::endl;
out << "# No. elements : " << nElmt << std::endl;
out << "# No. partitions: " << nPart << std::endl;
out << "# ID nElmt nLocDof nBndDof" << std::endl;
BoostVertexIterator vertit, vertit_end;
std::vector<int> partElmtCount(nPart, 0);
......@@ -775,6 +774,7 @@ namespace Nektar
void MeshPartition::PartitionGraph(BoostSubGraph& pGraph,
int nParts,
std::vector<BoostSubGraph>& pLocalPartition)
int i;
......@@ -819,7 +819,6 @@ namespace Nektar
// Call Metis and partition graph
int npart = m_comm->GetRowComm()->GetSize();
int vol = 0;
......@@ -831,10 +830,10 @@ namespace Nektar
// Attempt partitioning using METIS.
int ncon = 1;
PartitionGraphImpl(nGraphVerts, ncon, xadj, adjncy, vwgt, vsize, npart, vol, part);
PartitionGraphImpl(nGraphVerts, ncon, xadj, adjncy, vwgt, vsize, nParts, vol, part);
// Check METIS produced a valid partition and fix if not.
CheckPartitions(nParts, part);
if (!m_shared)
// distribute among columns
......@@ -872,7 +871,7 @@ namespace Nektar
// Create boost subgraph for this process's partitions
int nCols = m_comm->GetRowComm()->GetSize();
int nCols = nParts;
for (i = 0; i < nCols; ++i)
......@@ -892,7 +891,7 @@ namespace Nektar
void MeshPartition::CheckPartitions(Array<OneD, int> &pPart)
void MeshPartition::CheckPartitions(int nParts, Array<OneD, int> &pPart)
unsigned int i = 0;
unsigned int cnt = 0;
......@@ -900,7 +899,7 @@ namespace Nektar
bool valid = true;
// Check that every process has at least one element assigned
for (i = 0; i < npart; ++i)
for (i = 0; i < nParts; ++i)
cnt = std::count(pPart.begin(), pPart.end(), i);
if (cnt == 0)
......@@ -918,7 +917,7 @@ namespace Nektar
for (i = 0; i < pPart.num_elements(); ++i)
pPart[i] = i % npart;
pPart[i] = i % nParts;
......@@ -65,7 +65,7 @@ namespace Nektar
LIB_UTILITIES_EXPORT MeshPartition(const SessionReaderSharedPtr& pSession);
LIB_UTILITIES_EXPORT virtual ~MeshPartition();
LIB_UTILITIES_EXPORT void PartitionMesh(bool shared = false);
LIB_UTILITIES_EXPORT void PartitionMesh(int nParts, bool shared = false);
LIB_UTILITIES_EXPORT void WriteLocalPartition(
SessionReaderSharedPtr& pSession);
LIB_UTILITIES_EXPORT void WriteAllPartitions(
......@@ -223,6 +223,7 @@ namespace Nektar
void WeightElements();
void CreateGraph(BoostSubGraph& pGraph);
void PartitionGraph(BoostSubGraph& pGraph,
int nParts,
std::vector<BoostSubGraph>& pLocalPartition);
virtual void PartitionGraphImpl(
......@@ -237,7 +238,7 @@ namespace Nektar
Nektar::Array<Nektar::OneD, int>& part) = 0;
void OutputPartition(SessionReaderSharedPtr& pSession, BoostSubGraph& pGraph, TiXmlElement* pGeometry);
void CheckPartitions(Array<OneD, int> &pPart);
void CheckPartitions(int nParts, Array<OneD, int> &pPart);
typedef boost::shared_ptr<MeshPartition> MeshPartitionSharedPtr;
......@@ -345,6 +345,8 @@ namespace Nektar
"number of procs in Y-dir")
("npz", po::value<int>(),
"number of procs in Z-dir")
("part-only", po::value<int>(),
"only partition mesh into N partitions.")
("part-info", "Output partition information")
......@@ -1330,6 +1332,19 @@ namespace Nektar
"Error: File '" + pFilename + "' is corrupt.");
else if (pFilename.size() > 4 &&
pFilename.substr(pFilename.size() - 4, 4) == "_xml")
fs::path pdirname(pFilename);
boost::format pad("P%1$07d.xml");
pad % m_comm->GetRank();
fs::path pRankFilename(pad.str());
fs::path fullpath = pdirname / pRankFilename;
ifstream file(PortablePath(fullpath).c_str());
ASSERTL0(file.good(), "Unable to open file: " + fullpath.string());
file >> (*pDoc);
ifstream file(pFilename.c_str());
......@@ -1458,28 +1473,97 @@ namespace Nektar
// Get row of comm, or the whole comm if not split
CommSharedPtr vCommMesh = m_comm->GetRowComm();
const bool isRoot = (m_comm->GetRank() == 0);