node.h

Go to the documentation of this file.

#ifndef node_h
#define node_h
 
//#include <memory>
#include <iostream>
#include <eigen3/Eigen/Dense>
#include "config.h"
 
namespace Nodes
    {
    const int DIM = 3;
    
    const int NB_DATANODE = 2;
 
    enum step
        {
        CURRENT = 0,
        NEXT = 1
        };
 
    enum index
        {
        IDX_UNDEF = -1,
        IDX_X = 0,
        IDX_Y = 1,
        IDX_Z = 2
        };
 
inline double sq(const double x) { return x * x; }
 
struct dataNode
    {
    Eigen::Vector3d u;  
    Eigen::Vector3d v;  
    double phi;         
    double phiv;        
    };
 
struct Node
    {
    Eigen::Vector3d p;  
    Eigen::Vector3d ep; 
    Eigen::Vector3d eq; 
    dataNode d[NB_DATANODE];
 
    inline void setBasis(const double r)
        {
        // Choose for an initial ep the direction, among (X, Y, Z), which is further away from u0.
        Eigen::Index minIdx;
        d[CURRENT].u.cwiseAbs().minCoeff(&minIdx);
        #if EIGEN_VERSION_AT_LEAST(3,4,0)
            ep.setUnit(minIdx);
        #else
            if (minIdx == IDX_X) { ep = Eigen::Vector3d::UnitX(); }
            else if (minIdx == IDX_Y) { ep = Eigen::Vector3d::UnitY(); }
            else { ep = Eigen::Vector3d::UnitZ(); }
        #endif
        // Gram-Schmidt orthonormalization of (u0, ep).
        ep -= ep.dot(d[CURRENT].u) * d[CURRENT].u;
        ep.normalize();
 
        // Complete the basis with a vector product.
        eq = d[CURRENT].u.cross(ep);
 
        // Rotate (ep, eq) by the random angle.
        Eigen::Vector3d new_ep = cos(r) * ep - sin(r) * eq;
        eq = sin(r) * ep + cos(r) * eq;
        ep = new_ep;
 
        // The basis (u0, ep, eq) should already be orthonormal. An extra orthonormalization could
        // reduce the rounding errors.
        if (PARANOID_ORTHONORMALIZATION)
            {
            // Modified Gram-Schmidt orthonormalization.
            ep -= ep.dot(d[CURRENT].u) * d[CURRENT].u;
            ep.normalize();
            eq -= eq.dot(d[CURRENT].u) * d[CURRENT].u;
            eq -= eq.dot(ep) * ep;
            eq.normalize();
            }
        }
 
    inline void evolution(void) { d[step::CURRENT] = d[step::NEXT]; }
 
    inline void make_evol(const double vp , const double vq ,
                          const double dt )
        {
        d[NEXT].v = vp * ep + vq * eq;
        d[NEXT].u = d[CURRENT].u + dt * d[NEXT].v;
        d[NEXT].u.normalize();
        }
 
    inline double proj_ep(void) const { return d[NEXT].v.dot(ep); }
 
    inline double proj_eq(void) const { return d[NEXT].v.dot(eq); }
 
    inline const Eigen::Vector3d get_u(step k ) const
    { return d[k].u; }
 
    inline const Eigen::Vector3d get_v(step k ) const 
    { return d[k].v; }
 
    };  // end struct node
 
template <step K>
Eigen::Vector3d get_u(Node const &n ) { return n.d[K].u; }
 
template <step K>
Eigen::Vector3d get_v(Node const &n ) { return n.d[K].v; }
 
template <step K>
double get_phi(Node const &n ) { return n.d[K].phi; }
 
template <step K>
double get_phiv(Node const &n ) { return n.d[K].phiv; }
 
inline double get_u_comp(Node const &n , index idx )
    { return n.d[NEXT].u(idx); }
 
inline double get_v_comp(Node const &n , index idx )
    { return n.d[NEXT].v(idx); }
 
inline void set_phi(Node &n, double val) { n.d[NEXT].phi = val; }
 
inline void set_phiv(Node &n, double val) { n.d[NEXT].phiv = val; }
 
    }  // end namespace Nodes
 
#endif /* node_h */