Commit 8b340823 authored by Dave Moxey's avatar Dave Moxey

Merge remote-tracking branch 'upstream/master' into feature/python-array-creation

parents e135a5c0 5a5ebfca
......@@ -58,6 +58,7 @@ v5.0.0
- Fix compilation issue with newer Boost versions and clang (!940)
- If only `NEKTAR_BUILD_LIBRARY` is enabled, only libraries up to and including
`MultiRegions` will be built by default (!945)
- Fix missing metadata import from Hdf5 files (!971)
**NekMesh**:
- Add feature to read basic 2D geo files as CAD (!731)
......@@ -89,6 +90,7 @@ v5.0.0
- Fix manifold face curvature nodes (!913)
- Fix writing 1D surfaces (!930)
- Fix surface string parsin in BL splitting (!937)
- Add support for Gmsh 4.0 mesh file format (!964)
**FieldConvert**:
- Add input module for Semtex field files (!777)
......@@ -105,6 +107,7 @@ v5.0.0
- Add module for modifying/adding fields from expressions (!889, !903)
- Add module for evaluating the mean of variables on the domain (!894)
- Add module for counting the total number of DOF (!948)
- Fixed wss module for compressible flows (!958)
**IncNavierStokesSolver**
- Replace steady-state check based on difference of norms by check based on
......@@ -122,13 +125,23 @@ v5.0.0
- Fix compressible solver with NUMMODES=1 (!868)
- Introduce equations of state to account for real gas effects (!880)
**APESolver:**
**AcousticSolver:**
- Added two new boundary conditions to the APE system: RiemannInvariantBC
and WhiteNoise (!782)
- Store base flow fields in a discontinuous projection (!918)
- Enabled 1D cases (!918)
- The APE system now uses u_i, c^2 and rho as base flow fields (!918)
- Added the Linearized Euler Equations (LEE) (!918)
**APESolver:**
- APESolver was replaced with AcousticSolver (!918)
**Documentation**:
- Added the developer-guide repository as a submodule (!751)
**Tester**
- Fix build with boost 1.67 (!947)
v4.4.2
------
**Library**
......@@ -160,6 +173,7 @@ v4.4.2
**FieldConvert**
- Allow passing input name with trailing separator (!879)
- Fix the interpcoord option of the interppointdatatofld module (!952)
**Utilities**
- Fix VtkToPng to account for deprecated VTK API for VTK version > 8.1 (!925)
......
......@@ -5,6 +5,7 @@ RelWithDebInfo MinSizeRel.")
PROJECT(Nektar++ C CXX)
CMAKE_POLICY(SET CMP0022 NEW)
# Nektar++ requires C++11. Try to infer this for older CMake versions (less than
# 3.1.0)
IF ("${CMAKE_VERSION}" VERSION_LESS "3.1")
......
......@@ -33,7 +33,8 @@ IF( NEKTAR_USE_VTK )
SET(VTK_USE_FILE ${VTK_DIR}/UseVTK.cmake)
INCLUDE (${VTK_DIR}/VTKConfig.cmake)
ELSE()
FIND_PACKAGE(VTK)
FIND_PACKAGE(VTK COMPONENTS
vtkFiltersGeometry vtkIOLegacy vtkIOXML vtkIOImage vtkRenderingCore)
IF (VTK_FOUND)
MESSAGE(STATUS "Found VTK: ${VTK_USE_FILE}")
IF (VTK_MAJOR_VERSION EQUAL 6 AND VTK_MINOR_VERSION EQUAL 0 AND VTK_BUILD_VERSION EQUAL 0)
......
......@@ -53,16 +53,24 @@ IF (THIRDPARTY_BUILD_ZLIB)
IF (WIN32)
THIRDPARTY_LIBRARY(ZLIB_LIBRARIES STATIC zlib DESCRIPTION "Zlib library")
THIRDPARTY_LIBRARY(ZLIB_LIBRARIES_DEBUG STATIC zlibd DESCRIPTION "Zlib library")
ELSE ()
THIRDPARTY_LIBRARY(ZLIB_LIBRARIES SHARED z DESCRIPTION "Zlib library")
THIRDPARTY_LIBRARY(ZLIB_LIBRARIES_DEBUG SHARED z DESCRIPTION "Zlib library")
ENDIF ()
MESSAGE(STATUS "Build Zlib: ${ZLIB_LIBRARIES}")
MESSAGE(STATUS "Build Zlib: ")
MESSAGE(STATUS " -- Optimized: ${ZLIB_LIBRARIES}")
MESSAGE(STATUS " -- Debug: ${ZLIB_LIBRARIES_DEBUG}")
SET(ZLIB_INCLUDE_DIR ${TPDIST}/include CACHE PATH "Zlib include" FORCE)
SET(ZLIB_CONFIG_INCLUDE_DIR ${TPINC})
ELSE (THIRDPARTY_BUILD_ZLIB)
ADD_CUSTOM_TARGET(zlib-1.2.7 ALL)
MESSAGE(STATUS "Found Zlib: ${ZLIB_LIBRARIES} (version ${ZLIB_VERSION_STRING})")
# We use the found library also for debug builds.
SET(ZLIB_LIBRARIES_DEBUG ${ZLIB_LIBRARIES})
SET(ZLIB_CONFIG_INCLUDE_DIR ${ZLIB_INCLUDE_DIRS})
ENDIF (THIRDPARTY_BUILD_ZLIB)
......
......@@ -507,10 +507,35 @@ year = {2016}
publisher={Elsevier}
}
@inproceedings{LaSaJa17,
author = {Lackhove, Kilian and Sadiki, Amsini and Janicka, Johannes},
booktitle = {Proceedings of ASME Turbo Expo 2007},
keywords = {GT2017-63050},
title = {{Efficient Three Dimensional Time-Domain Combustion Noise Simulation of a Premixed Flame Using Acoustic Perturbation Equations and Incompressible LES}},
year = {2017}
@article{EwSc03,
author = {Ewert, Roland and Schr{\"{o}}der, Wolfgang},
title = {{Acoustic perturbation equations based on flow decomposition via source filtering}},
journal = {Journal of Computational Physics},
year = {2003},
volume = {188},
number = {2},
month = {7},
pages = {365--398},
issn = {00219991},
doi = {10.1016/S0021-9991(03)00168-2},
url = {http://linkinghub.elsevier.com/retrieve/pii/S0021999103001682},
keywords = {APE,acoustic perturbation equations,caa,flow decomposition,les,source filtering},
}
@InProceedings{Ge14,
author = {Geiser, Georg and Nawroth, Holger and Hosseinzadeh, Arash and Zhang, Feichi and Bockhorn, Henning and Habisreuther, Peter and Janicka, Johannes and Paschereit, Christian O. and Schröder, Wolfgang},
title = {Thermoacoustics of a turbulent premixed flame},
booktitle = {20th AIAA/CEAS Aeroacoustics Conference},
date = {2014-06},
publisher = {American Institute of Aeronautics and Astronautics},
doi = {10.2514/6.2014-2476},
}
@phdthesis{La18,
author = {Lackhove, Kilian},
title = {{Hybrid Noise Simulation for Enclosed Configurations}},
year = {2018},
type = {Doctoral Thesis},
school = {Technische Universität Darmstadt},
url = {http://tuprints.ulb.tu-darmstadt.de/id/eprint/7611},
}
This diff is collapsed.
\chapter{Acoustic Perturbation Equations Solver}
\section{Synopsis}
The aim of APESolver is to predict aerodynamic sound generation. Through
the application of a splitting technique, the flow-induced acoustic field is
totally decoupled from the underlying incompressible hydrodynamic field. The
acoustic perturbation equations (APE-1/APE-4) proposed by Ewert and Shroeder are employed as
the governing equations of the acoustic field and they assure stable
aeroacoustic simulation due to the suppression of the term related to the
production of perturbed vorticity. These equations are similar to the linearised
perturbed compressible equations, but while in the original formulation the flow
decomposition is based on solenoidal vortical perturbations as well as
irrotational acoustic perturbations, in this case perturbations are assumed to
be exclusively of acoustic nature.
\begin{subequations}
\begin{align*}
\frac{\partial p'}{\partial t}
+ \frac{\partial \gamma \overline{p} u'_i}{\partial x_i}
+ \frac{\partial \overline{u}_i p'}{\partial x_i}
&= \overline{c}^2 q_c
\\
\frac{\partial u'_i}{\partial t}
+ \frac{\partial \overline{u}_j u'_j}{\partial x_i}
+ \frac{\partial p' / \overline{\rho}}{\partial x_i}
&= U_i
\end{align*}
\end{subequations}
where $(\overline{u}_i,\overline{p}, \overline{\rho}, \overline{c}^2 = \gamma \overline{p} / \overline{\rho} )$ represents the base flow and $(u'_i,p')$ the perturbations.
$\overline{c}^2 q_c$ is the acoustic source term.
\section{Usage}
\begin{lstlisting}[style=BashInputStyle]
APESolver session.xml
\end{lstlisting}
\section{Session file configuration}
\subsection{Solver Info}
\begin{itemize}
\item \inltt{Eqtype} Specifies the equation to solve. This should be set to
\inltt{APE}.
\item \inltt{UpwindType} Specifies the numerical interface flux scheme.
Currently, only \inltt{APEUpwind} supported.
\end{itemize}
\subsection{Parameters}
\begin{itemize}
\item \inltt{Gamma}: Ratio of specific heats
\end{itemize}
\subsection{Functions}
\begin{itemize}
\item \inltt{BaseFlow} Baseflow $(\overline{u}_i, \overline{p}, \overline{\rho})$ defined by the variables \inltt{u0,v0,w0,p0,rho0}
\item \inltt{InitialConditions}
\end{itemize}
\subsection{Boundary Conditions}
In addition to plain Dirichlet and Neumann boundary conditions, the APESolver features a wall boundray condition, a non-reflecting boundary and a white noise boundary condition.
\begin{itemize}
\item Rigid (Slip-) Wall Boundary Condition
\begin{lstlisting}[style=XmlStyle]
<BOUNDARYCONDITIONS>
<REGION REF="0">
<D VAR="p" USERDEFINEDTYPE="Wall" VALUE="0" />
<D VAR="u" USERDEFINEDTYPE="Wall" VALUE="0" />
<D VAR="v" USERDEFINEDTYPE="Wall" VALUE="0" />
<D VAR="w" USERDEFINEDTYPE="Wall" VALUE="0" />
</REGION>
</BOUNDARYCONDITIONS>
\end{lstlisting}
This BC imposes zero wall-normal perturbation velocity in a way that is more robust than using a Dirichlet boundary condition directly.
\item Non-Reflecting Boundary Condition
\begin{lstlisting}[style=XmlStyle]
<BOUNDARYCONDITIONS>
<REGION REF="0">
<D VAR="p" USERDEFINEDTYPE="RiemannInvariantBC"/>
<D VAR="u" USERDEFINEDTYPE="RiemannInvariantBC"/>
<D VAR="v" USERDEFINEDTYPE="RiemannInvariantBC"/>
<D VAR="w" USERDEFINEDTYPE="RiemannInvariantBC"/>
</REGION>
</BOUNDARYCONDITIONS>
\end{lstlisting}
The Riemann-Invariant BC approximates a non-reflecting (r.g. Farfield) boundary condition by setting incoming invariants to zero.
\item White Noise Boundary Condition
\begin{lstlisting}[style=XmlStyle]
<BOUNDARYCONDITIONS>
<REGION REF="0">
<D VAR="p" USERDEFINEDTYPE="Wall" VALUE="10" />
<D VAR="u" USERDEFINEDTYPE="Wall" VALUE="10" />
<D VAR="v" USERDEFINEDTYPE="Wall" VALUE="10" />
<D VAR="w" USERDEFINEDTYPE="Wall" VALUE="10" />
</REGION>
</BOUNDARYCONDITIONS>
\end{lstlisting}
The white noise BC imposes a stochastic, uniform pressure at the boundary. The implementation uses a Mersenne-Twister pseudo random number generator to generate white Gaussian noise.
The standard deviation $\sigma$ of the pressure is specified by the \inltt{VALUE} attribute.
\end{itemize}
\section{Examples}
\subsection{Aeroacoustic Wave Propagation}
In this section we explain how to set up a simple simulation of aeroacoustics in
Nektar++. We will study the propagation of an acoustic wave in the simple case
where the base flow is $\overline{u}_i = 0, \, \overline{p}=p_{\infty}=10^6, \, \overline{\rho} = \rho_0 = 1.204$. The geometry consists
of $64$ quadrilateral elements, as shown in Fig.~\ref{f:apesolver:geometry}.
\begin{figure}
\centering
\includegraphics[width=0.5\linewidth]{img/APE_Geometry.png}
\caption{Geometry used for the example case of modelling propagation of
acoustic waves where $\overline{u}_i = 0, \, \overline{p}=p_{\infty}=10^6, \, \overline{\rho} = \rho_0 = 1.204$}
\label{f:apesolver:geometry}
\end{figure}
\subsubsection{Input file}
We require a discontinuous Galerkin projection and use an explicit
fourth-order Runge-Kutta time integration scheme. We therefore set the following
solver information:
\begin{lstlisting}[style=XmlStyle]
<I PROPERTY="EQType" VALUE="APE"/>
<I PROPERTY="Projection" VALUE="DisContinuous"/>
<I PROPERTY="TimeIntegrationMethod" VALUE="ClassicalRungeKutta4"/>
<I PROPERTY="UpwindType" VALUE="APEUpwind"/>
\end{lstlisting}
To maintain numerical stability we must use a small time-step. The total
simulation time is $150$ time units. Finally, we set the density, heat ratio and
ambient pressure.
\begin{lstlisting}[style=XMLStyle]
<P> TimeStep = 0.00001 </P>
<P> NumSteps = 150 </P>
<P> FinTime = TimeStep*NumSteps </P>
<P> Rho0 = 1.204 </P> <!-- Incompressible density -->
<P> Gamma = 1.4 </P> <!-- Ratio of specific heats -->
<P> Pinfinity = 100000 </P> <!-- Ambient pressure -->
\end{lstlisting}
Let us note that to solve efficiently this problem a discontinuous Garlerkin
approach was used. The
system is excited via the initial conditions putting a Gaussian pulse for pulse
fluctuations. Finally, it is necessary to specify the base flow and the
eventual source terms using the following functions:
\begin{lstlisting}[style=XmlStyle]
<FUNCTION NAME="Baseflow">
<E VAR="u0" VALUE="0" />
<E VAR="v0" VALUE="0" />
<E VAR="p0" VALUE="Pinfinity" />
<E VAR="rho0" VALUE="Rho0" />
</FUNCTION>
<!-- Gaussian pulse located at the origin -->
<FUNCTION NAME="InitialConditions">
<E VAR="p" VALUE="100*exp(-32*(((x)*(x))+((y)*(y))))" />
<E VAR="u" VALUE="0" />
<E VAR="v" VALUE="0" />
</FUNCTION>
\end{lstlisting}
\subsubsection{Running the code}
\begin{lstlisting}[style=BashInputStyle]
APESolver Test_pulse.xml
\end{lstlisting}
\subsubsection{Results}
Fig.~\ref{f:apesolver:results} shows the pressure profile at different
time steps, showing the acoustic propagation.
\begin{figure}
\centering
\includegraphics[width=0.3\linewidth]{img/Prop_1.png}
\includegraphics[width=0.3\linewidth]{img/Prop_2.png}
\includegraphics[width=0.3\linewidth]{img/Prop_3.png}
\caption{}
\label{f:apesolver:results}
\end{figure}
It is possible to show the profile of the pressure perturbations with respect to
the spatial coordinate. The pressure fluctuations, that are concentrated in a
specific locations at the beginning (as specified by the initial conditions),
propagate with time and for sufficiently large time the decay is exponential as
predicted by literature \cite{DoFf83}.
\begin{figure}
\centering
\includegraphics[width=0.7\linewidth]{img/prog_4.png}
\caption{}
\end{figure}
%\chapter{Solvers}
%\label{s:solvers}
\input{acoustic}
\input{adr}
\input{ape}
\input{cardiac-ep}
\input{compressible-flow}
\input{dummy}
......
......@@ -362,7 +362,7 @@ To count the number of DOF in a solution file,
one can run the following command
%
\begin{lstlisting}[style=BashInputStyle]
FieldConvert file.xml file.fld out.stdout
FieldConvert -m dof file.xml file.fld out.stdout
\end{lstlisting}
%
%
......
......@@ -227,7 +227,7 @@ where it exists. The table below indicates support currently implemented.
\toprule
\textbf{Format} & \textbf{Extension} & \textbf{High-order} & \textbf{Notes}\\
\midrule
Gmsh & \texttt{msh} & \cmark & Only reads nodes, elements and physical groups (which are mapped to composites).\\
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.. \\
......
......@@ -82,7 +82,7 @@ The coupling waits for the file given in the receive function to appear.
The Cwipi coupling is only available when Nektar++ is compiled with OpenMPI and CWIPI
\end{notebox}
The Cwipi coupling uses CWIPI\footnote{http://sites.onera.fr/cwipi/} to facilitate real time data exchange over MPI.
See \cite{LaSaJa17} for details.
See \cite{La18} for details.
All data transfers are non-blocking to minimize the computational overhead.
The interface must be enabled with the command line option \inltt{--cwipi} and a unique application name, e.g:\begin{lstlisting}[style=BashInputStyle]
DummySolver --cwipi 'Dummy1' Dummy_3DCubeCwipi_1.xml
......
......@@ -129,7 +129,7 @@ void ProcessInterpPointDataToFld::Process(po::variables_map &vm)
MemoryManager<LibUtilities::PtsField>::AllocateSharedPtr(3, intFields);
int coord_id = m_config["interpcoord"].as<int>();
ASSERTL0(coord_id <= fieldPts->GetDim() - 1,
ASSERTL0(coord_id <= outPts->GetDim() - 1,
"interpcoord is bigger than the Pts files dimension");
Interpolator interp(LibUtilities::eNoMethod, coord_id);
......
......@@ -95,10 +95,6 @@ void ProcessWSS::Process(po::variables_map &vm)
"be computed");
}
// Load viscosity coefficients
NekDouble kinvis, lambda;
GetViscosity( kinvis, lambda);
// Declare arrays
int nshear = m_spacedim + 1;
int nstress = m_spacedim * m_spacedim;
......@@ -183,6 +179,11 @@ void ProcessWSS::Process(po::variables_map &vm)
// Extract Velocities
GetVelocity( BndElmtExp, velocity);
// Extract viscosity coefficients
NekDouble lambda;
Array<OneD, NekDouble> mu(nqe, 0.0);
GetViscosity( BndElmtExp, mu, lambda);
// Compute gradients
for (i = 0; i < m_spacedim; ++i)
{
......@@ -204,7 +205,8 @@ void ProcessWSS::Process(po::variables_map &vm)
}
// Velocity divergence scaled by lambda * mu
Vmath::Smul(nqe, lambda*kinvis, div, 1, div, 1);
Vmath::Smul(nqe, lambda, div, 1, div, 1);
Vmath::Vmul(nqe, mu, 1, div, 1, div, 1);
// Compute stress component terms
// tau_ij = mu*(u_i,j + u_j,i) + mu*lambda*delta_ij*div(u)
......@@ -216,7 +218,7 @@ void ProcessWSS::Process(po::variables_map &vm)
grad[j * m_spacedim + i], 1,
stress[i * m_spacedim + j], 1);
Vmath::Smul(nqe, kinvis,
Vmath::Vmul(nqe, mu, 1,
stress[i * m_spacedim + j], 1,
stress[i * m_spacedim + j], 1);
......@@ -294,20 +296,93 @@ void ProcessWSS::Process(po::variables_map &vm)
}
}
void ProcessWSS::GetViscosity(NekDouble &kinvis, NekDouble &lambda)
void ProcessWSS::GetViscosity(
const Array<OneD, MultiRegions::ExpListSharedPtr> exp,
Array<OneD, NekDouble> &mu,
NekDouble &lambda)
{
NekDouble m_mu;
int npoints = exp[0]->GetNpoints();
if(boost::iequals(m_f->m_variables[0], "u"))
{
// IncNavierStokesSolver
kinvis = m_f->m_session->GetParameter("Kinvis");
m_mu = m_f->m_session->GetParameter("Kinvis");
Vmath::Fill(npoints, m_mu, mu, 1);
lambda = 0;
}
else if(boost::iequals(m_f->m_variables[0], "rho") &&
boost::iequals(m_f->m_variables[1], "rhou"))
{
// CompressibleFlowSolver
m_f->m_session->LoadParameter ("mu", kinvis, 1.78e-05);
std::string m_ViscosityType;
m_f->m_session->LoadParameter ("mu", m_mu, 1.78e-05);
m_f->m_session->LoadParameter ("lambda", lambda, -2.0/3.0);
m_f->m_session->LoadSolverInfo("ViscosityType", m_ViscosityType
, "Constant");
if (m_ViscosityType == "Variable")
{
// Check equation of state
std::string eosType;
bool m_idealGas;
m_f->m_session->LoadSolverInfo("EquationOfState", eosType,
"IdealGas");
m_idealGas = boost::iequals(eosType,"IdealGas");
ASSERTL0(m_idealGas,
"Only IdealGas EOS implemented for Variable ViscosityType");
// Get relevant parameters
NekDouble m_gamma;
NekDouble m_pInf;
NekDouble m_rhoInf;
NekDouble m_gasConstant;
NekDouble cv_inv;
m_f->m_session->LoadParameter("Gamma", m_gamma, 1.4);
m_f->m_session->LoadParameter("pInf", m_pInf, 101325);
m_f->m_session->LoadParameter("rhoInf", m_rhoInf, 1.225);
m_f->m_session->LoadParameter("GasConstant", m_gasConstant
, 287.058);
// Get temperature from flowfield
cv_inv = (m_gamma - 1.0) / m_gasConstant;
// e = 1/rho ( E - 1/2 ( rhou^2/rho + ... ) )
Array<OneD, NekDouble> tmp(npoints, 0.0);
Array<OneD, NekDouble> energy(npoints, 0.0);
Array<OneD, NekDouble> temperature(npoints, 0.0);
Vmath::Vcopy(npoints, exp[m_spacedim+1]->GetPhys(), 1, energy, 1);
for (int i = 0; i < m_spacedim; i++)
{
// rhou^2
Vmath::Vmul(npoints, exp[i + 1]->GetPhys(), 1
, exp[i + 1]->GetPhys(), 1, tmp, 1);
// rhou^2/rho
Vmath::Vdiv(npoints, tmp, 1, exp[0]->GetPhys(), 1, tmp, 1);
// 0.5 rhou^2/rho
Vmath::Smul(npoints, 0.5, tmp, 1, tmp, 1);
// E - 0.5 rhou^2/rho - ...
Vmath::Vsub(npoints, energy, 1, tmp, 1, energy, 1);
}
// rhoe/rho
Vmath::Vdiv(npoints, energy, 1, exp[0]->GetPhys(), 1, energy, 1);
// T = e/Cv
Vmath::Smul(npoints, cv_inv, energy, 1, temperature, 1 );
// Variable viscosity through the Sutherland's law
NekDouble mu_star = m_mu;
NekDouble T_star = m_pInf / (m_rhoInf * m_gasConstant);
NekDouble ratio;
for (int i = 0; i < npoints; ++i)
{
ratio = temperature[i] / T_star;
mu[i] = mu_star * ratio * sqrt(ratio) *
(T_star + 110.0) / (temperature[i] + 110.0);
}
}
else
{
Vmath::Fill(npoints, m_mu, mu, 1);
}
}
else
{
......
......@@ -74,7 +74,8 @@ public:
}
protected:
void GetViscosity(NekDouble &kinvis, NekDouble &lambda);
void GetViscosity(const Array<OneD, MultiRegions::ExpListSharedPtr> exp,
Array<OneD, NekDouble > &mu, NekDouble &lambda);
void GetVelocity(const Array<OneD, MultiRegions::ExpListSharedPtr> exp,
Array<OneD, Array<OneD, NekDouble> > &vel);
......
......@@ -1063,6 +1063,7 @@ void FieldIOHdf5::v_Import(const std::string &infilename,
}
}
ImportHDF5FieldMetaData(dataSource, fieldinfomap);
m_comm->Block();
}
......
......@@ -123,11 +123,6 @@ void Interpolator::CalcWeights(const LibUtilities::PtsFieldSharedPtr ptsInField,
m_weights = Array<TwoD, float>(nOutPts, 3, 0.0);
m_neighInds = Array<TwoD, unsigned int>(nOutPts, 3, (unsigned int) 0);
if (m_ptsInField->GetDim() == 1)
{
m_coordId = 0;
}
for (int i = 0; i < nOutPts; ++i)
{
Array<OneD, NekDouble> tmp(m_dim, 0.0);
......
......@@ -378,7 +378,7 @@ TARGET_LINK_LIBRARIES(LibUtilities LINK_PUBLIC
${Boost_PROGRAM_OPTIONS_LIBRARY}
${Boost_FILESYSTEM_LIBRARY}
${Boost_SYSTEM_LIBRARY}
${ZLIB_LIBRARIES}
optimized ${ZLIB_LIBRARIES} debug ${ZLIB_LIBRARIES_DEBUG}
)
# TinyXML
......
......@@ -306,6 +306,170 @@ namespace Nektar
}
}
void Expansion::ComputeGmatcdotMF(
const Array<TwoD, const NekDouble> &df,
const Array<OneD, const NekDouble> &direction,
Array<OneD, Array<OneD, NekDouble> > &dfdir)
{
int shapedim = dfdir.num_elements();
int coordim = GetCoordim();
int nqtot = direction.num_elements()/coordim;
for(int j = 0; j < shapedim; j++)
{
dfdir[j] = Array<OneD, NekDouble>(nqtot,0.0);
for (int k = 0; k < coordim; k++)
{
if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
{
Vmath::Vvtvp(nqtot,
&df[shapedim*k+j][0], 1,
&direction[k*nqtot], 1,
&dfdir[j][0], 1,
&dfdir[j][0], 1);
}
else
{
Vmath::Svtvp(nqtot,
df[shapedim*k+j][0],
&direction[k*nqtot], 1,
&dfdir[j][0], 1,