Commit 8954c9aa authored by Dave Moxey's avatar Dave Moxey
Browse files

Finish documentation, add tests

parent c8bb1e1a
......@@ -34,3 +34,12 @@ TARGET_LINK_LIBRARIES(NodalDemo LibUtilities)
ADD_NEKTAR_EXECUTABLE(TimeIntegrationDemo demos TimeIntegrationDemoSources)
TARGET_LINK_LIBRARIES(TimeIntegrationDemo LibUtilities)
ADD_NEKTAR_TEST(NodalDemo_Tri_Deriv_P8)
ADD_NEKTAR_TEST(NodalDemo_Tri_Integral_P6)
ADD_NEKTAR_TEST(NodalDemo_Tri_Interp_P7)
ADD_NEKTAR_TEST(NodalDemo_Prism_Deriv_P8)
ADD_NEKTAR_TEST(NodalDemo_Prism_Integral_P6)
ADD_NEKTAR_TEST(NodalDemo_Prism_Interp_P7)
ADD_NEKTAR_TEST(NodalDemo_Tet_Deriv_P8)
ADD_NEKTAR_TEST(NodalDemo_Tet_Integral_P6)
ADD_NEKTAR_TEST(NodalDemo_Tet_Interp_P7)
......@@ -181,6 +181,9 @@ int main(int argc, char *argv[])
case ePrism:
exact = M_E - 1.0 / M_E / M_E / M_E;
break;
default:
exact = 0.0;
break;
}
cout << "L inf error : " << fabs(exact - integral) << endl;
......
<?xml version="1.0" encoding="utf-8" ?>
<test>
<description>Nodal prism integration, evenly spaced points, P = 8</description>
<executable>NodalDemo</executable>
<parameters>--order 8 --type 27 --deriv</parameters>
<metrics>
<metric type="L2" id="1">
<value variable="u" tolerance="1e-12">0.000197703</value>
<value variable="v" tolerance="1e-12">5.30087e-07</value>
<value variable="w" tolerance="1e-12">0.000197703</value>
</metric>
</metrics>
</test>
<?xml version="1.0" encoding="utf-8" ?>
<test>
<description>Nodal prism integration, evenly spaced points, P = 6</description>
<executable>NodalDemo</executable>
<parameters>--order 6 --type 27 --integral</parameters>
<metrics>
<metric type="Linf" id="1">
<value tolerance="1e-12">3.00544e-06</value>
</metric>
</metrics>
</test>
<?xml version="1.0" encoding="utf-8" ?>
<test>
<description>Nodal prism interpolation, evenly spaced points, P = 7</description>
<executable>NodalDemo</executable>
<parameters>--order 7 --type 27 --interp -0.3,-0.636,-0.9</parameters>
<metrics>
<metric type="Linf" id="1">
<value tolerance="1e-12">3.18397e-08</value>
</metric>
</metrics>
</test>
<?xml version="1.0" encoding="utf-8" ?>
<test>
<description>Nodal tetrahedron integration, electrostatic points, P = 8</description>
<executable>NodalDemo</executable>
<parameters>--order 8 --type 26 --deriv</parameters>
<metrics>
<metric type="L2" id="1">
<value variable="u" tolerance="1e-12">0.000519629</value>
<value variable="v" tolerance="1e-12">0.000577377</value>
<value variable="w" tolerance="1e-12">0.000519629</value>
</metric>
</metrics>
</test>
<?xml version="1.0" encoding="utf-8" ?>
<test>
<description>Nodal tetrahedron integration, electrostatic points, P = 6</description>
<executable>NodalDemo</executable>
<parameters>--order 6 --type 26 --integral</parameters>
<metrics>
<metric type="Linf" id="1">
<value tolerance="1e-12">2.69479e-07</value>
</metric>
</metrics>
</test>
<?xml version="1.0" encoding="utf-8" ?>
<test>
<description>Nodal tetrahedron interpolation, evenly spaced points, P = 7</description>
<executable>NodalDemo</executable>
<parameters>--order 7 --type 25 --interp -0.3,-0.636,-0.9</parameters>
<metrics>
<metric type="Linf" id="1">
<value tolerance="1e-12">1.66676e-09</value>
</metric>
</metrics>
</test>
<?xml version="1.0" encoding="utf-8" ?>
<test>
<description>Nodal triangle integration, electrostatic points, P = 8</description>
<executable>NodalDemo</executable>
<parameters>--order 8 --type 22 --deriv</parameters>
<metrics>
<metric type="L2" id="1">
<value variable="u" tolerance="1e-12">0.000192327</value>
<value variable="v" tolerance="1e-12">0.000192208</value>
</metric>
</metrics>
</test>
<?xml version="1.0" encoding="utf-8" ?>
<test>
<description>Nodal triangle integration, Fekete points, P = 6</description>
<executable>NodalDemo</executable>
<parameters>--order 6 --type 23 --integral</parameters>
<metrics>
<metric type="Linf" id="1">
<value tolerance="1e-12">7.28803e-07</value>
</metric>
</metrics>
</test>
<?xml version="1.0" encoding="utf-8" ?>
<test>
<description>Nodal triangle interpolation, evenly spaced points, P = 7</description>
<executable>NodalDemo</executable>
<parameters>--order 7 --type 24 --interp -0.3,-0.636</parameters>
<metrics>
<metric type="Linf" id="1">
<value tolerance="1e-12">1.38282e-08</value>
</metric>
</metrics>
</test>
......@@ -172,13 +172,19 @@ SharedMatrix NodalUtil::GetDerivMatrix(int dir)
* @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
* This routine returns a matrix \f$ \mathbf{I}(\mathbf{a}) \f$ that can be used
* to evaluate the nodal basis at the points defined by the parameter @p xi,
* which is denoted by \f$ \mathbf{a} = (a_1, \dots, a_N) \f$ and \f$ N \f$ is
* the number of points in @p xi.
*
* In particular, if the array \f$ \mathbf{u} \f$ with components \f$
* \mathbf{u}_i = u^\delta(\xi_i) \f$ represents the polynomial approximation of
* a function \f$ u \f$ evaluated at the nodal points NodalUtil::m_xi, then the
* evaluation of \f$ u^\delta \f$ evaluated at the input points \f$ \mathbf{a}
* \f$ is given by \f$ \mathbf{I}(\mathbf{a})\mathbf{u} \f$.
*
* @param xi An array of first size number of spatial dimensions \f$ d \f$ and
* secondary size the number of points to
* secondary size the number of points to interpolate.
*/
SharedMatrix NodalUtil::GetInterpolationMatrix(
Array<OneD, Array<OneD, NekDouble> > &xi)
......@@ -198,21 +204,24 @@ SharedMatrix NodalUtil::GetInterpolationMatrix(
/**
* @brief Construct the nodal utility class for a triangle.
*
* The constructor of this class sets up two important member variables:
* The constructor of this class sets up two member variables used in the
* evaluation of the orthogonal basis:
*
* - 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$
* \f$[-1,1]^2\f$ on which the orthogonal basis functions are defined.
* - 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
* \f$ (i,j) \f$ pair corresponding 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.
* @param r \f$ \xi_1 \f$-coordinates of nodal points in the standard
* element.
* @param s \f$ \xi_2 \f$-coordinates of nodal points in the standard
* element.
*/
NodalUtilTriangle::NodalUtilTriangle(int degree,
Array<OneD, NekDouble> r,
......@@ -263,7 +272,8 @@ NodalUtilTriangle::NodalUtilTriangle(int degree,
* \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.
* where \f$ m(ij) \f$ is the mapping defined in NodalUtilTriangle::m_ordering
* and \f$ J_n^{(\alpha,\beta)}(z) \f$ denotes the standard Jacobi polynomial.
*
* @param mode The mode of the orthogonal basis to evaluate.
*/
......@@ -295,14 +305,9 @@ NekVector<NekDouble> NodalUtilTriangle::v_OrthoBasis(const int mode)
* @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
* Note that this routine must use the chain rule combined with the collapsed
* coordinate derivatives as described in Sherwin & Karniadakis (2nd edition),
* pg 150.
*
* @param mode The mode of the orthogonal basis to evaluate.
*/
......@@ -333,6 +338,7 @@ NekVector<NekDouble> NodalUtilTriangle::v_OrthoBasisDeriv(
if (dir == 0)
{
// d/d(\xi_1) = 2/(1-\eta_2) d/d(\eta_1)
for (int i = 0; i < m_numPoints; ++i)
{
ret[i] = 2.0 * sqrt2 * jacobi_di[i] * jacobi_j[i];
......@@ -344,6 +350,7 @@ NekVector<NekDouble> NodalUtilTriangle::v_OrthoBasisDeriv(
}
else
{
// d/d(\xi_2) = 2(1+\eta_1)/(1-\eta_2) d/d(\eta_1) + d/d(eta_2)
for (int i = 0; i < m_numPoints; ++i)
{
ret[i] = (1 + m_eta[0][i]) * sqrt2 * jacobi_di[i] * jacobi_j[i];
......@@ -369,21 +376,27 @@ NekVector<NekDouble> NodalUtilTriangle::v_OrthoBasisDeriv(
/**
* @brief Construct the nodal utility class for a tetrahedron.
*
* The constructor of this class sets up two important member variables:
* The constructor of this class sets up two member variables used in the
* evaluation of the orthogonal basis:
*
* - 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.
* - NodalUtilTetrahedron::m_eta is used to construct the collapsed coordinate
* locations of the nodal points \f$ (\eta_1, \eta_2, \eta_3) \f$ inside the
* cube \f$[-1,1]^3\f$ on which the orthogonal basis functions are defined.
* - NodalUtilTetrahedron::m_ordering constructs a mapping from the index set
* \f$ I = \{ (i,j,k)\ |\ 0\leq i,j,k \leq P, i+j \leq P, i+j+k \leq P \}\f$
* to an ordering \f$ 0 \leq m(ijk) \leq (P+1)(P+2)(P+3)/6 \f$ that defines
* the monomials \f$ \xi_1^i \xi_2^j \xi_3^k \f$ that span the tetrahedral
* space. This is then used to calculate which \f$ (i,j,k) \f$ triple
* (represented as a boost tuple) corresponding 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.
* @param degree Polynomial order of this nodal tetrahedron
* @param r \f$ \xi_1 \f$-coordinates of nodal points in the standard
* element.
* @param s \f$ \xi_2 \f$-coordinates of nodal points in the standard
* element.
* @param t \f$ \xi_3 \f$-coordinates of nodal points in the standard
* element.
*/
NodalUtilTetrahedron::NodalUtilTetrahedron(int degree,
Array<OneD, NekDouble> r,
......@@ -444,6 +457,20 @@ NodalUtilTetrahedron::NodalUtilTetrahedron(int degree,
}
}
/**
* @brief Return the value of the modal functions for the tetrahedral element at
* the nodal points #m_xi for a given mode.
*
* In a tetrahedron, we use the orthogonal basis
*
* \f[ \psi_{m(ijk)} = \sqrt{8} P^{(0,0)}_i(\xi_1) P_j^{(2i+1,0)}(\xi_2)
* P_k^{(2i+2j+2,0)}(\xi_3) (1-\xi_2)^i (1-\xi_3)^{i+j} \f]
*
* where \f$ m(ijk) \f$ is the mapping defined in #m_ordering and \f$
* J_n^{(\alpha,\beta)}(z) \f$ denotes the standard Jacobi polynomial.
*
* @param mode The mode of the orthogonal basis to evaluate.
*/
NekVector<NekDouble> NodalUtilTetrahedron::v_OrthoBasis(const int mode)
{
std::vector<NekDouble> jacA(m_numPoints), jacB(m_numPoints);
......@@ -472,6 +499,16 @@ NekVector<NekDouble> NodalUtilTetrahedron::v_OrthoBasis(const int mode)
return ret;
}
/**
* @brief Return the value of the derivative of the modal functions for the
* tetrahedral element at the nodal points #m_xi for a given mode.
*
* Note that this routine must use the chain rule combined with the collapsed
* coordinate derivatives as described in Sherwin & Karniadakis (2nd edition),
* pg 152.
*
* @param mode The mode of the orthogonal basis to evaluate.
*/
NekVector<NekDouble> NodalUtilTetrahedron::v_OrthoBasisDeriv(
const int dir, const int mode)
{
......@@ -565,6 +602,31 @@ NekVector<NekDouble> NodalUtilTetrahedron::v_OrthoBasisDeriv(
return ret;
}
/**
* @brief Construct the nodal utility class for a prism.
*
* The constructor of this class sets up two member variables used in the
* evaluation of the orthogonal basis:
*
* - NodalUtilPrism::m_eta is used to construct the collapsed coordinate
* locations of the nodal points \f$ (\eta_1, \eta_2, \eta_3) \f$ inside the
* cube \f$[-1,1]^3\f$ on which the orthogonal basis functions are defined.
* - NodalUtilPrism::m_ordering constructs a mapping from the index set
* \f$ I = \{ (i,j,k)\ |\ 0\leq i,j,k \leq P, i+k \leq P \}\f$ to an ordering
* \f$ 0 \leq m(ijk) \leq (P+1)(P+1)(P+2)/2 \f$ that defines the monomials \f$
* \xi_1^i \xi_2^j \xi_3^k \f$ that span the prismatic space. This is then
* used to calculate which \f$ (i,j,k) \f$ triple (represented as a boost
* tuple) corresponding to a column of the Vandermonde matrix when calculating
* the orthogonal polynomials.
*
* @param degree Polynomial order of this nodal tetrahedron
* @param r \f$ \xi_1 \f$-coordinates of nodal points in the standard
* element.
* @param s \f$ \xi_2 \f$-coordinates of nodal points in the standard
* element.
* @param t \f$ \xi_3 \f$-coordinates of nodal points in the standard
* element.
*/
NodalUtilPrism::NodalUtilPrism(int degree,
Array<OneD, NekDouble> r,
Array<OneD, NekDouble> s,
......@@ -612,6 +674,20 @@ NodalUtilPrism::NodalUtilPrism(int degree,
}
}
/**
* @brief Return the value of the modal functions for the prismatic element at
* the nodal points #m_xi for a given mode.
*
* In a prism, we use the orthogonal basis
*
* \f[ \psi_{m(ijk)} = \sqrt{2} P^{(0,0)}_i(\xi_1) P_j^{(0,0)}(\xi_2)
* P_k^{(2i+1,0)}(\xi_3) (1-\xi_3)^i \f]
*
* where \f$ m(ijk) \f$ is the mapping defined in #m_ordering and \f$
* J_n^{(\alpha,\beta)}(z) \f$ denotes the standard Jacobi polynomial.
*
* @param mode The mode of the orthogonal basis to evaluate.
*/
NekVector<NekDouble> NodalUtilPrism::v_OrthoBasis(const int mode)
{
std::vector<NekDouble> jacA(m_numPoints), jacB(m_numPoints);
......@@ -640,6 +716,16 @@ NekVector<NekDouble> NodalUtilPrism::v_OrthoBasis(const int mode)
return ret;
}
/**
* @brief Return the value of the derivative of the modal functions for the
* prismatic element at the nodal points #m_xi for a given mode.
*
* Note that this routine must use the chain rule combined with the collapsed
* coordinate derivatives as described in Sherwin & Karniadakis (2nd edition),
* pg 152.
*
* @param mode The mode of the orthogonal basis to evaluate.
*/
NekVector<NekDouble> NodalUtilPrism::v_OrthoBasisDeriv(
const int dir, const int mode)
{
......
......@@ -64,18 +64,28 @@ typedef boost::shared_ptr<NekMatrix<NekDouble> > SharedMatrix;
* The NodalUtil class and its subclasses are designed to take care of some
* common issues that arise when considering triangles, tetrahedra and prismatic
* elements that are equipped with a nodal Lagrangian basis, defined using a set
* of nodal points that we store in the array NodalUtil::m_xi. Since one cannot
* write this basis analytically, we instead construct the Vandermonde matrix
* of nodal points \f$ \xi_i \f$ that we store in the array
* NodalUtil::m_xi. Since one cannot write this basis analytically, we instead
* construct the Vandermonde matrix
*
* \f[ \mathbf{V}_{ij} \f]
* \f[ \mathbf{V}_{ij} = \psi_j(\xi_i) \f]
*
* where \f$ \psi_j \f$ is a basis that spans the polynomial space of the
* element. The Vandermonde matrix can then be used to construct the integration
* weights, derivative and interpolation matrices. Although this can be any
* basis, such as the monomial basis \f$ x^i y^j z^k \f$, in practice this is
* numerically unstable at high polynomial orders. Elements are therefore
* expected to use the 'traditional' modal orthogonal basis. See Sherwin &
* Karniadakis or Hesthaven & Warburton for further details of this basis and
* the surrounding numerical issues.
*
* This class therefore contains the generic logic needed to construct various
* matrices, and subclasses override virtual functions that define the
* orthogonal basis and its derivatives for a particular element type.
*/
class NodalUtil
{
public:
NodalUtil(int degree, int dim) : m_dim(dim), m_degree(degree), m_xi(dim)
{
}
LIB_UTILITIES_EXPORT NekVector<NekDouble> GetWeights();
LIB_UTILITIES_EXPORT SharedMatrix GetVandermonde();
LIB_UTILITIES_EXPORT SharedMatrix GetVandermondeForDeriv(int dir);
......@@ -84,20 +94,78 @@ public:
Array<OneD, Array<OneD, NekDouble> > &xi);
protected:
/**
* @brief Set up the NodalUtil object.
*
* @param dim Dimension of the element.
* @param degree Polynomial degree of the element.
*/
NodalUtil(int degree, int dim) : m_dim(dim), m_degree(degree), m_xi(dim)
{
}
/// Dimension of the nodal element
int m_dim;
/// Degree of the nodal element
int m_degree;
/// Total number of nodal points
int m_numPoints;
/// Coordinates of the nodal points defining the basis
Array<OneD, Array<OneD, NekDouble> > m_xi;
/**
* @brief Return the values of the orthogonal basis at the nodal points for
* a given mode.
*
* @param dir Coordinate direction of derivative.
* @param mode Mode number, which is between 0 and NodalUtil::v_NumModes()
* - 1.
*/
virtual NekVector<NekDouble> v_OrthoBasis(const int mode) = 0;
/**
* @brief Return the values of the derivative of the orthogonal basis at the
* nodal points for a given mode.
*
* @param dir Coordinate direction of derivative.
* @param mode Mode number, which is between 0 and NodalUtil::v_NumModes()
* - 1.
*/
virtual NekVector<NekDouble> v_OrthoBasisDeriv(
const int dir, const int mode) = 0;
/**
* @brief Construct a NodalUtil object of the appropriate element type for a
* given set of points.
*
* This function is used inside NodalUtil::GetInterpolationMatrix so that
* the (potentially non-square) Vandermonde matrix can be constructed to
* create the interpolation matrix at an arbitrary set of points in the
* domain.
*/
virtual boost::shared_ptr<NodalUtil> v_CreateUtil(
Array<OneD, Array<OneD, NekDouble> > &xi) = 0;
/**
* @brief Return the value of the integral of the zero-th mode for this
* element.
*
* Note that for the orthogonal basis under consideration, all modes
* integrate to zero asides from the zero-th mode. This function is used in
* NodalUtil::GetWeights to determine integration weights.
*/
virtual NekDouble v_ModeZeroIntegral() = 0;
/**
* @brief Calculate the number of degrees of freedom for this element.
*/
virtual int v_NumModes() = 0;
};
/**
* @brief Specialisation of the NodalUtil class to support nodal triangular
* elements.
*/
class NodalUtilTriangle : public NodalUtil
{
public:
......@@ -106,7 +174,11 @@ public:
Array<OneD, NekDouble> s);
protected:
/// Mapping from the \f$ (i,j) \f$ indexing of the basis to a continuous
/// ordering.
std::vector<std::pair<int, int> > m_ordering;
/// Collapsed coordinates \f$ (\eta_1, \eta_2) \f$ of the nodal points.
Array<OneD, Array<OneD, NekDouble> > m_eta;
virtual NekVector<NekDouble> v_OrthoBasis(const int mode);
......@@ -131,6 +203,10 @@ protected:
}
};
/**
* @brief Specialisation of the NodalUtil class to support nodal tetrahedral
* elements.
*/
class NodalUtilTetrahedron : public NodalUtil
{
typedef boost::tuple<int, int, int> Mode;
......@@ -142,7 +218,12 @@ public:
Array<OneD, NekDouble> t);
protected:
/// Mapping from the \f$ (i,j,k) \f$ indexing of the basis to a continuous
/// ordering.
std::vector<Mode> m_ordering;
/// Collapsed coordinates \f$ (\eta_1, \eta_2, \eta_3) \f$ of the nodal
/// points.
Array<OneD, Array<OneD, NekDouble> > m_eta;
virtual NekVector<NekDouble> v_OrthoBasis(const int mode);
......@@ -167,6 +248,10 @@ protected:
}
};
/**
* @brief Specialisation of the NodalUtil class to support nodal prismatic
* elements.
*/
class NodalUtilPrism : public NodalUtil
{
typedef boost::tuple<int, int, int> Mode;
......@@ -178,7 +263,12 @@ public:
Array<OneD, NekDouble> t);
protected:
/// Mapping from the \f$ (i,j) \f$ indexing of the basis to a continuous
/// ordering.
std::vector<Mode> m_ordering;
/// Collapsed coordinates \f$ (\eta_1, \eta_2, \eta_3) \f$ of the nodal
/// points.
Array<OneD, Array<OneD, NekDouble> > m_eta;
virtual NekVector<NekDouble> v_OrthoBasis(const int mode);
......@@ -203,7 +293,7 @@ protected:
}
};
} // end of LibUtilities namespace
} // end of Nektar namespace
}
}
#endif //NODALUTIL_H
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment