3
\$\begingroup\$

Introduction

I am developing a dynamical system simulation/animation framework, which provides three abstract base classes, DynamicalSystem, Visualisation and Animation, from which clients can derive concrete classes for the purposes of animating dynamical (physical) systems in 3D (and 2D).

The DynamicalSystem abstract class strips away all unnecessary details from a real physical system. In essence it blueprints a representation of the ordinary differential equation system for the physical system. It has two member variables: a dimension_ and a state_ vector. The pure virtual function stateDerivative(...) is to be overridden in the derived class whose duty it is to provide all the physical detail necessary to model the dynamical system.

The Visualisation abstract base class blueprints the interface for a graphical visualisation of the physical system (using OpenGL, typically). It has member variables representing a rendering target, and near and far clip plane distances (absolute values). The pure virtual functions display(...) and step(...) are to be overridden in the derived class whose duty it is to provide all the geometric detail necessary to model the dynamical system.

The Animation base class blueprints the interface for an animation. It contains pointers to abstract base classes DynamicalSystem and Visualisation, and a pointer to a GraphicsWindow object (I have written a GraphicsWindow class providing a thin veneer over SDL2. Code available on my GitHub -- link below.). It contains pure virtual functions handleKeys(...) and handleMouse(...) which the deriving class overrides to process user (window) events. It handles the deallocation (destruction) of the objects pointed to by the DynamicalSystem, Visualisation and GraphicsWindow base class pointers much like a smart triple pointer.

The code for these classes is provided below, as well as a concrete implementation for a simple molecular dynamics animation for example and completeness. The complete code with build instructions can be found at my GitHub https://github.com/Richard-Mace/MolecularDynamics.

Since this is perhaps my first attempt at "medium scale" object-oriented design, I would appreciate comments and suggestions for improvements. In particular, I know that the (derived?) animation object should be a singleton to avoid multiple instances. To my mind, the (abstract) base class Animation should ensure this behaviour by being a singleton itself. However, I cannot see an elegant way to achieve this through inheritance, and I would appreciate your opinions on how or whether this could be achieved, and whether it is worth it in terms of effort and potential code "uglification". (I know that singletons are a contentious issue. I am asking out of curiosity and to learn.)

The code is fairly long and, as already mentioned, can be downloaded from my GitHub at the link above (build instructions provided there). Below I have provided the code for the abstract base classes, then below that the code for an example molecular dynamics animation using these classes. I include the latter as an example of intended usage, for completeness. A full review would be an onerous task, but I would appreciate comments on the design of the base classes and how they might be improved.

Thank you in advance for your time.

Example main function

main.cpp

#include "GraphicsWindow.h"
#include "MoleculeSystem.h"
#include "MoleculeVisualisation.h"
#include "MoleculeAnimation.h"


int main() {
    MoleculeSystem::Extents  extents {-20.0, 20.0, -20.0, 20.0, -20.0, 20.0};
    MoleculeSystem::LoadType loadType   {MoleculeSystem::LoadType::SOLID};
    GraphicsWindow::Type     windowType {GraphicsWindow::Type::FULLSCREEN};

    MoleculeAnimation molecules(100, extents, loadType, windowType);

    molecules.run();

    return 0;
}

Abstract base classes (framework)

Animation

Animation.h

#ifndef ANIMATION_H
#define ANIMATION_H

#include "DynamicalSystem.h"
#include "Visualisation.h"

class Animation {
public:
    Animation(DynamicalSystem*, Visualisation*, GraphicsWindow*);
    virtual ~Animation();

    void setInitialTime(double);
    double getInitialTime() const;
    virtual void step(double, double);
    virtual void run();
    virtual bool handleKeys(const GraphicsWindow::Event&) = 0;
    virtual bool handleMouse(const GraphicsWindow::Event&) = 0;

protected:
    DynamicalSystem* ptrSystem_;
    Visualisation*   ptrVisualisation_;
    GraphicsWindow*  ptrWindow_;
    double           timeInitial_;
};

#endif /* ANIMATION_H */

Animation.cpp

#include "Animation.h"
#include <chrono>
#include <iostream>

Animation::Animation(DynamicalSystem* system, Visualisation* visualiser, 
        GraphicsWindow* window) 

    : ptrSystem_(system)
    , ptrVisualisation_(visualiser)
    , ptrWindow_(window)
    , timeInitial_(0.0) {

    // empty
}

// ~Animation:
//
// Destructor for the Animation class. 
//
// NOTE: The Animation object is responsible for the DELETION of the 
// GraphicsWindow object, the Visualisation object and the DynamicalSystem
// object. In effect is behaves like a smart triple pointer!
//

Animation::~Animation() {
    delete ptrWindow_;
    delete ptrVisualisation_;
    delete ptrSystem_;
}

/*
 * setInitialTime
 *
 * Sets the initial time, t_0,  for the animation.
 * 
 */

void Animation::setInitialTime(double time) {
    timeInitial_ = time;
}

/*
 * getInitialTime()
 *
 * Returns the initial time, t_0.
 * 
 */

double Animation::getInitialTime() const {
    return timeInitial_;
}

// step()
//
// The main sequence of calls in run(). Can be overridden in derived classes
// to provide additional functionality.

void Animation::step(double t1, double t2) {
    ptrVisualisation_->display(ptrSystem_);
    ptrSystem_->step(t1, t2);
    ptrVisualisation_->step(t1, t2);
}

/*
 * run()
 *
 * The main Animation loop.
 * 
 * Will return if the user closes the window, or if either handleKeys() or 
 * handleMouse() returns false.
 * 
 */

void Animation::run() {

    std::chrono::time_point<std::chrono::steady_clock> timeNow;
    std::chrono::time_point<std::chrono::steady_clock> timeStart;

    timeStart = timeNow = std::chrono::steady_clock::now();

    double timeLast = timeInitial_;

    while (true) {

        GraphicsWindow::Event event;
        while (ptrWindow_->pollEvent(&event)) {
            switch (event.type) {
                case SDL_QUIT:
                    return;                  // exit run loop
                    break;

                case SDL_KEYDOWN:
                case SDL_KEYUP:
                    if (!handleKeys(event)) { // return false to exit run loop
                        return;
                    }
                    break;

                case SDL_MOUSEMOTION:
                    if (!handleMouse(event)) { // return false to exit run loop
                        return;
                    }
                    break;

                default:
                    break;
            }
        }

        std::chrono::duration<double> elapsed_seconds = timeNow - timeStart;
        double time = elapsed_seconds.count() + timeInitial_;

        step(timeLast, time);

        timeLast = time;
        timeNow = std::chrono::steady_clock::now();
    }
}

