Commit 162e4be3 authored by Chris Cantwell's avatar Chris Cantwell

Merge branch 'feature/mesh-deform' into 'master'

Add linear elastic solver

This MR adds a solver for the linear elasticity equations and the library and pre/post-processing infrastructure that is needed to support it. In particular, the primary application for the solver within Nektar++ is for mesh generation purposes, whereby the mesh is treated as a solid body and deformed at the boundary to align the curvilinear elements with a given surface.

To do this robustly, as the linear elastic equations are only valid for small deformations, the implementation presented here adopts a multi-step approach, solving the equations after each step. As such, routines are needed in `MultiRegions` and `SpatialDomains` to reset the geometry information and rebuild the matrix system at each step.

The following changes have been made to accommodate this:

- `Geometry` objects now have a `Reset` method, which clears the geometry coefficients of itself and any derived edges/faces, rebuilding this from the vertex and curvature information. To facilitate this, most elements now have a `SetUpXmap` routine which regenerates the `m_xmap` object as curvature may have been added to the element.
- Similarly, `ExpList` also has a `Reset` method which clears the matrix managers and causes matrices to be rebuild on the next call to `HelmSolve` or similar.
- Added support for 1D curved segments in 2D space to `MeshConvert` and `MeshGraph`.
- `FieldConvert` has three new modules:
  - a displacement processing module which calculates the displacement between an XML manifold and a 2D/3D domain's boundary region and outputs a file suitable for a boundary condition in the elasticity solver;
  - a deformation processing module which takes the displacement field generated by the elasticity solver and applies it to the mesh;
  - and an XML output module which can be used in conjunction with the previous module. Since this is likely to be useful elsewhere, most of the code for this has been added as a `MeshGraph` routine.
- `XmlToVtk` has a new command line argument to output the distribution of scaled Jacobians of a given mesh for quality assessment purposes.
- The elasticity code itself which includes the thermal terms discussed in publications.
- Documentation for all of the above.

