Commit 0492cc31 authored by Dave Moxey's avatar Dave Moxey

Merge branch 'feature/barycentricInterpolation' into 'master'

Add barycentric interpolation to PhysEvaluate

See merge request !1126
parents 0d4a83a8 29c69911
ADD_NEKTAR_EXECUTABLE(StdProject
COMPONENT demos DEPENDS StdRegions SOURCES StdProject.cpp)
ADD_NEKTAR_EXECUTABLE(StdInterp
COMPONENT demos DEPENDS StdRegions SOURCES StdInterp.cpp)
ADD_NEKTAR_EXECUTABLE(StdInterpBasis
COMPONENT demos DEPENDS StdRegions SOURCES StdInterpBasis.cpp)
ADD_NEKTAR_EXECUTABLE(StdEquiToCoeff2D
COMPONENT demos DEPENDS StdRegions SOURCES StdEquiToCoeff2D.cpp)
# Interpolation tests.
ADD_NEKTAR_TEST(StdInterp_Hex_Mod_P6_Q7)
ADD_NEKTAR_TEST(StdInterp_Quad_Lagrange_P6_Q7)
ADD_NEKTAR_TEST(StdInterp_Prism_Mod_P6_Q7)
ADD_NEKTAR_TEST(StdInterp_Pyr_Mod_P6_Q7)
ADD_NEKTAR_TEST(StdInterp_Tet_Mod_P6_Q7)
ADD_NEKTAR_TEST(StdInterp_Tri_Mod_P6_Q7)
# Basis interpolation tests.
ADD_NEKTAR_TEST(StdInterpBasis_Prism_Orth_P6_Q7)
ADD_NEKTAR_TEST(StdInterpBasis_Quad_Lagrange_P6_Q7)
ADD_NEKTAR_TEST(StdInterpBasis_Tet_Mod_P4_Q5)
ADD_NEKTAR_TEST(StdInterpBasis_Tri_Mod_P7_Q8)
ADD_NEKTAR_TEST(StdInterpBasis_Hex_Mod_P7_Q8)
ADD_NEKTAR_TEST(StdInterpBasis_Pyr_Mod_P7_Q8)
# Projection tests.
ADD_NEKTAR_TEST(StdProject1D_Seg_Orth_P6_Q7)
ADD_NEKTAR_TEST(StdProject1D_Seg_Mod_P6_Q7)
ADD_NEKTAR_TEST(StdProject1D_Single_Mode_Fourier_P2_Q2)
......
This diff is collapsed.
///////////////////////////////////////////////////////////////////////////////
//
// File: StdInterp.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).
//
// 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: Demo for testing functionality of PhysEvaluate
//
///////////////////////////////////////////////////////////////////////////////
#include <iostream>
#include <LibUtilities/BasicUtils/Timer.h>
#include "StdDemoSupport.hpp"
// Evaluate polynomial for testing and save in ret (size same as pts[0]) if
// tensorp = 0, we need tensorprod else just eval at pts
Array<OneD, NekDouble> EvalPoly(Array<OneD, Array<OneD, NekDouble>> &pts)
{
Array<OneD, NekDouble> ret(pts[0].size());
unsigned dim = pts.size();
// check if pts[0] and pts[1] have same size
// polynomial = x^2 + y^2 - 3x - 4
for (int i = 0; i < pts[0].size(); i++)
{
ret[i] = pow(pts[0][i],2) - 3*pts[0][i] - 4.0
+ (dim >= 2 ? pow(pts[1][i], 2) : 0.0)
+ (dim >= 3 ? pow(pts[2][i], 2) : 0.0);
}
return ret;
}
int main(int argc, char *argv[])
{
DemoSupport demo;
demo.ParseArguments(argc, argv);
StdExpansion *E = demo.CreateStdExpansion();
const auto totPoints = (unsigned) E->GetTotPoints();
const auto dimension = (unsigned) E->GetShapeDimension();
// Create a new element but with the evenly-spaced points type, so that we
// perform a PhysEvaluate at a different set of nodal points
// (i.e. non-collocated interpolation).
vector<string> &ptypes = demo.GetPointsType();
for (int i = 0; i < dimension; ++i)
{
ptypes[i] = "PolyEvenlySpaced";
}
StdExpansion *F = demo.CreateStdExpansion();
Array<OneD, Array<OneD, NekDouble>> coordsE = demo.GetCoords(E);
Array<OneD, Array<OneD, NekDouble>> coordsF = demo.GetCoords(F);
Array<OneD, NekDouble> physIn(totPoints), physOut(totPoints);
Array<OneD, NekDouble> tmpIn(dimension), sol(totPoints);
// Evaluate polynomial at the set of elemental solution points.
physIn = EvalPoly(coordsE);
for (int i = 0; i < totPoints; ++i)
{
for (int d = 0; d < dimension; ++d)
{
tmpIn[d] = coordsF[d][i];
}
physOut[i] = E->PhysEvaluate(tmpIn, physIn);
}
sol = EvalPoly(coordsF);
cout << "L infinity error : " << scientific << E->Linf(physOut, sol) << endl;
cout << "L 2 error : " << scientific << E->L2 (physOut, sol) << endl;
return 0;
}
///////////////////////////////////////////////////////////////////////////////
//
// File: StdInterpBasis.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).
//
// 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: Demo for testing functionality of PhysEvaluateBasis
//
///////////////////////////////////////////////////////////////////////////////
#include <iostream>
#include <LibUtilities/BasicUtils/Timer.h>
#include "StdDemoSupport.hpp"
int main(int argc, char *argv[])
{
DemoSupport demo;
demo.ParseArguments(argc, argv);
StdExpansion *E = demo.CreateStdExpansion();
if (E == nullptr)
{
return 1;
}
int nCoeffs = E->GetNcoeffs(), nPts = E->GetTotPoints();
int nTot = nCoeffs * nPts, dimension = E->GetShapeDimension();
Array<OneD, Array<OneD, NekDouble>> coords = demo.GetCoords(E);
Array<OneD, NekDouble> sol(nTot), phys(nTot), tmpIn(dimension);
// For each mode, we follow two approaches:
//
// 1) Evaluate the basis function at each quadrature point using the
// StdExpansion::PhysEvaluateBasis function.
// 2) Evaluate the basis function at all quadrature points using FillMode.
//
// These are then compared to ensure they give the same result.
for (int k = 0; k < nCoeffs; ++k)
{
// Evaluate each mode at the quadrature points.
for (int i = 0; i < nPts; ++i)
{
for (int d = 0; d < dimension; ++d)
{
tmpIn[d] = coords[d][i];
}
phys[k * nPts + i] = E->PhysEvaluateBasis(tmpIn, k);
}
// Fill the 'solution' field with each of the modes using FillMode.
Array<OneD, NekDouble> tmp = sol + k * nPts;
E->FillMode(k, tmp);
}
cout << "L infinity error : " << scientific << E->Linf(phys, sol) << endl;
cout << "L 2 error : " << scientific << E->L2(phys, sol) << endl;
return 0;
}
This diff is collapsed.
<?xml version="1.0" encoding="utf-8"?>
<test>
<description>StdInterpBasis Hexahedron Modified basis P=7 Q=8</description>
<executable>StdInterpBasis</executable>
<parameters>-s hexahedron -b Modified_A Modified_A Modified_A -o 7 7 7 -p 8 8 8</parameters>
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-12">0</value>
</metric>
<metric type="Linf" id="2">
<value tolerance="1e-12">0</value>
</metric>
</metrics>
</test>
<?xml version="1.0" encoding="utf-8"?>
<test>
<description>StdInterpBasis Prism Orthonormal basis P=6 Q=7</description>
<executable>StdInterpBasis</executable>
<parameters>-s prism -b Ortho_A Ortho_A Ortho_B -o 6 6 6 -p 7 7 7</parameters>
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-12">1.89694e-15</value>
</metric>
<metric type="Linf" id="2">
<value tolerance="1e-12">1.77636e-14</value>
</metric>
</metrics>
</test>
<?xml version="1.0" encoding="utf-8"?>
<test>
<description>StdInterpBasis Triangle Modified basis P=7 Q=8</description>
<executable>StdInterpBasis</executable>
<parameters>-s pyramid -b Modified_A Modified_A ModifiedPyr_C -o 7 7 7 -p 8 8 8</parameters>
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-12">0</value>
</metric>
<metric type="Linf" id="2">
<value tolerance="1e-12">0</value>
</metric>
</metrics>
</test>
<?xml version="1.0" encoding="utf-8"?>
<test>
<description>StdInterpBasis Quadrilateral Lagrange basis P=6 Q=7</description>
<executable>StdInterpBasis</executable>
<parameters>-s quadrilateral -b GLL_Lagrange GLL_Lagrange -o 6 6 -p 7 7</parameters>
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-12">0</value>
</metric>
<metric type="Linf" id="2">
<value tolerance="1e-12">0</value>
</metric>
</metrics>
</test>
<?xml version="1.0" encoding="utf-8"?>
<test>
<description>StdInterpBasis Tet Modified basis P=4 Q=5</description>
<executable>StdInterpBasis</executable>
<parameters> -s tetrahedron -b Modified_A Modified_B Modified_C -o 4 4 4 -p 5 5 5</parameters>
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-12">1.89694e-15</value>
</metric>
<metric type="Linf" id="2">
<value tolerance="1e-12">1.77636e-14</value>
</metric>
</metrics>
</test>
<?xml version="1.0" encoding="utf-8"?>
<test>
<description>StdInterpBasis Triangle Modified basis P=7 Q=8</description>
<executable>StdInterpBasis</executable>
<parameters>-s triangle -b Modified_A Modified_B -o 7 7 -p 8 8</parameters>
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-12">0</value>
</metric>
<metric type="Linf" id="2">
<value tolerance="1e-12">0</value>
</metric>
</metrics>
</test>
<?xml version="1.0" encoding="utf-8"?>
<test>
<description>StdInterp Hex Modified basis P=6 Q=7</description>
<executable>StdInterp</executable>
<parameters>-s hexahedron -b Modified_A Modified_A Modified_A -o 6 6 6 -p 7 7 7</parameters>
<metrics>
<metric type="Linf" id="1">
<value tolerance="1e-12">0</value>
</metric>
<metric type="L2" id="2">
<value tolerance="1e-12">0</value>
</metric>
</metrics>
</test>
<?xml version="1.0" encoding="utf-8"?>
<test>
<description>StdInterp Prism Modified basis P=6 Q=7</description>
<executable>StdInterp</executable>
<parameters>-s prism -b Modified_A Modified_A Modified_B -o 6 6 6 -p 7 7 7</parameters>
<metrics>
<metric type="Linf" id="1">
<value tolerance="1e-12">0</value>
</metric>
<metric type="L2" id="2">
<value tolerance="1e-12">0</value>
</metric>
</metrics>
</test>
<?xml version="1.0" encoding="utf-8"?>
<test>
<description>StdInterpBasis Prism Orthonormal basis P=6 Q=7</description>
<executable>StdInterp</executable>
<parameters>-s prism -b Ortho_A Ortho_A Ortho_B -o 6 6 6 -p 7 7 7</parameters>
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-12">1.89694e-15</value>
</metric>
<metric type="Linf" id="2">
<value tolerance="1e-12">1.77636e-14</value>
</metric>
</metrics>
</test>
<?xml version="1.0" encoding="utf-8"?>
<test>
<description>StdInterp Pyramid Modified basis P=6 Q=7</description>
<executable>StdInterp</executable>
<parameters>-s pyramid -b Modified_A Modified_A ModifiedPyr_C -o 6 6 6 -p 7 7 7</parameters>
<metrics>
<metric type="Linf" id="1">
<value tolerance="1e-12">0</value>
</metric>
<metric type="L2" id="2">
<value tolerance="1e-12">0</value>
</metric>
</metrics>
</test>
<?xml version="1.0" encoding="utf-8"?>
<test>
<description>StdInterp Quadrilateral Lagrange basis P=6 Q=7</description>
<executable>StdInterp</executable>
<parameters>-s quadrilateral -b GLL_Lagrange GLL_Lagrange -o 6 6 -p 7 7</parameters>
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-12">0</value>
</metric>
<metric type="Linf" id="2">
<value tolerance="1e-12">0</value>
</metric>
</metrics>
</test>
<?xml version="1.0" encoding="utf-8"?>
<test>
<description>StdInterp Tetrahedron Modified basis P=6 Q=7</description>
<executable>StdInterp</executable>
<parameters>-s tetrahedron -b Modified_A Modified_B Modified_C -o 6 6 6 -p 7 7 7</parameters>
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-12">1.89776e-13</value>
</metric>
<metric type="Linf" id="2">
<value tolerance="2e-12">4.32543e-13</value>
</metric>
</metrics>
</test>
<?xml version="1.0" encoding="utf-8"?>
<test>
<description>StdInterp Triangle Modified basis P=6 Q=7</description>
<executable>StdInterp</executable>
<parameters>-s triangle -b Modified_A Modified_B -o 6 6 -p 7 7</parameters>
<metrics>
<metric type="Linf" id="1">
<value tolerance="1e-12">0</value>
</metric>
<metric type="L2" id="2">
<value tolerance="1e-12">0</value>
</metric>
</metrics>
</test>
......@@ -190,7 +190,7 @@ namespace Nektar
const BasisKey &rhs) const;
protected:
int m_nummodes; ///< Expansion order.
unsigned int m_nummodes; ///< Expansion order.
BasisType m_basistype; ///< Expansion type.
PointsKey m_pointsKey; ///< Points specification.
......@@ -277,6 +277,11 @@ namespace Nektar
m_points->GetZW(z,w);
}
inline const Array<OneD, const NekDouble>& GetBaryWeights() const
{
return m_points->GetBaryWeights();
}
inline const std::shared_ptr<NekMatrix<NekDouble> > & GetD(
Direction dir = xDir) const
{
......
......@@ -261,6 +261,7 @@ namespace Nektar
{
CalculatePoints();
CalculateWeights();
CalculateBaryWeights();
CalculateDerivMatrix();
}
......@@ -301,6 +302,11 @@ namespace Nektar
w = m_weights;
}
inline const Array<OneD, const NekDouble>& GetBaryWeights() const
{
return m_bcweights;
}
inline void GetPoints(Array<OneD, const DataType> &x) const
{
x = m_points[0];
......@@ -377,9 +383,16 @@ namespace Nektar
}
protected:
/// Points type for this points distributions.
PointsKey m_pointsKey;
/// Storage for the point locations, allowing for up to a 3D points
/// storage.
Array<OneD, DataType> m_points[3];
/// Quadrature weights for the weights.
Array<OneD, DataType> m_weights;
/// Barycentric weights.
Array<OneD, DataType> m_bcweights;
/// Derivative matrices.
MatrixSharedPtrType m_derivmatrix[3];
NekManager<PointsKey, NekMatrix<DataType>, PointsKey::opLess> m_InterpManager;
NekManager<PointsKey, NekMatrix<DataType>, PointsKey::opLess> m_GalerkinProjectionManager;
......@@ -400,6 +413,40 @@ namespace Nektar
m_weights = Array<OneD, DataType>(GetTotNumPoints());
}
/**
* @brief This function calculates the barycentric weights used for
* enhanced interpolation speed.
*
* For the points distribution \f$ z_i \f$ with \f% 1\leq z_i \leq N
* \f$, the barycentric weights are computed as:
*
* \f[
* b_i=\prod_{\substack{1\leq j\leq N\\ i\neq j}} \frac{1}{z_i-z_j}
* \f]
*/
virtual void CalculateBaryWeights()
{
const unsigned int totNumPoints = m_pointsKey.GetNumPoints();
m_bcweights = Array<OneD, DataType>(totNumPoints, 1.0);
Array<OneD, DataType> z = m_points[0];
for (unsigned int i = 0; i < totNumPoints; ++i)
{
for (unsigned int j = 0; j < totNumPoints; ++j)
{
if (i == j)
{
continue;
}
m_bcweights[i] *= (z[i] - z[j]);
}
m_bcweights[i] = 1.0 / m_bcweights[i];
}
}
virtual void CalculateDerivMatrix()
{
int totNumPoints = GetTotNumPoints();
......
This diff is collapsed.
......@@ -123,7 +123,6 @@ namespace Nektar
}
val = Vmath::Vamax(ntot, wsp, 1);
return val;
}
......@@ -143,7 +142,7 @@ namespace Nektar
Vmath::Vsub(ntot, sol, 1, phys, 1, wsp, 1);
Vmath::Vmul(ntot, wsp, 1, wsp, 1, wsp, 1);
}
val = v_Integral(wsp);
// if val too small, sqrt returns nan.
......@@ -1417,6 +1416,12 @@ namespace Nektar
return 0;
}
NekDouble StdExpansion::v_PhysEvaluateBasis(const Array<OneD, const NekDouble>& coords, int mode)
{
boost::ignore_unused(coords, mode);
NEKERROR(ErrorUtil::efatal, "Method does not exist for this shape");
return 0;
}
void StdExpansion::v_FillMode(const int mode, Array<OneD, NekDouble> &outarray)
{
......
......@@ -1198,6 +1198,24 @@ namespace Nektar
return v_PhysEvaluate(I,physvals);
}
/**
* @brief This function evaluates the basis function mode @p mode at a
* point @p coords of the domain.
*
* This function uses barycentric interpolation with the tensor
* product separation of the basis function to improve performance.
*
* @param coord The coordinate inside the standard region.
* @param mode The mode number to be evaluated.
*
* @return The value of the basis function @p mode at @p coords.
*/
NekDouble PhysEvaluateBasis(
const Array<OneD, const NekDouble>& coords,
int mode)
{
return v_PhysEvaluateBasis(coords, mode);
}
/**
* \brief Convert local cartesian coordinate \a xi into local
......@@ -1559,8 +1577,86 @@ namespace Nektar
StdRegions::Orientation dir);
STD_REGIONS_EXPORT virtual NekDouble v_StdPhysEvaluate(
const Array<OneD, const NekDouble> &Lcoord,
const Array<OneD, const NekDouble> &physvals);
const Array<OneD, const NekDouble> &coord,
const Array<OneD, const NekDouble> &physvals);
/**
* @brief This function performs the barycentric interpolation of
* the polynomial stored in @p coord at a point @p physvals using
* barycentric interpolation weights in direction @tparam DIR.
*
* This method is intended to be used a helper function for
* StdExpansion::PhysEvaluate and its elemental instances, so that
* the calling method should provide @p coord for x, y and z
* sequentially and the appropriate @p physvals and @p weights for
* that particular direction.
*
* @param coord The coordinate of the single point.
* @param physvals The polynomial stored at each quadrature point.
* @tparam DIR The direction of evaluation.
*
* @return The value of @p physvals at @p coord in direction @p dir.
*/
template<int DIR>
inline NekDouble BaryEvaluate(
const NekDouble &coord,
const NekDouble *physvals)
{
NekDouble numer = 0.0, denom = 0.0;
ASSERTL2(DIR < m_base.size(),
"Direction should be less than shape dimension.");
const Array<OneD, const NekDouble> &z = m_base[DIR]->GetZ();
const Array<OneD, const NekDouble> &bw =
m_base[DIR]->GetBaryWeights();
const auto nquad = z.size();
for (int i = 0; i < nquad; ++i)
{
NekDouble xdiff = z[i] - coord;
NekDouble pval = physvals[i];
/*
* (in this specific case) you actually
* want to do the comparison exactly
* (believe it or not!) See chapter 7 of
* the paper here:
*https://people.maths.ox.ac.uk/trefethen/barycentric.pdf
*/
if (xdiff == 0.0)
{
return pval;
}
NekDouble tmp = bw[i] / xdiff;
numer += tmp * pval;
denom += tmp;
}