Visualisation

Visualisation.h

    #ifndef VISUALISATION_H
    #define VISUALISATION_H

    #include "GraphicsWindow.h"
    #include "DynamicalSystem.h"

    class Visualisation {
    public:
        Visualisation(GraphicsWindow*, double, double);
        virtual ~Visualisation();

        virtual void display(const DynamicalSystem*) const = 0;
        virtual void step(double, double) = 0;
        virtual void setRenderTarget(GraphicsWindow*);
        virtual GraphicsWindow* getRenderTarget() const;

    protected:
        GraphicsWindow* ptrRenderTarget_;
        double          nearClip_;
        double          farClip_;

        // utility function
        bool init_opengl();
    };

    #endif /* VISUALISATION_H */

Visualisation.cpp

#include "Visualisation.h"
#include <GL/glu.h>
#include <iostream>
#include <stdexcept>


Visualisation::Visualisation(GraphicsWindow* ptrRenderTarget, double nearClip, 
                             double farClip) 

    : ptrRenderTarget_(ptrRenderTarget)
    , nearClip_(nearClip)
    , farClip_(farClip) {

    if (ptrRenderTarget_ != nullptr) {
        bool status = init_opengl();

        if (status == false) {
            throw std::runtime_error{"Visualisation: could not initialise OpenGL"};
        }
    }
}

Visualisation::~Visualisation() {
    ptrRenderTarget_ = nullptr;
    nearClip_ = farClip_ = 0.0;
}

GraphicsWindow* Visualisation::getRenderTarget() const {
    return ptrRenderTarget_;
}

void Visualisation::setRenderTarget(GraphicsWindow* ptrRenderTarget) {
    if (ptrRenderTarget != nullptr) {
        ptrRenderTarget_ = ptrRenderTarget;
        bool status = init_opengl();
        if (status == false) {
            throw std::runtime_error{"setRenderTarget: could not initialise OpenGL!\n"};
        }
    }
}

///////////////////////////////////////////////////////////////////////////////
// protected member functions                                                //
///////////////////////////////////////////////////////////////////////////////

bool Visualisation::init_opengl() {
    glClearColor(0.0f, 0.0f, 0.0f, 0.0f);
    glClearDepth(1.0f);
    glEnable(GL_DEPTH_TEST);
    glDepthFunc(GL_LESS);
    glHint(GL_PERSPECTIVE_CORRECTION_HINT, GL_NICEST);

    glViewport(0, 0, static_cast<GLsizei>(ptrRenderTarget_->widthPixels()), 
            static_cast<GLsizei>(ptrRenderTarget_->heightPixels()));

    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    gluPerspective(45.0f, 
                   static_cast<GLfloat>(ptrRenderTarget_->widthPixels()) / 
                       static_cast<GLfloat>(ptrRenderTarget_->heightPixels()), 
                   static_cast<GLfloat>(nearClip_), 
                   static_cast<GLfloat>(farClip_));

    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();

    glEnable(GL_LIGHTING);

    return true;
}

DynamicalSystem

DynamicalSystem.h

#ifndef DYNAMICALSYSTEM_H
#define DYNAMICALSYSTEM_H

#include <cstddef>
#include <vector>

class DynamicalSystem {
public:
    explicit DynamicalSystem(std::size_t);
    virtual ~DynamicalSystem() = default;

    std::size_t getDimension() const;
    const double* getStatePtr() const;
    virtual std::vector<double> 
        stateDerivative(double, std::vector<double>&) = 0;
    virtual void step(double, double);
    void setAccuracy(double);
    double getAccuracy() const;

protected:
    static constexpr double kDefaultAccuracyGoal_ = 1.0e-8; // integration 
                                                            // accuracy goal

    std::size_t         dimension_; // dimension (deg. of freedom) of system
    std::vector<double> state_;     // state vector (q_1,..,q_n,p_1,..,p_n)
    double              stepper_accuracy_; // integration accuracy

    // helper functions for system ODE integration
    void RungeKuttaCashKarp(double t, double deltat, 
                            const std::vector<double>& currDerivative, 
                            std::vector<double>& newState, 
                            std::vector<double>& stateError);

    void RungeKuttaQualityControlled(double t, 
                                     double deltaTry,
                                     double epsAcc, 
                                     const std::vector<double>& stateDerivative,
                                     const std::vector<double>& stateScale,
                                     double& deltaDid,
                                     double& deltaNext);

    void ODEIntegrate(double t1, double t2, double accuracy);
};

#endif /* DYNAMICALSYSTEM_H */

DynamicalSystem.cpp

#include "DynamicalSystem.h"
#include <iostream>
#include <cmath>

// Constructor
//

DynamicalSystem::DynamicalSystem(std::size_t dimension) 
    : dimension_(dimension)
    , stepper_accuracy_(kDefaultAccuracyGoal_) {
    state_.resize(2 * dimension_);
}

// Return the dimension of the dynamical system.
//
// The dimension is equivalent to the number of degrees of freedom of the 
// system. The order of the ordinary differential equation system corresponding
// to a DynamicalSystem is 2 * dimension.
//

std::size_t DynamicalSystem::getDimension() const {
    return dimension_;
}

// Return a read-only pointer to the beginning of the system state vector.

const double* DynamicalSystem::getStatePtr() const {
    return state_.data();
}

// Integrate the system differential equations from time t1 to time t2.

void DynamicalSystem::step(double t1, double t2) {
    ODEIntegrate(t1, t2, stepper_accuracy_);
}

// Set integration accuracy.

void DynamicalSystem::setAccuracy(double accuracy) {
    stepper_accuracy_ = accuracy;
}

// Get integration accuracy.

double DynamicalSystem::getAccuracy() const {
    return stepper_accuracy_;
}

///////////////////////////////////////////////////////////////////////////////
// private helper/utility member functions                                   //
///////////////////////////////////////////////////////////////////////////////