See merge request !400
parents dbb6cb15 93112dc2
......@@ -342,3 +342,14 @@ year={2011}
pages={95--136},
year={2014}
}
@article{Ko07,
title = {Vectorized Matlab codes for linear two-dimensional elasticity},
author = {Koko, Jonas},
journal = {Scientific Programming},
volume = 15,
number = 3,
pages = {157--172},
year = 2007,
publisher = {IOS Press}
}
\ No newline at end of file
......@@ -353,7 +353,7 @@ the correct solution is obtained by evaluating the $L_2$ and $L_{inf}$ errors.
\subsubsection{Running the code}
\begin{lstlisting}[style=BashInputStyle]
ADRSolver Test\_Helmholtz2D\_modal.xml
ADRSolver Test_Helmholtz2D_modal.xml
\end{lstlisting}
This execution should print out a summary of input file, the $L_2$ and
......@@ -363,7 +363,7 @@ $L_{inf}$ errors and the time spent on the calculation.
Simulation results are written in the file Helmholtz2D\_modal.fld. We can choose
to visualise the output in Gmsh
\begin{lstlisting}[style=BashInputStyle]
FldToGmsh Helmholtz2D\_modal.xml Helmholtz2D\_modal.fld
FldToGmsh Helmholtz2D_modal.xml Helmholtz2D_modal.fld
\end{lstlisting}
which generates the file Helmholtz2D\_modal\_u.pos as shown in
Fig.~\ref{f:adrsolver:helmholtz2D}
......
\chapter{Linear elasticity solver}
\label{s:elasticity}
\section{Synopsis}
The LinearElasticSolver is a solver for solving the linear elasticity equations
in two and three dimensions. Whilst this may be suitable for simple solid
mechanics problems, its main purpose is for use for mesh deformation and
high-order mesh generation, whereby the finite element mesh is treated as a
solid body, and the deformation is applied at the boundary in order to curve the
interior of the mesh.
Currently the following equation systems are supported:
%
\begin{center}
\begin{tabular}{lp{8cm}}
\toprule
Value & Description \\
\midrule
\inltt{LinearElasticSystem} & Solves the linear elastic equations. \\
\inltt{IterativeElasticSystem} & A multi-step variant of the elasticity solver,
which breaks a given deformation into multiple
steps, and applies the deformation to a mesh. \\
\bottomrule
\end{tabular}
\end{center}
\subsection{The linear elasticity equations}
The linear elasticity equations model how a solid body deforms under the
application of a `small' deformation or strain. The formulation starts with the
equilibrium of forces represented by the equation
%
\begin{equation}
\nabla \cdot \mathbf{S} + \mathbf{f} = \mathbf{0} \quad \textrm{in} \quad \Omega
\label{eq:strong}
\end{equation}
%
where $\mathbf{S}$ is the stress tensor and $\mathbf{f}$ denotes a
spatially-varying force. We further assume that the stress tensor $\mathbf{S}$
incorporates elastic and, optionally, thermal stresses that can be switched on
to assist in mesh deformation applications. We assume these can be decomposed so
that $\mathbf{S}$ is written as
%
\[
\mathbf{S} = \mathbf{S}_e + \mathbf{S}_t,
\]
%
where the subscripts $e$ and $t$ denote the elastic and thermal terms
respectively. We adopt the usual linear form of the elastic stress tensor as
%
\[
\mathbf{S}_e = \lambda\mbox{Tr}(\mathbf{E}) \, \mathbf{I} +\mu \mathbf{E},
\]
%
where $\lambda$ and $\mu$ are the Lam\'e constants, $\mathbf{E}$ represents the
strain tensor, and $\mathbf{I}$ is the identity tensor. For small deformations,
the strain tensor $\mathbf{E}$ is given as
%
\begin{equation}
\mathbf{E} =\frac{1}{2} \left ( \nabla \mathbf{u}+ \nabla \mathbf{u}^t \right )
\end{equation}
%
where $\mathbf{u}$ is the two- or three-dimensional vector of displacements. The
boundary conditions required to close the problem consist of prescribed
displacements at the boundary $\partial \Omega$, i.e.
\begin{equation}
\mathbf{u} = \hat{\mathbf{u}} \quad \textrm{in}\ \partial \Omega.
\end{equation}
We further express the Lam\'e constants in terms of the Young's modulus $E$ and
Poisson ratio $\nu$ as
%
\[
\lambda = \frac{\nu E}{(1+\nu)(1-2\nu)}, \qquad \mu = \frac{E}{2(1+\nu)}.
\]
%
The Poisson ratio, valid in the range $\nu < \tfrac{1}{2}$, is a measure of the
compressibility of the body, and the Young's modulus $E > 0$ is a measure of its
stiffness.
\section{Usage}
\begin{lstlisting}[style=BashInputStyle]
LinearElasticSolver [arguments] session.xml [another.xml] ...
\end{lstlisting}
\section{Session file configuration}
\subsection{Solver Info}
\begin{itemize}
\item \inltt{EqType} Specifies the PDE system to solve, based on the choices
in the table above.
\item \inltt{Temperature} Specifies the form of the thermal stresses to
use. The choices are:
\begin{itemize}
\item \inltt{None}: No stresses (default).
\item \inltt{Jacobian}: Sets $\mathbf{S}_t = \beta JI$, where $\beta$ is a
parameter defined in the parameters section, $J$ is the elemental Jacobian
determinant and $I$ is the identity matrix.
\item \inltt{Metric}: A more complex term, based on the eigenvalues of the
metric tensor. This can only be used for simplex elements (triangles and
tetrahedra). Controlled again by the parameter $\beta$.
\end{itemize}
\item \inltt{BCType} Specifies the type of boundary condition to apply when
the \inltt{IterativeElasticSystem} is being used.
\begin{itemize}
\item \inltt{Normal}: The boundary conditions are split into
\inltt{NumSteps} steps, as defined by a parameter in the session file
(default).
\item \inltt{Repeat}: As the geometry is updated, re-evaluate the boundary
conditions. This enables, for example, a cirlce to be rotated continuously.
\end{itemize}
\end{itemize}
\subsection{Parameters}
The following parameters can be specified in the \inltt{PARAMETERS} section of
the session file:
\begin{itemize}
\item \inltt{nu}: sets the Poisson ratio $\nu$.\\
\textit{Default value}: 0.25.
\item \inltt{E}: sets the Young's modulus $E$.\\
\textit{Default value}: 1.
\item \inltt{beta}: sets the thermal stress coefficient $\beta$.\\
\textit{Default value}: 1.
\item \inltt{NumSteps}: sets the number of steps to use in the case that the
iterative elastic system is enabled. Should be greater than 0.\\
\textit{Default value}: 0.
\end{itemize}
\section{Examples}
\subsection{L-shaped domain}
The first example is the classic L-shaped domain, in which an exact solution is
known, which makes it an ideal test case~\cite{Ko07}. The domain is the polygon
formed from the vertices
%
\[
(-1,-1), (0,-2), (2,0), (0,2), (-1,-1), (0,0).
\]
%
The exact solution for the displacements is known in polar co-ordinates
$(r,\theta)$ as
%
\begin{align*}
u_r(r,\theta) &= \frac{r^\alpha}{2\mu} \left[
C_1(C_2 - \alpha - 1)\cos((\alpha-1)\theta) - (\alpha+1)\cos((\alpha+1)\theta)
\right]\\
u_\theta(r,\theta) &= \frac{r^\alpha}{2\mu} \left[
(\alpha+1)\sin((\alpha+1)\theta) + C_1(C_2+\alpha-1)\sin((\alpha-1)\theta)
\right]
\end{align*}
%
where $\alpha\approx 0.544483737\dots$ is the solution of the equation
$\alpha\sin(2\omega) + \sin(2\omega\alpha) = 0$,
%
\[
C_1 = -\frac{\cos((\alpha+1)\omega)}{\cos((\alpha-1)\omega)},\qquad
C_2 = 2\frac{\lambda + 2\mu}{\lambda+\mu}
\]
with $\lambda$ and $\mu$ being the Lam\'e constants and $\omega = 3\pi/4$.
Boundary conditions are set to be the exact solution and
$\mathbf{f} = \mathbf{0}$. The solution has a singularity at the origin, and so
in order to test convergence $h$-refinement is required.
\begin{figure}
\begin{center}
\includegraphics[width=0.5\textwidth]{Figures/l-shape}
\end{center}
\caption{Solution of the $u$ displacement field for the L-shaped domain.}
\label{fig:elas:ldomain}
\end{figure}
A simple example of how the linear elastic solver can be set up can be found in
the \inlsh{Tests/L-shaped.xml} session file in the linear elastic solver
directory. A more refined domain with the obtained $u$ solution is shown in
figure~\ref{fig:elas:ldomain}. The solver can be run using the command:
\begin{lstlisting}[style=BashInputStyle]
LinearElasticSolver L-domain.xml
\end{lstlisting}
The obtained solution \inlsh{L-domain.fld} can be applied to the mesh to obtain
a deformed XML file using the \inltt{deform} module in \inltt{FieldConvert}:
\begin{lstlisting}[style=BashInputStyle]
FieldConvert -m deform L-domain.xml L-domain.fld L-domain-deformed.xml
\end{lstlisting}
\subsection{Boundary layer deformation}
In this example we use the iterative elastic system to apply a large deformation
to a triangular boundary layer mesh of a square mesh $\Omega = [0,1]^2$. At the
bottom edge, we apply a Dirichlet condition $g=\tfrac{1}{2}\sin(\pi x)$ that is
enforced by splitting it into $N$ substeps, so that at each step we solve the
system with the boundary condition $g^n(x) = g(x)/N$. The process is depicted in
figure~\ref{fig:elas:bl}.
\begin{figure}
\begin{center}
\includegraphics[width=0.3\textwidth]{Figures/bl-0}
\includegraphics[width=0.3\textwidth]{Figures/bl-1}
\includegraphics[width=0.3\textwidth]{Figures/bl-2}
\end{center}
\caption{Figures that show the initial domain (left), after 50 steps (middle)
and final deformation of the domain (right).}
\label{fig:elas:bl}
\end{figure}
The setup is very straightforward. The geometry can be found inside the file
\inlsh{Examples/bl-mesh.xml} and the conditions inside
\inlsh{Examples/bl-conditions.xml}. The solver can be set up using the following
parameters, with \inltt{NumSteps} denoting $N$:
\begin{lstlisting}[style=XMLStyle]
<SOLVERINFO>
<I PROPERTY="EQTYPE" VALUE="IterativeElasticSystem" />
</SOLVERINFO>
<PARAMETERS>
<P> nu = 0.3 </P>
<P> E = 1.0 </P>
<P> NumSteps = 100 </P>
</PARAMETERS>
\end{lstlisting}
To identify the boundary that we intend to split up into substeps, we must
assign the \inltt{WALL} tag to our boundary regions:
\begin{lstlisting}[style=XMLStyle]
<BOUNDARYCONDITIONS>
<REGION REF="0">
<D VAR="u" VALUE="0" USERDEFINEDTYPE="Wall" />
<D VAR="v" VALUE="0.5*sin(PI*x)" USERDEFINEDTYPE="Wall" />
</REGION>
<REGION REF="1">
<D VAR="u" VALUE="0" />
<D VAR="v" VALUE="0" />
</REGION>
</BOUNDARYCONDITIONS>
\end{lstlisting}
The solver can then be run using the command:
\begin{lstlisting}[style=BashInputStyle]
LinearElasticSolver bl-mesh.xml bl-conditions.xml
\end{lstlisting}
This will produce a series of meshes \inlsh{bl-mesh-\%d.xml}, where \%d is an
integer running between 0 and 100. If at any point the mesh becomes invalid,
that is, a negative Jacobian is detected, execution will cease.
%%% Local Variables:
%%% mode: latex
%%% TeX-master: "../user-guide"
%%% End:
......@@ -527,7 +527,7 @@ Finally the domain is specified by the first composite by
</SOLVERINFO>
\end{lstlisting}
\paragraph{Parameters:~} Parameters used for the simulation are taken from \cite{Sherwin2003}
\paragraph{Parameters:~} Parameters used for the simulation are taken from \cite{ShFoPeFr03}
\begin{lstlisting}[style=XMLStyle]
<PARAMETERS>
<P> TimeStep = 2e-6 </P>
......
......@@ -6,5 +6,6 @@
\input{cardiac-ep}
\input{compressible-flow}
\input{incompressible-ns}
\input{linear-elastic}
\input{pulse-wave}
\input{shallow-water}
\section{Field Convert}
\section{FieldConvert}
\label{s:utilities:fieldconvert}
FieldConvert is a utility embedded in \nekpp with the primary
aim of allowing the user to convert the \nekpp output binary files
......@@ -53,7 +53,7 @@ in its compressed format \inltt{test.xml.gz}.
%
%
%
\subsection{Field Convert range option \textit{-r}}
\subsection{Range option \textit{-r}}
The Fieldconvert range option \inltt{-r} allows the user to specify
a sub-range of the mesh (computational domain) by using an
additional flag, \inltt{-r} (which stands for \inltt{r}ange and either
......@@ -77,7 +77,7 @@ FieldConvert -r 2,3,-1,2 test.xml test.fld test.dat
\end{lstlisting}
%
\end{itemize}
where \inltt{-r} defines the range option of the field convert
where \inltt{-r} defines the range option of the FieldConvert
utility, the two first numbers define the range in $x$ direction
and the the third and fourth number specify the $y$ range.
A sub-range of a 3D domain can also be specified.
......@@ -86,7 +86,7 @@ to define the $z$ range.
%
%
%
\subsection{Field Convert modules \textit{-m}}
\subsection{FieldConvert modules \textit{-m}}
FieldConvert allows the user to manipulate the \nekpp output
binary files (.chk and .fld) by using the flag \inltt{-m} (which
stands for \inltt{m}odule)..
......@@ -423,7 +423,50 @@ way as described in section \ref{s:utilities:fieldconvert:sub:convert}.
%
%
%
\subsection{Field Convert in parallel}
\subsubsection{Manipulating meshes with FieldConvert}
FieldConvert has support for two modules that can be used in conjunction with
the linear elastic solver, as shown in chapter~\ref{s:elasticity}. To do this,
FieldConvert has an XML output module, in addition to the Tecplot and VTK
formats.
The \inltt{deform} module, which takes no options, takes a displacement field
and applies it to the geometry, producing a deformed mesh:
\begin{lstlisting}[style=BashInputStyle]
FieldConvert -m deform input.xml input.fld deformed.xml
\end{lstlisting}
The \inltt{displacement} module is designed to create a boundary condition field
file. Its intended use is for mesh generation purposes. It can be used to
calculate the displacement between the linear mesh and a high-order surface, and
then produce a \inltt{fld} file, prescribing the displacement at the boundary,
that can be used in the linear elasticity solver.
Presently the process is somewhat convoluted and must be used in conjunction
with MeshConvert to create the surface file. However the bash input below
describes the procedure. Assume the high-order mesh is in a file called
\inlsh{mesh.xml}, the linear mesh is \inlsh{mesh-linear.xml} that can be
generated by removing the \inltt{CURVED} section from \inlsh{mesh.xml}, and that
we are interested in the surface with ID 123.
\begin{lstlisting}[style=BashInputStyle]
# Extract high order surface
MeshConvert -m extract:surf=123 mesh.xml mesh-surf-curved.xml
# Use FieldConvert to calculate displacement between two surfaces
FieldConvert -m displacement:id=123:to=mesh-surf-curved.xml \
mesh-linear.xml mesh-deformation.fld
# mesh-deformation.fld is used as a boundary condition inside the
# solver to prescribe the deformation conditions.xml contains
# appropriate Nektar++ parameters (mu, E, other BCs, ...)
LinearElasticSolver mesh-linear.xml conditions.xml
# This produces the final field mesh-linear.fld which is the
# displacement field, use FieldConvert to apply it:
FieldConvert-g -m deform mesh-linear.xml mesh-linear.fld mesh-deformed.xml
\end{lstlisting}
\subsection{FieldConvert in parallel}
To run FieldConvert in parallel the user needs to compile
\nekpp with MPI support and can employ the following
command
......@@ -477,3 +520,8 @@ Note the final \inltt{:info} extension on the last argument is necessary to
tell FieldConvert that you wish to generate an info file but the extension
to this file is .xml. This syntax allows the routine not to get confused with
the input/output xml files.
%%% Local Variables:
%%% mode: latex
%%% TeX-master: "../user-guide"
%%% End:
......@@ -390,7 +390,7 @@ MeshConvert -m peralign:surf1=11:surf2=12:dir=y:orient input.dat output.xml
Often it is the case that one can generate a coarse boundary layer grid of a
mesh. \mc has a method for splitting prismatic and hexahedral elements into
finer elements based on the work presented in~\cite{MoHaJoSh14}
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.
......
......@@ -192,15 +192,15 @@ SQRT1\_2 & $\frac{1}{\sqrt{2}}$ & 0.70710678118654752440 \\
\end{tabular}
\end{center}
\item parameters: alphanumeric names with underscores, e.g. \inltt{GAMMA\_123,
GaM123\_45a\_, \_gamma123} are perfectly acceptable parameter names. However
parameter name cannot start with a numeral. Parameters must be defined with
\inltt{<PARAMETERS>...</PARAMETERS>}. Parameters play the role of constants
that may change their values in between of expression evaluations.
\item parameters: alphanumeric names with underscores, e.g. \inltt{GAMMA\_123},
\inltt{GaM123\_45a\_}, \inltt{\_gamma123} are perfectly acceptable parameter
names. However parameter name cannot start with a numeral. Parameters must be
defined with \inltt{<PARAMETERS>...</PARAMETERS>}. Parameters play the role of
constants that may change their values in between of expression evaluations.
\item variables (i.e., \inlsh{x, y, z} and \inlsh{t})
\item unary minus operator (e.g. \inltt{-x}
\item binary arithmetic operators \inltt{+, -, *, /, \^}
\item binary arithmetic operators \inltt{+, -, *, /, \^{}}
Powering operator allows using real exponents (it is implemented with
\inlsh{std::pow()} function)
\item boolean comparison operations \inlsh{<, <=, >, >=, ==} evaluate their
......@@ -239,7 +239,7 @@ These functions are implemented by means of the cmath library:
\url{http://www.cplusplus.com/reference/clibrary/cmath/}. Underlying data type
is \inltt{double} at each stage of expression evaluation. As consequence,
complex-valued expressions (e.g. $(-2)^0.123$) get value \inlsh{nan} (not a
number). The operator \inlsh{\^} is implemented via call to \inlsh{std::pow()}
number). The operator \inlsh{\^{}} is implemented via call to \inlsh{std::pow()}
function and accepts arbitrary real exponents.
\item random noise generation functions. Currently implemented is
......@@ -265,9 +265,9 @@ costs.
\subsubsection{Pre-evaluation details}
Partial evaluation of all constant sub-expressions makes no sense in using
derived constants from table above. This means, either make use of pre-defined
constant \inlsh{LN10\^2} or straightforward expression \inlsh{log10(2)\^2}
results in constant \inlsh{5.3018981104783980105} being stored internally
after pre-processing. The rules of pre-evaluation are as follows:
constant \inlsh{LN10\^{}2} or straightforward expression \inlsh{log10(2)\^{}2}
results in constant \inlsh{5.3018981104783980105} being stored internally after
pre-processing. The rules of pre-evaluation are as follows:
\begin{itemize}
\item constants, numbers and their combinations with arithmetic, analytic and
comparison operators are pre-evaluated,
......@@ -304,21 +304,21 @@ increasing complexity:
\begin{itemize}
\item \inlsh{+, -, <, >, <=, >=, ==, }
\item \inlsh{*, /, abs, fabs, ceil, floor,}
\item \inlsh{\^, sqrt, exp, log, log10, sin, cos, tan, sinh, cosh, tanh, asin,
\item \inlsh{\^{}, sqrt, exp, log, log10, sin, cos, tan, sinh, cosh, tanh, asin,
acos, atan}.
\end{itemize}
For example,
\begin{itemize}
\item \inlsh{x*x} is faster than \inlsh{x\^2} --- it is one double
\item \inlsh{x*x} is faster than \inlsh{x\^{}2} --- it is one double
multiplication vs generic calculation of arbitrary power with real exponents.
\item \inlsh{(x+sin(y))\^2} is faster than \inlsh{(x+sin(y))*(x+sin(y))} -
\item \inlsh{(x+sin(y))\^{}2} is faster than \inlsh{(x+sin(y))*(x+sin(y))} -
sine is an expensive operation. It is cheaper to square complicated expression rather than
compute it twice and add one multiplication.
\item An expression
\inltt{exp(-41*( (x+(0.3*cos(2*PI*t)))\^2 + (0.3*sin(2*PI*t))\^2 ))}
\inltt{exp(-41*( (x+(0.3*cos(2*PI*t)))\^{}2 + (0.3*sin(2*PI*t))\^{}2 ))}
makes use of 5 expensive operations (\inlsh{exp}, \inlsh{sin}, \inlsh{cos}
and power \inlsh{\^} twice) while an equivalent expression
and power \inlsh{\^{}} twice) while an equivalent expression
\inltt{exp(-41*( x*x+0.6*x*cos(2*PI*t) + 0.09 ))}
uses only 2 expensive operations.
\end{itemize}
......
......@@ -286,7 +286,6 @@ namespace Nektar
{
TiXmlElement* x;
TiXmlElement *vGeometry, *vSubElement;
int i;
vGeometry = pSession->GetElement("Nektar/Geometry");
m_dim = atoi(vGeometry->Attribute("DIM"));
......@@ -297,7 +296,7 @@ namespace Nektar
// Retrieve any VERTEX attributes specifying mesh transforms
std::string attr[] = {"XSCALE", "YSCALE", "ZSCALE",
"XMOVE", "YMOVE", "ZMOVE" };
for (i = 0; i < 6; ++i)
for (int i = 0; i < 6; ++i)
{
const char *val = vSubElement->Attribute(attr[i].c_str());
if (val)
......@@ -307,18 +306,13 @@ namespace Nektar
}
x = vSubElement->FirstChildElement();
i = 0;
if (x->FirstAttribute())
{
i = x->FirstAttribute()->IntValue();
}
while(x)
{
TiXmlAttribute* y = x->FirstAttribute();
ASSERTL0(y, "Failed to get attribute.");
MeshVertex v;
v.id = y->IntValue();
ASSERTL0(v.id == i++, "Vertex IDs not sequential.");
std::vector<std::string> vCoords;
std::string vCoordStr = x->FirstChild()->ToText()->Value();
boost::split(vCoords, vCoordStr, boost::is_any_of("\t "));
......@@ -335,7 +329,6 @@ namespace Nektar
vSubElement = pSession->GetElement("Nektar/Geometry/Edge");
ASSERTL0(vSubElement, "Cannot read edges");
x = vSubElement->FirstChildElement();
i = 0;
while(x)
{
TiXmlAttribute* y = x->FirstAttribute();
......@@ -343,7 +336,6 @@ namespace Nektar
MeshEntity e;
e.id = y->IntValue();
e.type = 'E';
ASSERTL0(e.id == i++, "Edge IDs not sequential.");
std::vector<std::string> vVertices;
std::string vVerticesString = x->FirstChild()->ToText()->Value();
boost::split(vVertices, vVerticesString, boost::is_any_of("\t "));
......@@ -360,7 +352,6 @@ namespace Nektar
vSubElement = pSession->GetElement("Nektar/Geometry/Face");
ASSERTL0(vSubElement, "Cannot read faces.");
x = vSubElement->FirstChildElement();
i = 0;
while(x)
{
TiXmlAttribute* y = x->FirstAttribute();
......@@ -368,7 +359,6 @@ namespace Nektar
MeshEntity f;
f.id = y->IntValue();
f.type = x->Value()[0];
ASSERTL0(f.id == i++, "Face IDs not sequential.");
std::vector<std::string> vEdges;
std::string vEdgeStr = x->FirstChild()->ToText()->Value();
boost::split(vEdges, vEdgeStr, boost::is_any_of("\t "));
......@@ -385,14 +375,12 @@ namespace Nektar
vSubElement = pSession->GetElement("Nektar/Geometry/Element");
ASSERTL0(vSubElement, "Cannot read elements.");
x = vSubElement->FirstChildElement();
i = 0;
while(x)
{
TiXmlAttribute* y = x->FirstAttribute();
ASSERTL0(y, "Failed to get attribute.");
MeshEntity e;
e.id = y->IntValue();
ASSERTL0(e.id == i++, "Element IDs not sequential.");
std::vector<std::string> vItems;
std::string vItemStr = x->FirstChild()->ToText()->Value();
boost::split(vItems, vItemStr, boost::is_any_of("\t "));
......@@ -410,7 +398,6 @@ namespace Nektar
{
vSubElement = pSession->GetElement("Nektar/Geometry/Curved");
x = vSubElement->FirstChildElement();
i = 0;
while(x)
{
MeshCurved c;
......@@ -446,7 +433,6 @@ namespace Nektar
vSubElement = pSession->GetElement("Nektar/Geometry/Composite");
ASSERTL0(vSubElement, "Cannot read composites.");
x = vSubElement->FirstChildElement();
i = 0;
while(x)
{
TiXmlAttribute* y = x->FirstAttribute();
......@@ -497,7 +483,7 @@ namespace Nektar
std::map<int, int> elmtSizes;
std::map<int, int> elmtBndSizes;
for (unsigned int i = 0; i < m_domain.size(); ++i)
{
int cId = m_domain[i];
......@@ -528,7 +514,7 @@ namespace Nektar
for (unsigned int j = 0; j < m_meshComposites[cId].list.size(); ++j)
{
int elid = m_meshComposites[cId].list[j];
int elid = m_meshComposites[cId].list[j];
elmtSizes[elid] = weight;
elmtBndSizes[elid] = bndWeight;
}
......@@ -622,9 +608,10 @@ namespace Nektar
void MeshPartition::WeightElements()
{