Commit 71f69611 authored by Chris Cantwell's avatar Chris Cantwell
Browse files

Merge branch 'master' into feature/NekDouble

parents 17a74f80 0d2c94c5
......@@ -9,11 +9,9 @@ library/Demos/StdRegions/ExtraDemos
library/Demos/MultiRegions/ExtraDemos
solvers/ADR2DManifoldSolver
solvers/FitzHughNagumoSolver
solvers/CardiacEPSolver
solvers/FluxReconstruction
solvers/VortexWaveInteraction
solvers/ImageWarpingSolver
solvers/PulseWaveSolver
solvers/VortexWaveInteraction
docs/*.doc
docs/arch
docs/emacs
......
......@@ -186,27 +186,21 @@ int main(int argc, char *argv[])
// /////////////////////////////////////////////////////////
// /////////////////////////////////////////////////////////
// Interpolation //
// /////////////////////////////////////////////////////////
// /////////////////////////////////////////////////////////
// Generate a list of interpolating nodes
int nNodes = 2*nPts; // Number of interpolating nodes
boost::shared_ptr<Points<NekDouble> > nodes = PointsManager()[PointsKey(nNodes, pointsType)];
// const ptr<Points<NekDouble> > nodes = PointsManager()[PointsKey(nNodes, pointsType)];
//SharedArray<const NekDouble> zNode = nodes->GetZ();
Array<OneD, const NekDouble> zNode = nodes->GetZ();
// Get the interpolation matrix I
// Note that I is 2N rows by N columns
const Points<NekDouble>::MatrixSharedPtrType Iptr = points->GetI(nNodes,zNode);
// const Points<NekDouble>::MatrixSharedPtrType Iptr = points->GetI(nNodes, zNode);
const Points<NekDouble>::MatrixSharedPtrType Iptr = points->GetI(nNodes,zNode);
const NekMatrix<NekDouble> & I = *Iptr;
// Interpolate the data values in the y vector using the interpolation matrix I
// SharedArray<NekDouble> u(I.GetRows());
Array<OneD, NekDouble> u(I.GetRows());
for(int i = 0; i < int(I.GetRows()); ++i)
{
......@@ -216,7 +210,7 @@ int main(int argc, char *argv[])
u[i] += I(i,j) * y[j];
}
}
// Display the original samples
cout << setprecision(3);
cout << "\nOriginal data: \nx = ";
......@@ -282,7 +276,92 @@ int main(int argc, char *argv[])
}
///////////////////////////////////////////////////////////
// Galkerin Projection //
///////////////////////////////////////////////////////////
// Generate a list of projection nodes
Array<OneD, NekDouble> yNode(zNode.num_elements());
for(int i = 0; i < nNodes; ++i)
{
yNode[i] = func( zNode[i], nNodes, pointsType );
}
PointsKey key1(nNodes, pointsType);
// Note that I is 2N rows by N columns
const Points<NekDouble>::MatrixSharedPtrType GPptr = points->GetGalerkinProjection(key1);
const NekMatrix<NekDouble> & GP = *GPptr;
// Project the data values in the yNode vector using the projection matrix GP
for(int i = 0; i < int(GP.GetRows()); ++i)
{
u[i] = 0;
for(int j = 0; j < int(GP.GetColumns()); ++j)
{
u[i] += GP(i,j) * yNode[j];
}
}
// Display the original samples
cout << setprecision(3);
cout << "\n\n\n **** Galerkin Project ****";
cout << "\n\nResults of Galkerin Project with " << kPointsTypeStr[pointsType] << ":";
cout << "\nOriginal data: \nx = ";
for(int i = 0; i < nPts; ++i)
{
cout << setw(6) << z[i] << " ";
}
cout << "\ny = ";
for(int i = 0; i < nPts; ++i)
{
cout << setw(6) << y[i] << " ";
}
cout << "\nproject = ";
for(int i = 0; i < nPts; ++i)
{
cout << setw(6) << u[i] << " ";
}
// Display the pointwise error
cout << setprecision(1);
cout << "\nerror = ";
Linf = 0;
RMS = 0;
for(int i = 0; i < int(GP.GetRows()); ++i)
{
long double exact = func(z[i], nPts, pointsType);
long double error = exact - u[i];
Linf = max(Linf, fabs(error));
RMS += error*error;
long double epsilon = 1e-2;
if( fabs(exact) > epsilon )
{
error /= exact;
}
cout << setw(6) << error << " ";
}
RMS = sqrt(RMS) / int(GP.GetRows());
cout << setprecision(6);
cout << "\nLinf = " << setw(6) << Linf;
cout << "\nRMS = " << setw(6) << RMS << endl;
// Show the projection matrix
cout << "\nI = " << endl;
for(int i = 0; i < int(GP.GetRows()); ++i)
{
cout << " ";
for(int j = 0; j < int(GP.GetColumns()); ++j)
{
printf("% 5.3f ", GP(i,j));
}
cout << "" << endl;
}
// /////////////////////////////////////////////////////////
......
......@@ -8,10 +8,10 @@
</files>
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-12">0.000418001</value>
<value tolerance="1e-12">0.000416575</value>
</metric>
<metric type="Linf" id="2">
<value tolerance="1e-12">0.000874515</value>
<value tolerance="1e-12">0.000871589</value>
</metric>
</metrics>
</test>
......
......@@ -8,10 +8,10 @@
</files>
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-8">0.000418001</value>
<value tolerance="1e-8">0.000416575</value>
</metric>
<metric type="Linf" id="2">
<value tolerance="1e-8">0.000874515</value>
<value tolerance="1e-8">0.000871589</value>
</metric>
</metrics>
</test>
......
......@@ -9,10 +9,10 @@
</files>
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-12">0.000418001</value>
<value tolerance="1e-12">0.000416575</value>
</metric>
<metric type="Linf" id="2">
<value tolerance="1e-12">0.000874515</value>
<value tolerance="1e-12">0.000871589</value>
</metric>
</metrics>
</test>
......
......@@ -9,10 +9,10 @@
</files>
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-12">0.000418001</value>
<value tolerance="1e-12">0.000416575</value>
</metric>
<metric type="Linf" id="2">
<value tolerance="1e-12">0.000874515</value>
<value tolerance="1e-12">0.000871589</value>
</metric>
</metrics>
</test>
......
......@@ -8,10 +8,10 @@
</files>
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-12">0.00042687</value>
<value tolerance="1e-12">0.000417674</value>
</metric>
<metric type="Linf" id="2">
<value tolerance="1e-12">0.00170172</value>
<value tolerance="1e-12">0.00168661</value>
</metric>
</metrics>
</test>
......@@ -9,10 +9,10 @@
</files>
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-12">0.00042687</value>
<value tolerance="1e-12">0.000417674</value>
</metric>
<metric type="Linf" id="2">
<value tolerance="1e-12">0.00170173</value>
<value tolerance="1e-06">0.00168661</value>
</metric>
</metrics>
</test>
......@@ -392,10 +392,10 @@
<N VAR="u" VALUE="-sin(PI/2*x)*cos(PI/2*y)*sin(PI/2*z)*PI/2" />
</REGION>
<REGION REF="2">
<R VAR="u" VALUE="cos(PI/2*x)*sin(PI/2*y)*sin(PI/2*z)*PI/2+sin(PI/2*x)*sin(PI/2*y)*sin(PI/2*z)" PRIMCOEFF="1" />
<R VAR="u" VALUE="cos(PI/2*x)*sin(PI/2*y)*sin(PI/2*z)*PI/2+2*sin(PI/2*x)*sin(PI/2*y)*sin(PI/2*z)" PRIMCOEFF="2" />
</REGION>
<REGION REF="3">
<R VAR="u" VALUE="sin(PI/2*x)*cos(PI/2*y)*sin(PI/2*z)*PI/2+sin(PI/2*x)*sin(PI/2*y)*sin(PI/2*z)" PRIMCOEFF="1" />
<R VAR="u" VALUE="sin(PI/2*x)*cos(PI/2*y)*sin(PI/2*z)*PI/2+3*sin(PI/2*x)*sin(PI/2*y)*sin(PI/2*z)" PRIMCOEFF="3" />
</REGION>
<REGION REF="4">
<D VAR="u" VALUE="sin(PI/2*x)*sin(PI/2*y)*sin(PI/2*z)" />
......
......@@ -31,6 +31,7 @@ int main(int argc, char *argv[])
fprintf(stderr,"\t Chebyshev = 10\n");
fprintf(stderr,"\t Monomial = 11\n");
fprintf(stderr,"\t FourierSingleMode = 12\n");
fprintf(stderr,"\t Gauss Lagrange = 15\n");
fprintf(stderr,"Note type = 1,2,4,5 are for higher dimensional basis\n");
......@@ -81,7 +82,7 @@ int main(int argc, char *argv[])
}
else
{
Qtype = LibUtilities::eGaussLobattoLegendre;
Qtype = LibUtilities::eGaussLobattoLegendre;
}
......
......@@ -44,12 +44,12 @@
namespace ErrorUtil
{
static boost::optional<std::ostream&> outStream;
static void SetErrorStream(std::ostream& o)
{
outStream = o;
}
static bool HasCustomErrorStream()
{
return outStream;
......@@ -63,8 +63,8 @@ namespace ErrorUtil
class NekError : public std::runtime_error
{
public:
NekError(const std::string& message) : std::runtime_error(message) {}
public:
NekError(const std::string& message) : std::runtime_error(message) {}
};
static void Error(ErrType type, const char *routine, int lineNumber, const char *msg, unsigned int level)
......
......@@ -69,6 +69,7 @@ SET(FoundationHeaders
./Foundations/NodalTriEvenlySpaced.h
./Foundations/NodalTetEvenlySpaced.h
./Foundations/NodalPrismEvenlySpaced.h
./Foundations/PhysGalerkinProject.h
)
SET(InterpreterHeaders
......@@ -167,6 +168,7 @@ SET(FoundationSources
./Foundations/NodalTriEvenlySpaced.cpp
./Foundations/NodalTetEvenlySpaced.cpp
./Foundations/NodalPrismEvenlySpaced.cpp
./Foundations/PhysGalerkinProject.cpp
)
SET(InterpreterSources
......
......@@ -654,6 +654,23 @@ namespace Nektar
}
}
void Transposition::SetSpecVanVisc(Array<OneD, NekDouble> visc)
{
m_specVanVisc = visc;
}
NekDouble Transposition::GetSpecVanVisc(const int k)
{
NekDouble returnval = 0.0;
if(m_specVanVisc.num_elements())
{
returnval = m_specVanVisc[k];
}
return returnval;
}
// TODO: Impelement 2D and 3D transposition routines
}
......
......@@ -94,6 +94,10 @@ namespace Nektar
bool UseNumMode = false,
TranspositionDir dir = eNoTrans);
LIB_UTILITIES_EXPORT void SetSpecVanVisc(Array<OneD, NekDouble> visc);
LIB_UTILITIES_EXPORT NekDouble GetSpecVanVisc(const int k);
protected:
CommSharedPtr m_hcomm;
......@@ -156,6 +160,9 @@ namespace Nektar
/// MPI_Alltoallv offset map of send/recv buffer in global vector.
Array<OneD,int> m_OffsetMap;
/// Spectral vanishing Viscosity coefficient for stabilisation
Array<OneD, NekDouble> m_specVanVisc;
};
typedef boost::shared_ptr<Transposition> TranspositionSharedPtr;
......
......@@ -307,6 +307,7 @@ namespace Nektar
return GetI(numpoints, xpoints);
}
const boost::shared_ptr<NekMatrix<NekDouble> > GaussPoints::GetI(const PointsKey &pkey)
{
ASSERTL0(pkey.GetPointsDim()==1, "Gauss Points can only interp to other 1d point distributions");
......@@ -391,6 +392,49 @@ namespace Nektar
}
return y;
}
const boost::shared_ptr<NekMatrix<NekDouble> > GaussPoints::GetGalerkinProjection(const PointsKey &pkey)
{
return m_GalerkinProjectionManager[pkey];
}
boost::shared_ptr< NekMatrix<NekDouble> > GaussPoints::CreateGPMatrix(const PointsKey &pkey)
{
boost::shared_ptr< NekMatrix<NekDouble> > returnval = CalculateGalerkinProjectionMatrix(pkey);
// Delegate to function below
return returnval;
}
boost::shared_ptr<NekMatrix<NekDouble> > GaussPoints::CalculateGalerkinProjectionMatrix(const PointsKey &pkey)
{
int numpointsfrom = pkey.GetNumPoints();
int numpointsto = GetNumPoints();
Array<OneD, const NekDouble> weightsfrom;
weightsfrom = PointsManager()[pkey]->GetW();
boost::shared_ptr< NekMatrix<NekDouble> > Interp = GetI(pkey);
Array<OneD, NekDouble> GalProj(numpointsfrom*numpointsto);
// set up inner product matrix and multiply by inverse of
// diagaonal mass matrix
for(int i = 0; i < numpointsto; ++i)
{
Vmath::Vmul(numpointsfrom,Interp->GetPtr().get() +i*numpointsfrom,1,
&weightsfrom[0],1,&GalProj[0] +i,numpointsto);
Vmath::Smul(numpointsfrom,1.0/m_weights[i],&GalProj[0]+i,numpointsto,
&GalProj[0]+i,numpointsto);
}
NekDouble* t = GalProj.data();
boost::shared_ptr< NekMatrix<NekDouble> > returnval(MemoryManager<NekMatrix<NekDouble> >::AllocateSharedPtr(numpointsto,numpointsfrom,t));
return returnval;
}
} // end of namespace LibUtilities
} // end of namespace Nektar
......
......@@ -62,6 +62,12 @@ namespace Nektar
LIB_UTILITIES_EXPORT const boost::shared_ptr<NekMatrix<NekDouble> > GetI(const Array<OneD, const NekDouble>& x);
LIB_UTILITIES_EXPORT const boost::shared_ptr<NekMatrix<NekDouble> > GetI(unsigned int numpoints, const Array<OneD, const NekDouble>& x);
LIB_UTILITIES_EXPORT boost::shared_ptr< NekMatrix<NekDouble> > CreateGPMatrix(const PointsKey &pkey);
LIB_UTILITIES_EXPORT const boost::shared_ptr<NekMatrix<NekDouble> > GetGalerkinProjection(const PointsKey &pkey);
GaussPoints(const PointsKey &pkey):PointsBaseType(pkey)
{
m_InterpManager.RegisterCreator(PointsKey(0, eGaussGaussLegendre),
......@@ -101,7 +107,46 @@ namespace Nektar
m_InterpManager.RegisterCreator(PointsKey(0, ePolyEvenlySpaced),
boost::bind(&GaussPoints::CreateMatrix, this, _1));
m_InterpManager.RegisterCreator(PointsKey(0, eBoundaryLayerPoints),
boost::bind(&GaussPoints::CreateMatrix, this, _1)); }
boost::bind(&GaussPoints::CreateMatrix, this, _1));
m_GalerkinProjectionManager.RegisterCreator(PointsKey(0, eGaussGaussLegendre),
boost::bind(&GaussPoints::CreateGPMatrix, this, _1));
m_GalerkinProjectionManager.RegisterCreator(PointsKey(0, eGaussRadauMLegendre),
boost::bind(&GaussPoints::CreateGPMatrix, this, _1));
m_GalerkinProjectionManager.RegisterCreator(PointsKey(0, eGaussRadauPLegendre),
boost::bind(&GaussPoints::CreateGPMatrix, this, _1));
m_GalerkinProjectionManager.RegisterCreator(PointsKey(0, eGaussLobattoLegendre),
boost::bind(&GaussPoints::CreateGPMatrix, this, _1));
m_GalerkinProjectionManager.RegisterCreator(PointsKey(0, eGaussGaussChebyshev),
boost::bind(&GaussPoints::CreateGPMatrix, this, _1));
m_GalerkinProjectionManager.RegisterCreator(PointsKey(0, eGaussRadauMChebyshev),
boost::bind(&GaussPoints::CreateGPMatrix, this, _1));
m_GalerkinProjectionManager.RegisterCreator(PointsKey(0, eGaussRadauPChebyshev),
boost::bind(&GaussPoints::CreateGPMatrix, this, _1));
m_GalerkinProjectionManager.RegisterCreator(PointsKey(0, eGaussLobattoChebyshev),
boost::bind(&GaussPoints::CreateGPMatrix, this, _1));
m_GalerkinProjectionManager.RegisterCreator(PointsKey(0, eGaussRadauMAlpha0Beta1),
boost::bind(&GaussPoints::CreateGPMatrix, this, _1));
m_GalerkinProjectionManager.RegisterCreator(PointsKey(0, eGaussRadauMAlpha0Beta2),
boost::bind(&GaussPoints::CreateGPMatrix, this, _1));
m_GalerkinProjectionManager.RegisterCreator(PointsKey(0, eGaussRadauMAlpha1Beta0),
boost::bind(&GaussPoints::CreateGPMatrix, this, _1));
m_GalerkinProjectionManager.RegisterCreator(PointsKey(0, eGaussRadauMAlpha2Beta0),
boost::bind(&GaussPoints::CreateGPMatrix, this, _1));
m_GalerkinProjectionManager.RegisterCreator(PointsKey(0, eGaussKronrodLegendre),
boost::bind(&GaussPoints::CreateGPMatrix, this, _1));
m_GalerkinProjectionManager.RegisterCreator(PointsKey(0, eGaussRadauKronrodMLegendre),
boost::bind(&GaussPoints::CreateGPMatrix, this, _1));
m_GalerkinProjectionManager.RegisterCreator(PointsKey(0, eGaussRadauKronrodMAlpha1Beta0),
boost::bind(&GaussPoints::CreateGPMatrix, this, _1));
m_GalerkinProjectionManager.RegisterCreator(PointsKey(0, eGaussLobattoKronrodLegendre),
boost::bind(&GaussPoints::CreateGPMatrix, this, _1));
m_GalerkinProjectionManager.RegisterCreator(PointsKey(0, eFourierEvenlySpaced),
boost::bind(&GaussPoints::CreateGPMatrix, this, _1));
m_GalerkinProjectionManager.RegisterCreator(PointsKey(0, ePolyEvenlySpaced),
boost::bind(&GaussPoints::CreateGPMatrix, this, _1));
m_GalerkinProjectionManager.RegisterCreator(PointsKey(0, eBoundaryLayerPoints),
boost::bind(&GaussPoints::CreateGPMatrix, this, _1)); }
private:
/// These should not be called. All creation is done
......@@ -115,6 +160,10 @@ namespace Nektar
void CalculateDerivMatrix();
void CalculateInterpMatrix(unsigned int npts, const Array<OneD, const NekDouble>& xpoints, Array<OneD, NekDouble>& interp);
boost::shared_ptr<NekMatrix<NekDouble> > CalculateGalerkinProjectionMatrix(const PointsKey &pkey);
/// functions used by the Kronrod points
NekDouble LagrangeInterpolant(NekDouble x, int npts, const Array<OneD, const NekDouble>& xpts, const Array<OneD, const NekDouble>& funcvals);
NekDouble LagrangePoly(NekDouble x, int pt, int npts, const Array<OneD, const NekDouble>& xpts);
......
///////////////////////////////////////////////////////////////////////////////
//
// File PhysGalerkinProject.cpp
//
// For more information, please see: http://www.nektar.info
//
// The MIT License
//
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// License for the specific language governing rights and limitations under
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: Definition of Physical Space Galerkin Projection methods
//
///////////////////////////////////////////////////////////////////////////////
#include <LibUtilities/Foundations/PhysGalerkinProject.h>
#include <LibUtilities/BasicUtils/Vmath.hpp>
#include <LibUtilities/BasicUtils/VmathArray.hpp>
#include <LibUtilities/Foundations/ManagerAccess.h>
#include <LibUtilities/LinearAlgebra/NekTypeDefs.hpp>
#include <LibUtilities/Foundations/Points.h>
#include <LibUtilities/Foundations/Basis.h>
namespace Nektar
{
namespace LibUtilities
{
// Physical Space Interpolation methods
// 1D Interpolation
void PhysGalerkinProject1D(const BasisKey &fbasis0,
const Array<OneD, const NekDouble>& from,
const BasisKey &tbasis0,
Array<OneD, NekDouble> &to)
{
PhysGalerkinProject1D(fbasis0.GetPointsKey(),from,tbasis0.GetPointsKey(),to);
}
void PhysGalerkinProject1D(const PointsKey &fpoints0,
const Array<OneD, const NekDouble>& from,
const PointsKey &tpoints0,
Array<OneD, NekDouble> &to)
{
if(fpoints0 == tpoints0) //check to see if the same
{
Vmath::Vcopy(fpoints0.GetNumPoints(),from,1,to,1);
}
else // interpolate
{
DNekMatSharedPtr GP0;
GP0 = PointsManager()[tpoints0]->GetGalerkinProjection(fpoints0);
NekVector<NekDouble> in(fpoints0.GetNumPoints(),from,eWrapper);
NekVector<NekDouble> out(tpoints0.GetNumPoints(),to,eWrapper);
GP0->Transpose();
out = (*GP0)*in;
}
}
void PhysGalerkinProject1D(const BasisKey &fbasis0,
const NekDouble *from,
const BasisKey &tbasis0,
NekDouble *to)
{
PhysGalerkinProject1D(fbasis0.GetPointsKey(),from,tbasis0.GetPointsKey(),to);
}
void PhysGalerkinProject1D(const PointsKey &fpoints0,
const NekDouble *from,
const PointsKey &tpoints0,
NekDouble *to)
{
if(fpoints0 == tpoints0) //check to see if the same
{
Vmath::Vcopy(fpoints0.GetNumPoints(),from,1,to,1);
}
else // interpolate
{
DNekMatSharedPtr GP0;