// RungeKuttaCashKarp:
//
// Given a vector currDerivative containing the time derivative of the system 
// state vector (state_) at time t, executes a low-level Runge-Kutta-Cash-Karp 
// step that reads the current member state_ (vector) of the dynamical system 
// at time t and returns in the vector argument newState the new state of the
// system at time t + deltat. The member variable state_ is NOT changed. An 
// estimate for the error in the solution is returned in argument stateError.
//
// This function uses the pure virtual member function stateDerivative(), which
// must be overriden in the derived class. The function 
// std::vector<double> stateDerivative(double t, std::vector<double>& state) 
// must return the time derivative at time t of an ARBITRARY given state vector, 
// state (at time t). The stateDerivative function effectively defines the 
// first order ordinary differential equation system governing the dynamics of
// the dynamical system (derived from abstract base class DynamicalSystem) being 
// modelled.
//
// REFERENCE: Numerical Recipes in C, p. 713.
//

void DynamicalSystem::RungeKuttaCashKarp(
                                    double t, 
                                    double deltat,
                                    const std::vector<double>& currDerivative,
                                    std::vector<double>& newState,
                                    std::vector<double>& stateError) {

    const double A2 = 0.2, A3 = 0.3, A4 = 0.6, A5 = 1.0,
                 A6 = 0.875, B21 = 0.2, B31 = 3.0 / 40.0,
                 B32 = 9.0 / 40.0, B41 = 0.3, B42 = -0.9,
                 B43 = 1.2, B51 = -11.0 / 54.0, B52 = 2.5,
                 B53 = -70.0 / 27.0, B54 = 35.0 / 27.0,
                 B61 = 1631.0 / 55296.0, B62 = 175.0 / 512.0,
                 B63 = 575.0 / 13824.0, B64 = 44275.0 / 110592.0,
                 B65 = 253.0 / 4096.0,
                 C1 = 37.0 / 378.0, C3 = 250.0 / 621.0,
                 C4 = 125.0 / 594.0, C6 = 512.0 / 1771.0,
                 DC1 = C1 - 2825.0 / 27648.0,
                 DC3 = C3 - 18575.0 / 48384.0,
                 DC4 = C4 - 13525.0 / 55296.0,
                 DC5 = -277.0 / 14336.0, DC6 = C6 - 0.25;

    std::size_t order = 2 * dimension_;

    // First step. We already have the derivs in currDerivative, so needn't
    // call stateDerivative here...
    std::vector<double> tempState(order);    
    for (std::size_t i = 0; i < order; ++i) {
        tempState[i] = state_[i] + B21 * deltat * currDerivative[i];
    }

    // Second step
    std::vector<double> ak2 = stateDerivative(t + A2 * deltat, tempState);
    for (std::size_t i = 0; i < order; ++i) {
        tempState[i] = state_[i] + deltat * (B31 * currDerivative[i] 
                + B32 * ak2[i]);
    }

    // Third step
    std::vector<double> ak3 = stateDerivative(t + A3 * deltat, tempState);
    for (std::size_t i = 0; i < order; ++i) {
        tempState[i] = state_[i] + deltat * (B41 * currDerivative[i] 
                + B42 * ak2[i] + B43 * ak3[i]);
    }

    // Fourth step
    std::vector<double> ak4 = stateDerivative(t + A4 * deltat, tempState);
    for (std::size_t i = 0; i < order; ++i) {
        tempState[i] = state_[i] + deltat * (B51 * currDerivative[i] 
                + B52 * ak2[i] + B53 * ak3[i] + B54 * ak4[i]);
    }

    // Fifth step
    std::vector<double> ak5 = stateDerivative(t + A5 * deltat, tempState);
    for (std::size_t i = 0; i < order; ++i) {
        tempState[i] = state_[i] + deltat * (B61 * currDerivative[i] 
                + B62 * ak2[i] + B63 * ak3[i] + B64 * ak4[i] + B65 * ak5[i]);
    }

    // Evaluate the new state
    std::vector<double> ak6 = stateDerivative(t + A6 * deltat, tempState);
    for (std::size_t i = 0; i < order; ++i) {
        newState[i] = state_[i] + deltat * (C1 * currDerivative[i] + C3 * ak3[i] 
                + C4 * ak4[i] + C6 * ak6[i]);
    }

    // The error is estimated as the difference between the
    // fourth and fifth order Runge-Kutta methods
    for (std::size_t i = 0; i < order; ++i) {
        stateError[i] = deltat * (DC1 * currDerivative[i] + DC3 * ak3[i] 
                + DC4 * ak4[i] + DC5 * ak5[i] + DC6 * ak6[i]);
    }
}


// Quality controlled Runge-Kutta-Cash-Karp step.
//
// Used by ODEIntegrate to integrate the ODE system defined in virtual member
// function stateDerivative(...) from time t to t + Delta, where Delta is 
// determined by accuracy and convergence conditions. The actual value of Delta
// used is returned in variable deltaDid, and a guess for Delta to be used in 
// the next integration step is returned in deltaNext. 
//
// This function updates the member variable state_, representing the current
// system state vector.

void DynamicalSystem::RungeKuttaQualityControlled(
                                    double t, 
                                    double deltaTry,
                                    double epsAcc, 
                                    const std::vector<double>& stateDerivative,
                                    const std::vector<double>& stateScale,
                                    double& deltaDid,
                                    double& deltaNext) {

    const double k_safety     = 0.9;
    const double k_pow_grow   = -0.2;
    const double k_pow_shrink = -0.25;
    const double k_error_cond = 1.89e-4;

    std::size_t order = 2 * dimension_;

    std::vector<double> stateError(order);
    std::vector<double> newState(order);

    double delta = deltaTry;
    double errorMax = 0.0;

    do {
        RungeKuttaCashKarp(t, delta, stateDerivative, newState, stateError);

        // Evaluate the error
        errorMax = 0.0;
        for (std::size_t i = 0; i < order; ++i) {
            errorMax = std::max(errorMax, fabs(stateError[i] / stateScale[i]));
        }

        errorMax /= epsAcc;

        if (errorMax > 1.0) {
            double deltaTemp = k_safety * delta * pow(errorMax, k_pow_shrink);
            delta = copysign(std::max(fabs(deltaTemp), 0.1 * fabs(delta)), delta);
            if (t + delta == t) {
                std::cerr <<"RungeKuttaQualityControlled: stepsize underflow!"
                          << '\n';
            }
        }
    } while (errorMax > 1.0);

    // We have a solution of the required accuracy. Set next stepsize.
    if (errorMax > k_error_cond) {
        deltaNext = k_safety * delta * pow(errorMax, k_pow_grow);
    } 
    else {
        deltaNext = 0.5 * delta;
    }

    // Transfer the new state vector to the member state_.
    deltaDid = delta;
    for (std::size_t i = 0; i < order; ++i) {
        state_[i] = newState[i];
    }
}


