Commit e1d92bb0 authored by Dave Moxey's avatar Dave Moxey
Browse files

Working prism and interpolation routines, added demo to test routines

parent 1e607194
......@@ -11,6 +11,8 @@ SET(PartitionAnalyseSources
SET(FoundationSources
FoundationDemo.cpp)
SET(NodalDemoSources NodalDemo.cpp)
SET(NodalTriFeketeSources
NodalTriFeketeDemo.cpp)
......@@ -41,6 +43,9 @@ TARGET_LINK_LIBRARIES(PartitionAnalyse LibUtilities)
ADD_NEKTAR_EXECUTABLE(FoundationDemo demos FoundationSources )
TARGET_LINK_LIBRARIES(FoundationDemo LibUtilities)
ADD_NEKTAR_EXECUTABLE(NodalDemo demos NodalDemoSources)
TARGET_LINK_LIBRARIES(NodalDemo LibUtilities)
ADD_NEKTAR_EXECUTABLE(NodalTriFeketeDemo demos NodalTriFeketeSources )
TARGET_LINK_LIBRARIES(NodalTriFeketeDemo LibUtilities)
......
///////////////////////////////////////////////////////////////////////////////
//
// File: NodalDemo.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: Demo for testing functionality of NodalUtil classes
//
///////////////////////////////////////////////////////////////////////////////
#include <iostream>
#include <string>
#include <boost/program_options.hpp>
#include <boost/algorithm/string.hpp>
#include <LibUtilities/Foundations/ManagerAccess.h>
#include <LibUtilities/Foundations/NodalUtil.h>
#include <LibUtilities/Foundations/Foundations.hpp>
#include <LibUtilities/Foundations/Points.h>
#include <LibUtilities/BasicUtils/ShapeType.hpp>
using namespace std;
using namespace Nektar;
using namespace Nektar::LibUtilities;
namespace po = boost::program_options;
int main(int argc, char *argv[])
{
po::options_description desc("Available options");
desc.add_options()
("help,h", "Produce this help message.")
("integral,i", "Evaluate the integral of a known function and "
"return error.")
("order,o", po::value<int>(),
"Order to evaluate at.")
("type,t", po::value<int>(),
"Points type of the nodal element.")
("deriv,d", "Evaluate the derivative of a known function and "
"return error.")
("interp,p", po::value<string>(),
"Interpolate function at a known point of the form"
"x,y,z and return error.");
po::variables_map vm;
try
{
po::store(po::command_line_parser(argc, argv).
options(desc).run(), vm);
po::notify(vm);
}
catch (const std::exception& e)
{
cerr << e.what() << endl;
cerr << desc;
return 1;
}
if (!vm.count("order"))
{
cerr << "Must supply an element order." << endl;
return 1;
}
int order = vm["order"].as<int>();
map<PointsType, ShapeType> nodalTypes;
map<PointsType, ShapeType>::iterator nIt;
nodalTypes[eNodalTriElec] = eTriangle;
nodalTypes[eNodalTriFekete] = eTriangle;
nodalTypes[eNodalTriEvenlySpaced] = eTriangle;
nodalTypes[eNodalTetEvenlySpaced] = eTetrahedron;
nodalTypes[eNodalTetElec] = eTetrahedron;
nodalTypes[eNodalPrismEvenlySpaced] = ePrism;
if (!vm.count("type"))
{
cerr << "Must supply a points type. Valid points types are:" << endl;
for (nIt = nodalTypes.begin(); nIt != nodalTypes.end(); ++nIt)
{
cerr << " - " << (int)nIt->first << " ("
<< kPointsTypeStr[nIt->first] << ")" << endl;
}
return 1;
}
PointsType pointsType = (PointsType)vm["type"].as<int>();
if (nodalTypes.find(pointsType) == nodalTypes.end())
{
cerr << "Invalid points type. Valid points types are:" << endl;
for (nIt = nodalTypes.begin(); nIt != nodalTypes.end(); ++nIt)
{
cerr << " - " << (int)nIt->first << " ("
<< kPointsTypeStr[nIt->first] << ")" << endl;
}
return 1;
}
// Generate nodal points.
PointsKey pointsKey(order + 1, pointsType);
PointsSharedPtr points = PointsManager()[pointsKey];
Array<OneD, NekDouble> r, s, t;
points->GetPoints(r, s, t);
// Generate nodal utility.
ShapeType shape = nodalTypes[pointsType];
NodalUtil *util = NULL;
if (shape == eTriangle)
{
util = new NodalUtilTriangle(order, r, s);
}
else if (shape == eTetrahedron)
{
util = new NodalUtilTetrahedron(order, r, s, t);
}
else if (shape == ePrism)
{
util = new NodalUtilPrism(order, r, s, t);
}
ASSERTL1(util, "Unknown shape type!");
const int nPoints = r.num_elements();
const int dim = shape == eTriangle ? 2 : 3;
if (vm.count("integral"))
{
NekVector<NekDouble> weights = util->GetWeights();
NekDouble integral = 0.0;
for (int i = 0; i < nPoints; ++i)
{
NekDouble integrand = dim == 2 ?
exp(r[i]) * sin(s[i]) : exp(r[i] + s[i] + t[i]);
integral += weights[i] * integrand;
}
NekDouble exact = 0.0;
switch(shape)
{
case eTriangle:
exact = -0.5 * (sin(1) + cos(1) + M_E * M_E * (sin(1) - cos(1)))
/ M_E;
break;
case eTetrahedron:
exact = 1.0 / M_E - 1.0 / M_E / M_E / M_E;
break;
case ePrism:
exact = M_E - 1.0 / M_E / M_E / M_E;
break;
}
cout << "L inf error : " << fabs(exact - integral) << endl;
}
else if (vm.count("deriv"))
{
Array<OneD, NekVector<NekDouble> > exact(dim);
NekVector<NekDouble> input(nPoints), output(nPoints);
for (int i = 0; i < dim; ++i)
{
exact[i] = NekVector<NekDouble>(nPoints);
}
if (dim == 2)
{
// Exact solution: u(x,y) = sin(x) * cos(y)
for (int i = 0; i < nPoints; ++i)
{
input[i] = sin(r[i]) * cos(s[i]);
exact[0][i] = cos(r[i]) * cos(s[i]);
exact[1][i] = -sin(r[i]) * sin(s[i]);
}
}
else
{
// Exact solution: u(x,y) = sin(x) * cos(y) * sin(z)
for (int i = 0; i < nPoints; ++i)
{
input[i] = sin(r[i]) * cos(s[i]) * sin(t[i]);
exact[0][i] = cos(r[i]) * cos(s[i]) * sin(t[i]);
exact[1][i] = -sin(r[i]) * sin(s[i]) * sin(t[i]);
exact[2][i] = sin(r[i]) * cos(s[i]) * cos(t[i]);
}
}
char *vars[3] = { "u", "v", "w" };
for (int i = 0; i < dim; ++i)
{
boost::shared_ptr<NekMatrix<NekDouble> > deriv =
util->GetDerivMatrix(i);
output = *deriv * input;
NekVector<NekDouble> tmp = output - exact[i];
cout << "L 2 error (variable " << vars[i] << ") : " << tmp.L2Norm()
<< endl;
}
}
else if (vm.count("interp"))
{
string pointStr = vm["interp"].as<string>();
vector<string> point;
boost::split(point, pointStr, boost::is_any_of(","));
if (point.size() != dim)
{
cerr << "Point " << pointStr << " does not have correct dimensions"
<< endl;
return 1;
}
Array<OneD, Array<OneD, NekDouble> > tmp(dim);
for (int i = 0; i < dim; ++i)
{
tmp[i] = Array<OneD, NekDouble>(1);
try
{
tmp[i][0] = boost::lexical_cast<NekDouble>(point[i]);
}
catch(boost::bad_lexical_cast &)
{
cerr << "Could not convert " << point[i] << " to a coordinate"
<< endl;
return 1;
}
}
NekVector<NekDouble> input(nPoints);
for (int i = 0; i < nPoints; ++i)
{
input[i] = dim == 2 ? exp(r[i]) * exp(s[i]) :
exp(r[i]) * exp(s[i]) * exp(t[i]);
}
boost::shared_ptr<NekMatrix<NekDouble> > interp =
util->GetInterpolationMatrix(tmp);
NekVector<NekDouble> output = *interp * input;
NekDouble exact = dim == 2 ? exp(tmp[0][0]) * exp(tmp[1][0]) :
exp(tmp[0][0]) * exp(tmp[1][0]) * exp(tmp[2][0]);
cout << "L inf error : " << fabs(exact - output[0]) << endl;
}
else
{
cerr << "You must specify one of --integral, --deriv or --interp"
<< endl;
cerr << desc;
return 1;
}
return 0;
}
......@@ -50,6 +50,17 @@ namespace Nektar
namespace LibUtilities
{
/**
* @brief Obtain the integration weights for the given nodal distribution.
*
* This routine constructs the integration weights for the given nodal
* distribution inside NodalUtil::m_xi. To do this we solve the linear system
* \f$ \mathbf{V}^\top \mathbf{w} = \mathbf{g} \f$ where \f$ \mathbf{g}_i =
* \int_\Omega \psi_i(x)\, dx \f$, and we use the fact that under the definition
* of the orthogonal basis for each element type, \f$ \mathbf{g}_i = 0 \f$ for
* \f$ i > 0 \f$. We use NodalUtil::v_ModeZeroIntegral to return the analytic
* value of \f$ \mathbf{g}_0 \f$.
*/
NekVector<NekDouble> NodalUtil::GetWeights()
{
// Get number of modes in orthogonal basis
......@@ -60,7 +71,7 @@ NekVector<NekDouble> NodalUtil::GetWeights()
if (numModes == m_numPoints)
{
NekVector<NekDouble> g(m_numPoints, 0.0);
g(0) = 2.0*sqrt(2.0);
g(0) = v_ModeZeroIntegral();
SharedMatrix V = GetVandermonde();
......@@ -76,6 +87,14 @@ NekVector<NekDouble> NodalUtil::GetWeights()
}
}
/**
* @brief Return the Vandermonde matrix for the nodal distribution.
*
* This routine constructs and returns the Vandermonde matrix \f$\mathbf{V}\f$,
* with each entry as \f$\mathbf{V}_{ij} = (\psi_i(\xi_j))\f$ where \f$ \psi_i
* \f$ is the orthogonal basis obtained through the abstract function
* NodalUtil::v_OrthoBasis.
*/
SharedMatrix NodalUtil::GetVandermonde()
{
int rows = m_numPoints, cols = v_NumModes();
......@@ -93,7 +112,19 @@ SharedMatrix NodalUtil::GetVandermonde()
return matV;
}
/**
* @brief Return the Vandermonde matrix of the derivative of the basis functions
* for the nodal distribution.
*
* This routine constructs and returns the Vandermonde matrix for the derivative
* of the basis functions \f$\mathbf{V}_d\f$ for coordinate directions \f$ 0
* \leq d \leq 2 \f$, with each entry as \f$\mathbf{V}_{ij} = (\partial_d
* \psi_i(\xi_j))\f$ where \f$ \partial_d\psi_i \f$ is the derivative of the
* orthogonal basis obtained through the abstract function
* NodalUtil::v_OrthoBasisDeriv.
*
* @param dir Direction of derivative in the standard element.
*/
SharedMatrix NodalUtil::GetVandermondeForDeriv(int dir)
{
int rows = m_numPoints, cols = v_NumModes();
......@@ -112,6 +143,17 @@ SharedMatrix NodalUtil::GetVandermondeForDeriv(int dir)
return matV;
}
/**
* @brief Return the derivative matrix for the nodal distribution.
*
* This routine constructs and returns the derivative matrices
* \f$\mathbf{D}_d\f$ for coordinate directions \f$ 0 \leq d \leq 2 \f$, which
* can be used to evaluate the derivative of a nodal expansion at the points
* defined by NodalUtil::m_xi. These are calculated as \f$ \mathbf{D}_d =
* \mathbf{V}_d \mathbf{V}^{-1} \f$, where \f$ \mathbf{V}_d \f$ is the
* derivative Vandermonde matrix and \f$ \mathbf{V} \f$ is the Vandermonde
* matrix.
*/
SharedMatrix NodalUtil::GetDerivMatrix(int dir)
{
SharedMatrix V = GetVandermonde();
......@@ -126,15 +168,64 @@ SharedMatrix NodalUtil::GetDerivMatrix(int dir)
return D;
}
NodalUtilTriangle::NodalUtilTriangle(int degree,
Array<OneD, NekDouble> &r,
Array<OneD, NekDouble> &s)
/**
* @brief Construct the interpolation matrix used to evaluate the basis at the
* points @p xi inside the element.
*
* This routine returns a matrix \f$ \mathbf{I} \f$ that can be used to evaluate
* the nodal basis at the points defined by the parameter @p xi, In particular
* if the array \f$ \mathbf{u} \f$ with components \f$ \mathbf{u}_i =
* u^\delta(\xi_i) \f$ represents the function
*
* @param xi An array of first size number of spatial dimensions \f$ d \f$ and
* secondary size the number of points to
*/
SharedMatrix NodalUtil::GetInterpolationMatrix(
Array<OneD, Array<OneD, NekDouble> > &xi)
{
boost::shared_ptr<NodalUtil> subUtil = v_CreateUtil(xi);
SharedMatrix matS = GetVandermonde();
SharedMatrix matT = subUtil->GetVandermonde();
SharedMatrix D = MemoryManager<NekMatrix<NekDouble> >::AllocateSharedPtr(
matT->GetRows(), matS->GetColumns(), 0.0);
matS->Invert();
*D = (*matT) * (*matS);
return D;
}
/**
* @brief Construct the nodal utility class for a triangle.
*
* The constructor of this class sets up two important member variables:
*
* - NodalUtilTriangle::m_eta is used to construct the collapsed coordinate
* locations of the nodal points \f$ (\eta_1, \eta_2) \f$ inside the square
* \f$[-1,1]\f$
* - NodalUtilTriangle::m_ordering constructs a mapping from the index set \f$ I
* = \{ (i,j)\ |\ 0\leq i,j \leq P, i+j \leq P \}\f$ to an ordering \f$ 0 \leq
* m(ij) \leq (P+1)(P+2)/2 \f$ that defines the monomials \f$ \xi_1^i \xi_2^j
* \f$ that span the triangular space. This is then used to calculate which
* \f$ (i,j) \f$ pair corresponds to a column of the Vandermonde matrix when
* calculating the orthogonal polynomials.
*
* @param degree Polynomial order of this nodal triangle.
* @param r \f$ r \f$-coordinates of nodal points in the standard element.
* @param s \f$ s \f$-coordinates of nodal points in the standard element.
*/
NodalUtilTriangle::NodalUtilTriangle(int degree,
Array<OneD, NekDouble> r,
Array<OneD, NekDouble> s)
: NodalUtil(degree, 2), m_eta(2)
{
// Set up parent variables.
m_numPoints = r.num_elements();
m_xi[0] = r;
m_xi[1] = s;
// Construct a mapping (i,j) -> m from the triangular tensor product space
// (i,j) to a single ordering m.
for (int i = 0; i <= m_degree; ++i)
{
for (int j = 0; j <= m_degree - i; ++j)
......@@ -162,6 +253,20 @@ NodalUtilTriangle::NodalUtilTriangle(int degree,
}
}
/**
* @brief Return the value of the modal functions for the triangular element at
* the nodal points #m_xi for a given mode.
*
* In a triangle, we use the orthogonal basis
*
* \f[
* \psi_{m(ij)} = \sqrt{2} P^{(0,0)}_i(\xi_1) P_j^{(2i+1,0)}(\xi_2) (1-\xi_2)^i
* \f]
*
* where \f$ m(ij) \f$ is the mapping defined in NodalUtilTriangle::m_ordering.
*
* @param mode The mode of the orthogonal basis to evaluate.
*/
NekVector<NekDouble> NodalUtilTriangle::v_OrthoBasis(const int mode)
{
std::vector<NekDouble> jacobi_i(m_numPoints), jacobi_j(m_numPoints);
......@@ -186,6 +291,21 @@ NekVector<NekDouble> NodalUtilTriangle::v_OrthoBasis(const int mode)
return ret;
}
/**
* @brief Return the value of the derivative of the modal functions for the
* triangular element at the nodal points #m_xi for a given mode.
*
* In a triangle, we use the orthogonal basis
*
* \f[
* \psi_{m(ij)} = \sqrt{2} P^{(0,0)}_i(\xi_1) P_j^{(2i+1,0)}(\xi_2) (1-\xi_2)^i
* \f]
*
* where \f$ m(ij) \f$ is the mapping defined in
* NodalUtilTriangle::m_ordering. The derivative is then evaluated
*
* @param mode The mode of the orthogonal basis to evaluate.
*/
NekVector<NekDouble> NodalUtilTriangle::v_OrthoBasisDeriv(
const int dir, const int mode)
{
......@@ -246,10 +366,29 @@ NekVector<NekDouble> NodalUtilTriangle::v_OrthoBasisDeriv(
return ret;
}
NodalUtilTetrahedron::NodalUtilTetrahedron(int degree,
Array<OneD, NekDouble> &r,
Array<OneD, NekDouble> &s,
Array<OneD, NekDouble> &t)
/**
* @brief Construct the nodal utility class for a tetrahedron.
*
* The constructor of this class sets up two important member variables:
*
* - NodalUtilTriangle::m_eta is used to construct the collapsed coordinate
* locations of the nodal points \f$ (\eta_1, \eta_2) \f$ inside the square
* \f$[-1,1]\f$
* - NodalUtilTriangle::m_ordering constructs a mapping from the index set \f$ I
* = \{ (i,j)\ |\ 0\leq i,j \leq P, i+j \leq P \}\f$ to an ordering \f$ 0 \leq
* m(ij) \leq (P+1)(P+2)/2 \f$ that defines the monomials \f$ \xi_1^i \xi_2^j
* \f$ that span the triangular space. This is then used to calculate which
* \f$ (i,j) \f$ pair corresponds to a column of the Vandermonde matrix when
* calculating the orthogonal polynomials.
*
* @param degree Polynomial order of this nodal triangle.
* @param r \f$ r \f$-coordinates of nodal points in the standard element.
* @param s \f$ s \f$-coordinates of nodal points in the standard element.
*/
NodalUtilTetrahedron::NodalUtilTetrahedron(int degree,
Array<OneD, NekDouble> r,
Array<OneD, NekDouble> s,
Array<OneD, NekDouble> t)
: NodalUtil(degree, 3), m_eta(3)
{
m_numPoints = r.num_elements();
......@@ -426,6 +565,152 @@ NekVector<NekDouble> NodalUtilTetrahedron::v_OrthoBasisDeriv(
return ret;
}
NodalUtilPrism::NodalUtilPrism(int degree,
Array<OneD, NekDouble> r,
Array<OneD, NekDouble> s,
Array<OneD, NekDouble> t)
: NodalUtil(degree, 3), m_eta(3)
{
m_numPoints = r.num_elements();
m_xi[0] = r;
m_xi[1] = s;
m_xi[2] = t;
for (int i = 0; i <= m_degree; ++i)
{
for (int j = 0; j <= m_degree; ++j)
{
for (int k = 0; k <= m_degree - i; ++k)
{
m_ordering.push_back(Mode(i, j, k));
}
}
}
// Calculate collapsed coordinates from r/s values
m_eta[0] = Array<OneD, NekDouble>(m_numPoints);
m_eta[1] = Array<OneD, NekDouble>(m_numPoints);
m_eta[2] = Array<OneD, NekDouble>(m_numPoints);
for (int i = 0; i < m_numPoints; ++i)
{
if (fabs(m_xi[2][i] - 1.0) < NekConstants::kNekZeroTol)
{
// Very top point of the prism
m_eta[0][i] = -1.0;
m_eta[1][i] = m_xi[1][i];
m_eta[2][i] = 1.0;
}
else
{
// Third basis function collapsed to "pr" direction instead of "qr"
// direction
m_eta[0][i] = 2.0*(1.0 + m_xi[0][i])/(1.0 - m_xi[2][i]) - 1.0;
m_eta[1][i] = m_xi[1][i];
m_eta[2][i] = m_xi[2][i];
}
}
}
NekVector<NekDouble> NodalUtilPrism::v_OrthoBasis(const int mode)
{
std::vector<NekDouble> jacA(m_numPoints), jacB(m_numPoints);
std::vector<NekDouble> jacC(m_numPoints);
Mode modes = m_ordering[mode];
const int I = modes.get<0>(), J = modes.get<1>(), K = modes.get<2>();
// Calculate Jacobi polynomials
Polylib::jacobfd(
m_numPoints, &m_eta[0][0], &jacA[0], NULL, I, 0.0, 0.0);
Polylib::jacobfd(
m_numPoints, &m_eta[1][0], &jacB[0], NULL, J, 0.0, 0.0);
Polylib::jacobfd(
m_numPoints, &m_eta[2][0], &jacC[0], NULL, K, 2.0 * I + 1.0, 0.0);
NekVector<NekDouble> ret(m_numPoints);
NekDouble sqrt2 = sqrt(2.0);
for (int i = 0; i < m_numPoints; ++i)
{
ret[i] = sqrt2 * jacA[i] * jacB[i] * jacC[i] *