Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Cleanup #5

Merged
merged 2 commits into from
May 12, 2017
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
147 changes: 135 additions & 12 deletions src/Basis1D.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -7,22 +7,42 @@

using namespace Eigen;

/**
* @brief Gives the value of pi. constexpr means it is evaluated at compile
* time.
*
* @return pi = std::atan(1)*4
*/
constexpr double pi() { return std::atan(1)*4; }

/**
* @brief Class for developing polynomial bases in 1D. Provides functions
* for evaluating the Legendre polynomials, which can be used to
* construct other bases.
*/
class Basis1D
{
public:
Basis1D() = default;
virtual ~Basis1D() = default;
protected:
// double evalLeg(double x, int n);
// double evalLegD(double x, int n);
RowVector2d evalLegLegD(double x, int n);
Vector2d evalLegm1Leg(double x, int n);
VectorXd evalLegendre(double x, int n);
MatrixX2d evalLegendreD(double x, int n);
};

/**
* @brief Evaluate the Legendre polynomials and their derivatives up to a
* certain degree.
*
* @param[in] x Coordinate value `-1 < x < 1`.
* @param[in] n Order of highest Legendre polynomial.
*
* @return `n+1`-by-`2` matrix with the first column holding the values of the
* Legendre polynomials and the second column holding the valuse of
* the derivatives of the Legendre polynomials.
*/
inline MatrixX2d Basis1D::evalLegendreD(double x, int n)
{
MatrixX2d p(n + 1, 2);
Expand All @@ -39,31 +59,52 @@ inline MatrixX2d Basis1D::evalLegendreD(double x, int n)
return p;
}

/**
* @brief Evaluate the Legendre polynomials up to a certain degree.
*
* @param[in] x Coordinate value `-1 < x < 1`.
* @param[in] n Order of highest Legendre polynomial.
*
* @return Vector of length `n+1` holding the values of the Legendre
* polynomials.
*/
inline VectorXd Basis1D::evalLegendre(double x, int n)
{
return evalLegendreD(x, n).leftCols<1>();
}

/**
* @brief Evaluate a Legendre polynomial of a certain degre and degree
* minus 1.
*
* @param[in] x Coordinate value `-1 < x < 1`.
* @param[in] n Order of highest Legendre polynomial.
*
* @return Vector of length `2` containing the two values.
*/
inline Vector2d Basis1D::evalLegm1Leg(double x, int n)
{
return evalLegendreD(x, n).bottomLeftCorner<2,1>();
}

/**
* @brief Evaluate a Legendre polynomial of a certain degree and its
* derivative.
*
* @param[in] x Coordinate value `-1 < x < 1`.
* @param[in] n Order of Legendre polynomial.
*
* @return Row vector of length `2` containing the two values.
*/
inline RowVector2d Basis1D::evalLegLegD(double x, int n)
{
return evalLegendreD(x, n).bottomRows<1>();
}

// inline double Basis1D::evalLeg(double x, int n)
// {
// return evalLegendreD(x, n)(n,0);
// }

// inline double Basis1D::evalLegD(double x, int n)
// {
// return evalLegendreD(x, n)(n,1);
// }

/**
* @brief Class for the Lagrange cardinal basis at the Gauss-Legendre (GL)
* points.
*/
class GaussLegendre : public Basis1D
{
public:
Expand All @@ -80,6 +121,17 @@ class GaussLegendre : public Basis1D
VectorXd node, weight;
};

/**
* @brief Constructs the object. During construction the GL nodes and
* weights are calculated. The values of the Legendre polynomials up
* to degree `k - 1` at each node are also stored in a matrix. To
* evaluate the cardinal functions, this matrix is inverted on the
* vector containing the evalutation of the Legendre polynomials at
* the desired point.
*
* @param[in] k The number of GL points to use.
* @param[in] eps The required accuracy of the GL points.
*/
GaussLegendre::GaussLegendre(int k, double eps) :
n_nodes(k), leg2lag(k,k), lu(k), node(k), weight(k)
{
Expand Down Expand Up @@ -107,26 +159,58 @@ GaussLegendre::GaussLegendre(int k, double eps) :
lu.compute(leg2lag);
}

/**
* @brief Evaluates the Lagrange GL cardinal fuctions and their derivatives
* at a desired location.
*
* @param[in] x Coordinate value `-1 < x < 1`.
*
* @return `k`-by-`2` matrix of the cardinal function values and their
* derivatives.
*/
inline MatrixX2d GaussLegendre::evalGL(double x)
{
return lu.solve(evalLegendreD(x, n_nodes - 1));
}

/**
* @brief Gets the node.
*
* @param[in] i Index of desired node.
*
* @return The node location.
*/
inline double GaussLegendre::getNode(int i)
{
return node[i];
}