// High-level adaptive ode integrator. Integrates the ODE system defined in 
// virtual member function stateDerivative(...) from time t1 to time t2, where 
// usually t2 > t1. The member variable state_ representing the current system 
// state vector is updated by this function (through calls to
// RungeKuttaQualityControlled(...)).

void DynamicalSystem::ODEIntegrate(double t1, double t2, double accuracy) {

    const double k_tiny = 1.0e-30;

    std::size_t order = 2 * dimension_;

    double delta = 0.1;
    double t = t1;
    double deltaDid;
    double deltaNext;

    std::vector<double> stateDot(order);
    std::vector<double> stateScale(order);

    do {
        stateDot = stateDerivative(t, state_);

        for (std::size_t i = 0; i < order; ++i) {
            stateScale[i] = fabs(state_[i]) + fabs(delta * stateDot[i]) + k_tiny;
        }

        RungeKuttaQualityControlled(t, delta, accuracy, stateDot, stateScale, 
                deltaDid, deltaNext);

        t += deltaDid;
        delta = deltaNext;

        // make sure we don't overshoot
        if (t + delta > t2) {
            delta = t2 - t;
        }
    } while (t < t2);
}

Derived concrete classes (molecular dynamics example)

MoleculeAnimation

MoleculeAnimation.h

#ifndef MOLECULEANIMATION_H
#define MOLECULEANIMATION_H

#include "Animation.h"
#include "MoleculeSystem.h"
#include "MoleculeVisualisation.h"

class MoleculeAnimation : public Animation {
public:        
    MoleculeAnimation(std::size_t, const MoleculeSystem::Extents&, 
            MoleculeSystem::LoadType, GraphicsWindow::Type);

    ~MoleculeAnimation();

    virtual bool handleKeys(const GraphicsWindow::Event&) override;
    virtual bool handleMouse(const GraphicsWindow::Event&) override;

    // Clients of MoleculeSystem
    void expand();
    void contract();
    void heat();
    void cool();
    void toggleGravity();

    // Clients of MoleculeVisualisation
    void toggleCameraMotion();
    void toggleCameraDirection();
    void toggleCameraZoom();
    void changeColour();

private:
    static const std::size_t kDoF_                 = 3;
    static const std::size_t kDefaultWidth_        = 1920;
    static const std::size_t kDefaultHeight_       = 1080; 
    static constexpr double  kDefaultNearClip_     = 1.0;
    static constexpr double  kDefaultFarClip_      = 200.0;
    static constexpr double  kDefaultAccuracyGoal_ = 5.0e-5;
};

#endif /* MOLECULEANIMATION_H */

MoleculeAnimation.cpp

#include "MoleculeAnimation.h"
#include <iostream>
#include <stdexcept>


/*
 * MoleculeAnimation:
 *
 * Constructor for the MoleculeAnimation class. 
 *
 * NOTE: the deallocation of the MoleculeSystem, MoleculeVisualisation and
 *        GraphicsWindow objects dynamically allocated here is handled (smartly)
 *        by the Animation class destructor.
 */

MoleculeAnimation::MoleculeAnimation(std::size_t numMolecules, 
                                     const MoleculeSystem::Extents& extents,
                                     MoleculeSystem::LoadType loadType,
                                     GraphicsWindow::Type windowType) 

    : Animation(new MoleculeSystem(numMolecules, kDoF_, extents, loadType),
                nullptr, 
                new GraphicsWindow) {

    if (!ptrWindow_->isInitialised()) {
        throw std::runtime_error{"MoleculeAnimation: could not init window"};
    }

    ptrWindow_->setTitle("Molecular Dynamics Animation");

    bool success = true;
    if (windowType == GraphicsWindow::Type::WINDOWED) {
        success = ptrWindow_->open(kDefaultWidth_, kDefaultHeight_);
    }
    else {
        success = ptrWindow_->open();
    }
    if (!success) {
        throw std::runtime_error{"MoleculeAnimation: could not open window"};
    }

    ptrSystem_->setAccuracy(kDefaultAccuracyGoal_);

    // Create a molecule visualisation for the system. This must be done
    // only after a GraphicsWindow has been successfully opened.
    ptrVisualisation_ = new MoleculeVisualisation(ptrWindow_, kDefaultNearClip_, 
            kDefaultFarClip_);
}

MoleculeAnimation::~MoleculeAnimation() {
    std::cout << "MoleculeAnimation: Destructor." << '\n';
}

/*
 * handleKeys:
 *
 * The keyboard event processing handler.
 *
 */
bool MoleculeAnimation::handleKeys(const GraphicsWindow::Event& event) {

    if (event.type == SDL_KEYDOWN) {
        switch (event.key.keysym.sym) {
            case SDLK_ESCAPE:
                return false;
                break;

            case SDLK_UP:
                heat();
                break;

            case SDLK_DOWN:
                cool();
                break;

            case SDLK_LEFT:
                contract();
                break;

            case SDLK_RIGHT:
                expand();
                break;

            case SDLK_g:
                toggleGravity();
                break;

            case SDLK_p:
                toggleCameraMotion();
                break;

            case SDLK_z:
                toggleCameraZoom();
                break;

            case SDLK_r:
                toggleCameraDirection();
                break;

            case SDLK_c:
                changeColour();
                break;

            default:
                break;
        }
    }

    return true;
}

/*
 * handleMouse:
 *
 * The mouse event processing handler. Presently not used.
 *
 */

bool MoleculeAnimation::handleMouse(const GraphicsWindow::Event& event) {
    return true;
}


///////////////////////////////////////////////////////////////////////////////
// Molecule System (physical) effects                                        //
///////////////////////////////////////////////////////////////////////////////

/*
 * expand()
 *
 * Expand the system container (indefinitely) without volume checks.
 * 
 */

void MoleculeAnimation::expand() {
    MoleculeSystem::Extents extents = 
                        dynamic_cast<MoleculeSystem*>(ptrSystem_)->getExtents();
    double speed = 0.1;

    extents.xmin -= speed * MoleculeVisualisation::kMoleculeRadius_;
    extents.xmax += speed * MoleculeVisualisation::kMoleculeRadius_;
    extents.ymin -= speed * MoleculeVisualisation::kMoleculeRadius_;
    extents.ymax += speed * MoleculeVisualisation::kMoleculeRadius_;
    extents.zmin -= speed * MoleculeVisualisation::kMoleculeRadius_;
    extents.zmax += speed * MoleculeVisualisation::kMoleculeRadius_;

    dynamic_cast<MoleculeSystem*>(ptrSystem_)->setExtents(extents);
}

/*
 * contract()
 *
 * Contract the system container until the volume equals the number of molecules
 * times the volume of a cube with side length equal to the molecule diameter.
 *
 * TODO: Small volumes create stress on the integration routines, whose accuracy
 * has been minimised for speed. Think about how better to deal with low volumes
 * and high energy density. Perhaps limit the contraction by keeping track of
 * total energy.   
 */
void MoleculeAnimation::contract() {

    MoleculeSystem* ptrMolySys = dynamic_cast<MoleculeSystem*>(ptrSystem_);

    MoleculeSystem::Extents extents = ptrMolySys->getExtents();

    double volume = (extents.xmax - extents.xmin) 
                    * (extents.ymax - extents.ymin)
                    * (extents.zmax - extents.zmin);

    double r = MoleculeVisualisation::kMoleculeRadius_;

    double num = static_cast<double>(ptrMolySys->getNumber());

    // Minimum volume corresponds to each molecule confined to a cube
    // of side length equal to 2 * r.
    double volume_min = 8.0 * r * r * r * num;

    // Don't allow unlimited squashing.
    if (volume <= volume_min) {
        return;
    } 

    double speed = -0.1;

    extents.xmin -= speed * MoleculeVisualisation::kMoleculeRadius_;
    extents.xmax += speed * MoleculeVisualisation::kMoleculeRadius_;
    extents.ymin -= speed * MoleculeVisualisation::kMoleculeRadius_;
    extents.ymax += speed * MoleculeVisualisation::kMoleculeRadius_;
    extents.zmin -= speed * MoleculeVisualisation::kMoleculeRadius_;
    extents.zmax += speed * MoleculeVisualisation::kMoleculeRadius_;

    ptrMolySys->setExtents(extents);
}

/*
 * toggleGravity()
 *
 * Alternatively turn on and turn off the gravitational field of an external
 * (planetary, say) system.
 */

void MoleculeAnimation::toggleGravity() {
    MoleculeSystem* ptrMolySys = dynamic_cast<MoleculeSystem*>(ptrSystem_);

    if (ptrMolySys->getGravitationalAcceleration() == 0.0) {
        ptrMolySys->setGravitationalAcceleration(0.01);
    }
    else {
        ptrMolySys->setGravitationalAcceleration(0.0);
    }
}

/*
 * heat()
 *
 * Simulate the effects of heating the molecule system, simply by scaling
 * the velocity of each molecule by a pre-determined factor.
 * 
 * This is a bit artificial and should be replaced by a process in which each
 * molecule's velocity is perturbed by an amount determined by deviates from a 
 * 3D Gaussian distribution with a given temperature.
 */

void MoleculeAnimation::heat() {
    dynamic_cast<MoleculeSystem*>(ptrSystem_)->scaleVelocities(1.05);
}

/*
 * cool()
 *
 * Simulate the effects of cooling the molecule system, simply by scaling 
 * the velocity of each molecule by a pre-determined factor. See comments for
 * member function heat().
 */

void MoleculeAnimation::cool() {
    dynamic_cast<MoleculeSystem*>(ptrSystem_)->scaleVelocities(0.95);
}


///////////////////////////////////////////////////////////////////////////////
// Visualisation effects                                                     //
///////////////////////////////////////////////////////////////////////////////

/*
 * toggleCameraMotion:
 *
 * Alternatively stop or start the camera orbit.
 */

void MoleculeAnimation::toggleCameraMotion() {
    dynamic_cast<MoleculeVisualisation*>(ptrVisualisation_)->toggleCameraMotion();
}

/*
 * toggleCameraDirection()
 *
 * Reverses the direction of motion of the camera in its orbit.
 *
 */

void MoleculeAnimation::toggleCameraDirection() {
    dynamic_cast<MoleculeVisualisation*>(ptrVisualisation_)->toggleCameraDirection();
}

/*
 * toggleCameraZoom()
 *
 * Alternatively zooms in and out of the scene.
 *
 */

void MoleculeAnimation::toggleCameraZoom() {
    dynamic_cast<MoleculeVisualisation*>(ptrVisualisation_)->toggleCameraZoom();
}

/*
 * changeColour()
 *
 * Cycles through a stack of pre-set materials (colours).
 *
 */

void MoleculeAnimation::changeColour() {
    MoleculeVisualisation *ptrMolyVis = 
        dynamic_cast<MoleculeVisualisation*>(ptrVisualisation_);

    std::size_t num = ptrMolyVis->getNumMaterials();
    std::size_t currIndex = ptrMolyVis->getMaterialIndex();
    std::size_t newIndex = (currIndex + 1) % num;
    ptrMolyVis->setMaterialIndex(newIndex);
}

MoleculeVisualisation

MoleculeVisualisation.h

#ifndef MOLECULEVISUALISATION_H
#define MOLECULEVISUALISATION_H

#include "GraphicsWindow.h"
#include "Visualisation.h"
#include "MoleculeSystem.h"

class MoleculeVisualisation : public Visualisation {
public:
    struct Light {
        GLfloat ambient_[4];
        GLfloat diffuse_[4];
        GLfloat specular_[4];
        GLfloat position_[4];
    };

    struct Material {
        GLfloat ambient_[4];   
        GLfloat diffuse_[4];   
        GLfloat specular_[4];  
        GLfloat shininess_;    
    };

    MoleculeVisualisation(GraphicsWindow*, double, double);
    virtual ~MoleculeVisualisation() override;

    void setFocalPoint(double, double, double);
    void getFocalPoint(double*, double*, double*) const;
    void setAmbientLight(GLfloat, GLfloat, GLfloat);
    void addLight(const Light&);
    void setMaterialIndex(std::size_t);
    std::size_t getMaterialIndex() const;
    std::size_t getNumMaterials() const;
    virtual void display(const DynamicalSystem*) const override;
    virtual void step(double, double) override;
    void toggleCameraMotion();
    void toggleCameraZoom();
    void toggleCameraDirection();

    static constexpr double kMoleculeRadius_ = 0.5; // Radius -- don't change. 

private:
    static const int        kSphereSegments_ = 20;  // For sphere tessellation.
    static const int        kNumLines_       = 60;  // For the surface plane.
    static constexpr double kZoomFactor_     = 0.3; // Simulates camera zoom.