/**
* @brief Gets the weight.
*
* @param[in] i Index of desired node.
*
* @return The weight.
*/
inline double GaussLegendre::getWeight(int i)
{
return weight[i];
}

/**
* @brief Gets the number of nodes.
*
* @return The number of nodes.
*/
inline int GaussLegendre::getN()
{
return n_nodes;
}

/**
* @brief Class for the Lagrange cardinal basis at the
* Gauss-Legendre-Lobatto (GLL) points.
*/
class GaussLobatto : public Basis1D
{
public:
Expand All @@ -143,6 +227,17 @@ class GaussLobatto : public Basis1D
VectorXd node, weight;
};

/**
* @brief Constructs the object. During construction the GLL nodes and
* weights are calculated. The values of the Legendre polynomials up
* to degree `k - 1` at each node are also stored in a matrix. To
* evaluate the cardinal functions, this matrix is inverted on the
* vector containing the evalutation of the Legendre polynomials at
* the desired point.
*
* @param[in] k The number of GLL points to use.
* @param[in] eps The required accuracy of the GLL points.
*/
GaussLobatto::GaussLobatto(int k, double eps) :
n_nodes(k), leg2lag(k,k), lu(k), node(k), weight(k)
{
Expand Down Expand Up @@ -177,21 +272,49 @@ GaussLobatto::GaussLobatto(int k, double eps) :
lu.compute(leg2lag);
}

/**
* @brief Evaluates the Lagrange GLL cardinal fuctions and their
* derivatives at a desired location.
*
* @param[in] x Coordinate value `-1 < x < 1`.
*
* @return `k`-by-`2` matrix of the cardinal function values and their
* derivatives.
*/
inline MatrixX2d GaussLobatto::evalGLL(double x)
{
return lu.solve(evalLegendreD(x, n_nodes - 1));
}

/**
* @brief Gets the node.
*
* @param[in] i Index of desired node.
*
* @return The node.
*/
inline double GaussLobatto::getNode(int i)
{
return node[i];
}

/**
* @brief Gets the weight.
*
* @param[in] i Index of desire node
*
* @return The weight.
*/
inline double GaussLobatto::getWeight(int i)
{
return weight[i];
}

/**
* @brief Gets the number of nodes.
*
* @return The number of nodes.
*/
inline int GaussLobatto::getN()
{
return n_nodes;
Expand Down
53 changes: 30 additions & 23 deletions src/ElementTransforms.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -6,29 +6,11 @@

using namespace Eigen;

// class Transform1D_Linear
// {
// public:
// Transform1D_Linear(double a = -1.0, double b = 1.0);
// ~Transform1D_Linear() = default;
// double forwardTransform(double x_hat);
// double jacobian();
// private:
// const RowVector2d x_range;
// };

// Transform1D_Linear::Transform1D_Linear(double a, double b) : x_range({a, b}) {}

// inline double Transform1D_Linear::forwardTransform(double x_hat)
// {
// return ((x_range[1] - x_range[0])*x_hat + x_range[0] + x_range[1])/2.0;
// }

// inline double Transform1D_Linear::jacobian()
// {
// return (x_range[1] - x_range[0])/2.0;
// }

/**
* @brief Class for one-dimensional transforms. Generates a coordinate
* transformation map between $[-1,1]$ and a general segment in $x$.
* This transformation can be of arbitrary polynomial order.
*/
class Transform1D
{
public:
Expand All @@ -42,19 +24,44 @@ class Transform1D
GaussLobatto gl;
};

/**
* @brief Constructs the object. Creates a linear transformation.
*
* @param[in] a Minimum $x$ value of the segment.
* @param[in] b Maximum $x$ value of the segment.
*/
Transform1D::Transform1D(double a, double b) : x_nodes(2), gl(2)
{
x_nodes[0] = a;
x_nodes[1] = b;
}

/**
* @brief Constructs the object. Creates an arbitrary transformation.
*
* @param[in] nodes The nodes corresponding to the transformed GLL points.
*/
Transform1D::Transform1D(RowVectorXd nodes) : x_nodes(nodes), gl(nodes.size()) {}

/**
* @brief Calculates the forward transformation $x=F(\hat{x})$.
*
* @param[in] x_hat $\hat{x}$
*
* @return $x$
*/
inline double Transform1D::forwardTransform(double x_hat)
{
return x_nodes*gl.evalGLL(x_hat).leftCols<1>();
}

/**
* @brief Calculates the Jacobian of the transform at $\hat{x}$.
*
* @param[in] x_hat \hat{x}$
*
* @return $J$
*/
inline double Transform1D::jacobian(double x_hat)
{
return x_nodes*gl.evalGLL(x_hat).rightCols<1>();
Expand Down
Loading