    bool                    cameraPaused_;
    double                  cameraDirection_;
    double                  cameraZoomFactor_;
    double                  cameraPosition_[3];
    double                  focalPoint_[3]; 
    std::vector<Material>   materials_;
    std::size_t             materialIndex_;
    unsigned int            numLights_;  

    // utility functions
    void draw_box(const MoleculeSystem::Extents&) const;
    void orbit_camera(double, double);
    void init_material_buffer();
};

#endif /* MOLECULEVISUALISATION_H */

MoleculeVisualisation.cpp

#include "GraphicsWindow.h"
#include "MoleculeVisualisation.h"
#include <iostream>
#include <GL/gl.h>
#include <GL/glu.h>

MoleculeVisualisation::MoleculeVisualisation(GraphicsWindow* ptrRenderTarget,
                                             double nearClip, double farClip) 

    : Visualisation(ptrRenderTarget, nearClip, farClip)
    , cameraPaused_(false)
    , cameraDirection_(1.0)
    , cameraZoomFactor_(1.0)
    , cameraPosition_{0.0, 18.0, 75.0}
    , focalPoint_{0.0, 0.0, 0.0}
    , materialIndex_(0)
    , numLights_(0) {

    setAmbientLight(0.4f, 0.4f, 0.4f);

    Light default_light { {0.5f, 0.5f, 0.5f, 1.0f},
                          {0.6f, 0.6f, 0.6f, 1.0f},
                          {0.6f, 0.6f, 0.6f, 1.0f},
                          {0.0f, 40.0f, 20.0f, 1.0f} };

    addLight(default_light);

    Light second_light { {0.5f, 0.5f, 0.5f, 1.0f},
                         {0.6f, 0.6f, 0.6f, 1.0f},
                         {0.6f, 0.6f, 0.6f, 1.0f},
                         {0.0f, 40.0f, -20.0f, 1.0f} };

    addLight(second_light);

    init_material_buffer();
}


MoleculeVisualisation::~MoleculeVisualisation() {    
    numLights_ = 0;

    std::cout << "MoleculeVisualisation: Destructor." << '\n';
}

void MoleculeVisualisation::setFocalPoint(double x, double y, double z) {
    focalPoint_[0] = x;
    focalPoint_[1] = y;
    focalPoint_[2] = z;
}

void MoleculeVisualisation::getFocalPoint(double* x, double* y, 
                                          double* z) const {
    *x = focalPoint_[0];
    *y = focalPoint_[1];
    *z = focalPoint_[2];
}

void MoleculeVisualisation::setAmbientLight(GLfloat r, GLfloat g, GLfloat b) {
    GLfloat ambient[] = {r, g, b, 1.0f};
    glLightModelfv(GL_LIGHT_MODEL_AMBIENT, ambient);
}

void MoleculeVisualisation::addLight(const Light& light) {
    bool fail = false;
    unsigned int flag = 0;

    switch (numLights_) {
        case 0:
            flag = GL_LIGHT0;
            break;

        case 1: 
            flag = GL_LIGHT1;
            break;

        case 2:
            flag = GL_LIGHT2;
            break;

        case 3:
            flag = GL_LIGHT3;
            break;

        case 4:
            flag = GL_LIGHT4;
            break;

        case 5:
            flag = GL_LIGHT5;
            break;

        case 6:
            flag = GL_LIGHT6;
            break;

        case 7:
            flag = GL_LIGHT7;
            break;

        default:
            std::cerr <<"MoleculeVisualisation: Cannot add light. Too many lights" 
                      << std::endl;
            fail = true;
            break;
    }

    if (!fail) {
        glLightfv(flag, GL_AMBIENT,  light.ambient_);
        glLightfv(flag, GL_DIFFUSE,  light.diffuse_);
        glLightfv(flag, GL_SPECULAR, light.specular_);
        glLightfv(flag, GL_POSITION, light.position_);
        glEnable(flag);

        ++numLights_;
    }
}

void MoleculeVisualisation::setMaterialIndex(std::size_t index) {
    if (index <= materials_.size()) {
        materialIndex_ = index;
    }
}

std::size_t MoleculeVisualisation::getMaterialIndex() const {
    return materialIndex_;
}

std::size_t MoleculeVisualisation::getNumMaterials() const {
    return materials_.size();
}

void MoleculeVisualisation::display(const DynamicalSystem* ptrSystem) const {

    // Clear (fill) the frame and depth (z-) buffers.
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

    glLoadIdentity();

    // Follow a molecule for fun. A bit jiggly. Experiment with it.
    //const double* ptrState = ptrSystem->getStatePtr();
    //setFocalPoint(ptrState[0], ptrState[1], ptrState[2]);

    // Set up eye (coordinate system) transformations.
    gluLookAt(cameraPosition_[0], cameraPosition_[1], cameraPosition_[2], 
              focalPoint_[0], focalPoint_[1], focalPoint_[2], 0.0, 1.0, 0.0);

    // Draw the confining box.
    MoleculeSystem::Extents extents = 
            dynamic_cast<const MoleculeSystem*>(ptrSystem)->getExtents();

    draw_box(extents);

    // Set the material properties for all molecules.
    glMaterialfv(GL_FRONT, GL_SPECULAR, materials_[materialIndex_].specular_);
    glMaterialfv(GL_FRONT, GL_AMBIENT,  materials_[materialIndex_].ambient_);
    glMaterialfv(GL_FRONT, GL_DIFFUSE,  materials_[materialIndex_].diffuse_);
    glMaterialf(GL_FRONT, GL_SHININESS, materials_[materialIndex_].shininess_);

    const double* ptrMolecule = ptrSystem->getStatePtr();

    const size_t degreesOfFreedom = 
        dynamic_cast<const MoleculeSystem*>(ptrSystem)->getDegreesOfFreedom();

    GLUquadric* ptr_sphere = gluNewQuadric();

    const size_t numMolecules = 
        dynamic_cast<const MoleculeSystem*>(ptrSystem)->getNumber();

    for (size_t i = 0; i < numMolecules; ++i) {
        GLfloat x = static_cast<GLfloat>(ptrMolecule[0]);
        GLfloat y = static_cast<GLfloat>(ptrMolecule[1]);
        GLfloat z = static_cast<GLfloat>(ptrMolecule[2]); 

        // Position and draw molecule.
        glPushMatrix();
        glTranslatef(y, z, x);
        gluSphere(ptr_sphere, kMoleculeRadius_, kSphereSegments_, 
                kSphereSegments_);
        glPopMatrix();

        ptrMolecule += degreesOfFreedom;
    }

    gluDeleteQuadric(ptr_sphere);

    // Swap backbuffer with frontbuffer.
    ptrRenderTarget_->swapBuffers();
}

void MoleculeVisualisation::step(double t1, double t2) {
    orbit_camera(t1, t2);
}

void MoleculeVisualisation::toggleCameraMotion() {
    cameraPaused_ = !cameraPaused_;
}

void MoleculeVisualisation::toggleCameraZoom() {
    if (cameraZoomFactor_ == 1.0) {
        cameraZoomFactor_ = kZoomFactor_;
    }
    else {
        cameraZoomFactor_ = 1.0;
    }
}

void MoleculeVisualisation::toggleCameraDirection() {
    if (cameraDirection_ == 1.0) {
        cameraDirection_ = -1.0;
    }
    else {
        cameraDirection_ = 1.0;
    }
}

///////////////////////////////////////////////////////////////////////////////
// private utility functions                                                 //
///////////////////////////////////////////////////////////////////////////////

void MoleculeVisualisation::draw_box(
                                const MoleculeSystem::Extents& extents) const {

    GLfloat spec1[] = {0.9f, 0.9f, 0.9f, 1.0f};
    GLfloat amb1[] = {1.0f, 1.0f, 1.0f, 1.0f};
    GLfloat diff1[] = {1.0f, 1.0f, 1.0f, 1.0f};
    GLfloat shin1[] = {32.0f};

    glMaterialfv(GL_FRONT, GL_SPECULAR, spec1);
    glMaterialfv(GL_FRONT, GL_SHININESS, shin1);
    glMaterialfv(GL_FRONT, GL_AMBIENT, amb1);
    glMaterialfv(GL_FRONT, GL_DIFFUSE, diff1);

    // Account for the finite size of the "point" molecules
    GLfloat xmin = static_cast<GLfloat>(extents.xmin) - kMoleculeRadius_;
    GLfloat xmax = static_cast<GLfloat>(extents.xmax) + kMoleculeRadius_;
    GLfloat ymin = static_cast<GLfloat>(extents.ymin) - kMoleculeRadius_;
    GLfloat ymax = static_cast<GLfloat>(extents.ymax) + kMoleculeRadius_;
    GLfloat zmin = static_cast<GLfloat>(extents.zmin) - kMoleculeRadius_;
    GLfloat zmax = static_cast<GLfloat>(extents.zmax) + kMoleculeRadius_;

    glBegin(GL_LINES);

    glVertex3f(xmin, ymin, zmin);
    glVertex3f(xmax, ymin, zmin);
    glVertex3f(xmax, ymin, zmin);
    glVertex3f(xmax, ymax, zmin);
    glVertex3f(xmax, ymax, zmin);
    glVertex3f(xmin, ymax, zmin);
    glVertex3f(xmin, ymax, zmin);
    glVertex3f(xmin, ymin, zmin);

    glVertex3f(xmin, ymin, zmax);
    glVertex3f(xmax, ymin, zmax);
    glVertex3f(xmax, ymin, zmax);
    glVertex3f(xmax, ymax, zmax);
    glVertex3f(xmax, ymax, zmax);
    glVertex3f(xmin, ymax, zmax);
    glVertex3f(xmin, ymax, zmax);
    glVertex3f(xmin, ymin, zmax);

    glVertex3f(xmin, ymin, zmin);
    glVertex3f(xmin, ymin, zmax);

    glVertex3f(xmax, ymin, zmin);
    glVertex3f(xmax, ymin, zmax);

    glVertex3f(xmax, ymax, zmin);
    glVertex3f(xmax, ymax, zmax);

    glVertex3f(xmin, ymax, zmin);
    glVertex3f(xmin, ymax, zmax);

    // draw floor
    GLfloat deltax = 20.5f / 4.0f;
    GLfloat deltay = 20.5f / 4.0f;

    GLfloat col_green[] {0.0f, 0.25f, 0.0f, 1.0f};

    glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, col_green);
    for (int i = 0; i < kNumLines_; ++i) {
        glVertex3f(i * deltax, zmin - 0.1f, -kNumLines_ * deltay);
        glVertex3f(i * deltax, zmin - 0.1f,  kNumLines_ * deltay);

        glVertex3f(-kNumLines_ * deltax, zmin - 0.1f, i * deltay);
        glVertex3f(kNumLines_ * deltax, zmin - 0.1f,  i * deltay);
    }

    for (int i = 1; i < kNumLines_; ++i) {
        glVertex3f(-i * deltax, zmin - 0.1f, -kNumLines_ * deltay);
        glVertex3f(-i * deltax, zmin - 0.1f,  kNumLines_ * deltay);

        glVertex3f(-kNumLines_ * deltax, zmin - 0.1f, -i * deltay);
        glVertex3f(kNumLines_ * deltax, zmin - 0.1f,  -i * deltay);
    }

    glEnd();
}

void MoleculeVisualisation::orbit_camera(double t1, double t2) {

    const double k_cam_max_r = 50.0;
    const double k_cam_min_r = 30.0;
    const double k_cam_max_z = 30.0;
    const double k_slow = 0.5;

    static double cam_time = 0.0;


    double delta_t;
    if (cameraPaused_) {
        delta_t = 0.0;
    }
    else {
        delta_t = t2 - t1;
    }

    if (cameraDirection_ == -1.0) {
        delta_t *= -1.0;
    }

    cam_time += delta_t;

    double r = k_cam_max_r + (k_cam_max_r - k_cam_min_r) * 
        cos(0.1 * cam_time * k_slow);
    cameraPosition_[0] = cameraZoomFactor_ * r * cos(0.6 * cam_time * k_slow);
    cameraPosition_[2] = cameraZoomFactor_ * r * sin(0.3 * cam_time * k_slow);
    cameraPosition_[1] = cameraZoomFactor_ * k_cam_max_z * 
            (sin(0.5 * cam_time * k_slow) + 1.02);
}

/*
 * init_material_buffer:
 *
 * Defines a few pre-set materials for the user to cycle through, if desired.
 *
 */

void MoleculeVisualisation::init_material_buffer() {

    Material def_mat { {0.5f, 0.1995f, 0.0745f, 1.0f},
                       {0.75164f, 0.60648f, 0.22648f, 1.0f},
                       {0.9f, 0.9f, 0.9f, 1.0f},
                        32.0f };

    materials_.push_back(def_mat);

    Material emerald { {0.0215f, 0.1745f, 0.0215f, 1.0f},
                       {0.07568f, 0.61424f, 0.07568f, 1.0f},
                       {0.633f, 0.727811f, 0.633f, 1.0f}, 
                        76.8f };

    materials_.push_back(emerald);

    Material obsidian { {0.05375f, 0.05f, 0.06625f, 1.0f},
                        {0.18275f, 0.17f, 0.22525f, 1.0f},
                        {0.332741f, 0.328634f, 0.346435f, 1.0f}, 
                         38.4f };

    materials_.push_back(obsidian);

    Material ruby { {0.1745f, 0.01175f, 0.01175f, 1.0f},
                    {0.61424f, 0.04136f, 0.04136f, 1.0f},
                    {0.727811f, 0.626959f, 0.626959f, 1.0f},
                     76.8f };

    materials_.push_back(ruby);

    Material chrome { { 0.25f, 0.25f, 0.25f, 1.0f}, 
                      {0.4f, 0.4f, 0.4f, 1.0f},
                      {0.774597f, 0.774597f, 0.774597f, 1.0f},
                       76.8f };

    materials_.push_back(chrome);

    Material cyan_plastic{ {0.0f, 0.1f, 0.06f, 1.0f},
                           {0.0f, 0.50980392f, 0.50980392f, 1.0f},
                           {0.50196078f, 0.50196078f, 0.50196078f, 1.0f},
                            32.0f };

    materials_.push_back(cyan_plastic);    

    Material copper { {0.19125f, 0.0735f, 0.0225f, 1.0f},
                      {0.7038f, 0.27048f, 0.0828f, 1.0f},
                      {0.25677f, 0.137622f, 0.086014f, 1.0f}, 
                       12.8f };

    materials_.push_back(copper);   
}

MoleculeSystem

MoleculeSystem.h

#ifndef MOLECULESYSTEM_H
#define MOLECULESYSTEM_H

#include "DynamicalSystem.h"

class MoleculeSystem : public DynamicalSystem {
public:        
    struct Extents {
        double xmin;
        double xmax;
        double ymin;
        double ymax;
        double zmin;
        double zmax;
    };

    enum class LoadType {SOLID, LIQUID, GAS};

    MoleculeSystem(std::size_t, std::size_t, const Extents&, LoadType);
    ~MoleculeSystem();

    virtual std::vector<double> 
        stateDerivative(double, std::vector<double>&) override;
    void step(double t1, double t2) override;
    void setGravitationalAcceleration(double);
    double getGravitationalAcceleration() const;
    std::size_t getNumber() const; 
    std::size_t getDegreesOfFreedom() const;
    Extents getExtents() const;
    void setExtents(const MoleculeSystem::Extents&);
    void scaleVelocities(double factor);

private:
    static constexpr double kKineticMax_                = 10.0;
    static constexpr double kCoordinationNumber_        = 2.0;
    static constexpr double kKToverEpsilonGas_          = 2.0;
    static constexpr double kKToverEpsilonLiquid_       = 0.8;
    static constexpr double kKToverEpsilonSolid_        = 1.0e-5;

    std::size_t numMolecules_;        // number of molecules
    std::size_t numDegreesOfFreedom_; // number of DoF per molecule
    Extents     extents_;        // container extents in x-, y-, z-directions
    double      gravitationalAcceleration_;  // external gravitational field g

    // private helper functions
    void initRandom(double);
    void initLatticeHCP(double);
    void applyBoundaryConditions();
};

#endif /* MOLECULESYSTEM_H */

MoleculeSystem.cpp (see GitHub)

\$\endgroup\$

1 Answer 1

3
\$\begingroup\$

Ownership

This:

// NOTE: The Animation object is responsible for the DELETION of the 
// GraphicsWindow object, the Visualisation object and the DynamicalSystem
// object. In effect is behaves like a smart triple pointer!
//

Animation::~Animation() {
    delete ptrWindow_;
    delete ptrVisualisation_;
    delete ptrSystem_;

sets off alarm bells. There are new methods for C++ to safely transfer ownership of a pointer; it used to be auto_ptr and now unique_ptr is recommended. Among other things, it will greatly simplify your destructor. Further reading: https://en.cppreference.com/w/cpp/memory/unique_ptr

Even if you didn't use an autopointer, these:

DynamicalSystem* ptrSystem_;
Visualisation*   ptrVisualisation_;
GraphicsWindow*  ptrWindow_;

can be made const (the pointer itself does not change), i.e.

DynamicalSystem *const ptrSystem_;
Visualisation *const ptrVisualisation_;
GraphicsWindow *const ptrWindow_;

Boolean expressions

if (status == false) {

should simply be

if (!status) {

Structured coordinates

This:

const double* ptrMolecule = ptrSystem->getStatePtr();
// ...

    GLfloat x = static_cast<GLfloat>(ptrMolecule[0]);
    GLfloat y = static_cast<GLfloat>(ptrMolecule[1]);
    GLfloat z = static_cast<GLfloat>(ptrMolecule[2]); 

seems like a misuse of unsafe pointer references. The moment that you return a pointer, you lose information about the size of the array. In this particular case, if you were to keep it as a const reference to a vector, you could still index into it - but the indexing would have bounds-checking safety.

Templating

You may want to consider making the degree of freedom of DynamicalSystem an integer template parameter. That way, children such as MoleculeSystem can quickly and simply declare the DOF in the inheritance clause, and you can use a simple array instead of a vector.

\$\endgroup\$
3
  • \$\begingroup\$ Thank you for your detailed review. Regarding ownership, I am aware of unique_ptr and RAII, but I think that my Animation constructor and Destructor work, in principle, precisely as a unique_ptr. In effect, the Animation object is a wrapper around three base class pointers (owning resources) such that, when the Animation object (the wrapper) goes out of scope, the resources are automatically deleted by invoking the destructor of the type to which the enclosed pointers point. \$\endgroup\$ Mar 27, 2020 at 9:13
  • \$\begingroup\$ I stand to be corrected, but (polymorphism aside) that is my understanding of how unique pointers are designed and work, perhaps with a bit more sophistication in the case of std::unique_ptr. I agree that const should have been used. Thank you for pointing this out. \$\endgroup\$ Mar 27, 2020 at 9:14
  • \$\begingroup\$ Regarding the returning of a const double* instead of const vector<double>&, I agree. My oversight; I shall correct this in the next version. Regarding templating and use of std::array, I think that this is a good idea. I shall consider the implications in the next version. Thank you once again! \$\endgroup\$ Mar 27, 2020 at 9:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.