Add a visualization

README.md

@@ -7,6 +7,7 @@ Physical simulation experiments
 - a modern C++ compiler (i.e. one that supports C++17 features, ex. gcc 7+)
 - cmake 3.13 or newer
 - Eigen3 (installable via apt, yum, or homebrew)
+- GLUT and OpenGL
 
 ## Building
 

simucene/CMakeLists.txt

+find_package(OpenGL REQUIRED)
+find_package(GLUT REQUIRED)
+
 add_library(simucene_lib
+    glut_grapher.cpp
+    glut_grapher.h
     matrix_builder.cpp
     matrix_builder.h
     vector.h
@@ -10,4 +15,8 @@ target_compile_features(simucene_lib PUBLIC cxx_std_17)
 add_executable(simucene
     main.cpp
 )
-target_link_libraries(simucene shared_settings simucene_lib)
+target_link_libraries(simucene shared_settings simucene_lib ${OPENGL_LIBRARIES} ${GLUT_LIBRARY})
+target_include_directories(simucene PUBLIC
+    ${OPENGL_INCLUDE_DIRS}
+    ${GLUT_INCLUDE_DIRS}
+)

simucene/bounding_box.h

+#ifndef SIMUCENE_BOUNDING_BOX_H
+#define SIMUCENE_BOUNDING_BOX_H
+
+#include <iostream>
+#include <stdint.h>
+
+namespace simucene {
+
+//--------------------------------------------------------------------------------------------------
+// An axis aligned bounding box. It stores a low corner and a high corner, each represented as a T.
+// The idea is that this makes BoundingBox dimension and representation agnostic, though the class
+// is designed with fixed size Eigen vectors in mind for the latter.
+template <typename T>
+class BoundingBox {
+public:
+    BoundingBox() = default;
+    BoundingBox(BoundingBox const&) = default;
+    BoundingBox(BoundingBox&&) = default;
+    BoundingBox(T const& x): low_(x), high_(x) {}
+    BoundingBox(T const& low, T const& high): low_(low), high_(high) {}
+
+    BoundingBox& operator=(BoundingBox const&) = default;
+    BoundingBox& operator=(BoundingBox&&) = default;
+    BoundingBox& operator=(T const& x);
+
+    T const& low() const { return low_; }
+    T const& high() const { return high_; }
+    T extents() const { return high_ - low_; }
+    typename T::Scalar extent(uint32_t i) const { return high_[i] - low_[i]; }
+    T center() const { return 0.5 * (low_ + high_); }
+
+    bool empty() const;
+    bool nonempty() const;
+    bool contains(T const& x) const;
+    bool contains(BoundingBox const& box) const;
+    bool intersects(BoundingBox const& box) const;
+    bool on_boundary(T const& x) const;
+
+    BoundingBox expansion(T const& x) const;
+    BoundingBox expansion(BoundingBox const& box) const;
+    BoundingBox intersection(BoundingBox const& box) const;
+    BoundingBox translation(T const& displacement) const;
+    BoundingBox offset(typename T::Scalar distance) const;
+    BoundingBox scaled(typename T::Scalar factor) const;
+
+    BoundingBox& expand(T const& x);
+    BoundingBox& expand(BoundingBox const& box);
+    BoundingBox& intersect(BoundingBox const& box);
+    BoundingBox& translate(T const& displacement);
+    BoundingBox& apply_offset(typename T::Scalar distance);
+    BoundingBox& scale(typename T::Scalar factor);
+
+private:
+    T low_;
+    T high_;
+};
+
+//..................................................................................................
+template <typename T>
+BoundingBox<T>& BoundingBox<T>::operator=(T const& x) {
+    low_ = x;
+    high_ = x;
+    return *this;
+}
+
+//..................................................................................................
+template <typename T>
+bool BoundingBox<T>::empty() const {
+    return (low_.array() > high_.array()).any();
+}
+
+//..................................................................................................
+template <typename T>
+bool BoundingBox<T>::nonempty() const {
+    return (low_.array() <= high_.array()).all();
+}
+
+//..................................................................................................
+template <typename T>
+bool BoundingBox<T>::contains(T const& x) const {
+    return (low_.array() <= x.array()).all() && (x.array() <= high_.array()).all();
+}
+
+//..................................................................................................
+template <typename T>
+bool BoundingBox<T>::contains(BoundingBox const& box) const {
+    return contains(box.low()) && contains(box.high());
+}
+
+//..................................................................................................
+template <typename T>
+bool BoundingBox<T>::intersects(BoundingBox const& box) const {
+    // An equivalent implementation would be
+    // return intersection(box).nonempty();
+    // But this version permits more expression template magic.
+    return (low_.cwiseMax(box.low()).array() <= high_.cwiseMin(box.high()).array()).all();
+}
+
+//..................................................................................................
+template <typename T>
+bool BoundingBox<T>::on_boundary(T const& x) const {
+    return (low_.array() == x.array()).any() || (x.array() == high_.array()).any();
+}
+
+//..................................................................................................
+template <typename T>
+BoundingBox<T> BoundingBox<T>::expansion(T const& x) const {
+    return BoundingBox(low_.cwiseMin(x), high_.cwiseMax(x));
+}
+
+//..................................................................................................
+template <typename T>
+BoundingBox<T> BoundingBox<T>::expansion(BoundingBox const& box) const {
+    return BoundingBox(low_.cwiseMin(box.low()), high_.cwiseMax(box.high()));
+}
+
+//..................................................................................................
+template <typename T>
+BoundingBox<T> BoundingBox<T>::intersection(BoundingBox const& box) const {
+    return BoundingBox(low_.cwiseMax(box.low()), high_.cwiseMin(box.high()));
+}
+
+//..................................................................................................
+template <typename T>
+BoundingBox<T> BoundingBox<T>::translation(T const& displacement) const {
+    return BoundingBox(low_ + displacement, high_ + displacement);
+}
+
+//..................................................................................................
+template <typename T>
+BoundingBox<T> BoundingBox<T>::offset(typename T::Scalar distance) const {
+    // Note: This only works for fixed size vectors. May need to update later.
+    return BoundingBox(low_ - distance * T::Ones(), high_ + distance * T::Ones());
+}
+
+//..................................................................................................
+template <typename T>
+BoundingBox<T> BoundingBox<T>::scaled(typename T::Scalar factor) const {
+    T const bb_center = center();
+    T const half_diagonal = factor * (low_ - bb_center);
+    return BoundingBox(bb_center + half_diagonal, bb_center - half_diagonal);
+}
+
+//..................................................................................................
+template <typename T>
+BoundingBox<T>& BoundingBox<T>::expand(T const& x) {
+    low_ = low_.cwiseMin(x);
+    high_ = high_.cwiseMax(x);
+    return *this;
+}
+
+//..................................................................................................
+template <typename T>
+BoundingBox<T>& BoundingBox<T>::expand(BoundingBox const& box) {
+    low_ = low_.cwiseMin(box.low());
+    high_ = high_.cwiseMax(box.high());
+    return *this;
+}
+
+//..................................................................................................
+template <typename T>
+BoundingBox<T>& BoundingBox<T>::intersect(BoundingBox const& box) {
+    low_ = low_.cwiseMax(box.low());
+    high_ = high_.cwiseMin(box.high());
+    return *this;
+}
+
+//..................................................................................................
+template <typename T>
+BoundingBox<T>& BoundingBox<T>::translate(T const& displacement) {
+    low_ += displacement;
+    high_ += displacement;
+    return *this;
+}
+
+//..................................................................................................
+template <typename T>
+BoundingBox<T>& BoundingBox<T>::apply_offset(typename T::Scalar distance) {
+    // Note: This only works for fixed size vectors. May need to update later.
+    low_ -= distance * T::Ones();
+    high_ += distance * T::Ones();
+    return *this;
+}
+
+//..................................................................................................
+template <typename T>
+BoundingBox<T>& BoundingBox<T>::scale(typename T::Scalar factor) {
+    T const bb_center = center();
+    T const half_diagonal = factor * (low_ - bb_center);
+    low_ = bb_center + half_diagonal;
+    high_ = bb_center - half_diagonal;
+    return *this;
+}
+
+//--------------------------------------------------------------------------------------------------
+template <typename T>
+std::ostream& operator<<(std::ostream& os, BoundingBox<T> const& bb) {
+    os << bb.low()[0] << ", " << bb.low()[1] << " to " << bb.high()[0] << ", " << bb.high()[1];
+    return os;
+}
+
+}
+
+#endif

simucene/glut_grapher.cpp

+#include "glut_grapher.h"
+#include <iostream>
+#include <iomanip>
+
+namespace simucene {
+
+//--------------------------------------------------------------------------------------------------
+template <typename T>
+void print_vector(VectorX<T> const& vec) {
+    std::cout << std::fixed << std::setprecision(3);
+    for (uint32_t i = 0; i < vec.size() - 1; ++i) {
+        std::cout << std::setw(5) << vec[i] << ", ";
+    }
+    std::cout << vec[vec.size() - 1] << '\n';
+}
+
+//..................................................................................................
+void GlutGrapher::initialize(
+    VectorX<Scalar>&& state,
+    SparseMatrix<Scalar>&& transport,
+    BoundingBox<Vector2<Scalar>> const& bounds
+) {
+    frame_ = 0;
+    state_ = std::move(state);
+    transport_ = std::move(transport);
+
+    auto const extents = bounds.extents();
+    d_ = extents[0] / state_.size();
+
+    aspect_ratio_ = extents[0] / extents[1];
+    Scalar window_aspect_ratio = static_cast<Scalar>(max_width_) / max_height_;
+    if (aspect_ratio_ >= window_aspect_ratio) {
+        // Width is the limiting factor.
+        window_size_[0] = max_width_;
+        window_size_[1] = max_width_ / aspect_ratio_;
+    } else {
+        // Height is the limiting factor.
+        window_size_[1] = max_height_;
+        window_size_[0] = max_height_ * aspect_ratio_;
+    }
+
+    // Set the correct zoom and translation to fit the simulation into the window.
+    translation_ = -bounds.low();
+    zoom_ = 1 / extents[0];
+}
+
+//..................................................................................................
+void GlutGrapher::init_glut() {
+    int argc = 0;
+    glutInit(&argc, nullptr);
+    glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH);
+    glutInitWindowSize(window_size_[0], window_size_[1]);
+    glutCreateWindow("glimage");
+    glShadeModel(GL_SMOOTH);
+    glMatrixMode(GL_PROJECTION);
+    glLoadIdentity();
+    glClearColor(1, 1, 1, 0);
+    gluOrtho2D(0, 1, 0, 1 / aspect_ratio_);
+}
+
+//..................................................................................................
+void GlutGrapher::display() const {
+    glClear(GL_COLOR_BUFFER_BIT);
+
+    // Draw the function
+    glColor3f(0, 0, 0);
+    glBegin(GL_LINE_STRIP);
+    for (int32_t i = 0; i < state_.size(); ++i) {
+        Scalar const x = (i + static_cast<Scalar>(0.5)) * d_;
+        Vector2<Scalar> const loc = zoom_ * (Vector2<Scalar>(x, state_[i]) + translation_);
+        glVertex2f(loc[0], loc[1]);
+    }
+    glEnd();
+
+    glFlush();
+}
+
+//..................................................................................................
+void GlutGrapher::idle() {
+    //std::cout << frame_ << ": ";
+    //print_vector(state_);
+    for (uint32_t i = 0; i < 100; ++i) {
+        state_ += transport_ * state_;
+        ++frame_;
+    }
+    glutSwapBuffers();
+    glutPostRedisplay();
+}
+
+}

simucene/glut_grapher.h

+#ifndef SIMUCENE_GLUT_GRAPHER_H
+#define SIMUCENE_GLUT_GRAPHER_H
+
+#include "scalar.h"
+#include "vector.h"
+#include "bounding_box.h"
+#include <GL/glut.h>
+
+namespace simucene {
+
+//--------------------------------------------------------------------------------------------------
+class GlutGrapher {
+public:
+    GlutGrapher() = default;
+    void initialize(
+        VectorX<Scalar>&& state,
+        SparseMatrix<Scalar>&& transport,
+        BoundingBox<Vector2<Scalar>> const& bounds
+    );
+
+    void init_glut();
+    void run() { glutMainLoop(); }
+    void display() const;
+    void idle();
+
+    Vector2<uint32_t> const& window_size() const { return window_size_; }
+    Scalar aspect_ratio() const { return aspect_ratio_; }
+
+public:
+    static constexpr uint32_t max_width_ = 1200;
+    static constexpr uint32_t max_height_ = 900;
+
+    Vector2<uint32_t> window_size_;
+    Scalar aspect_ratio_;
+    Vector2<Scalar> translation_;
+    Scalar zoom_;
+
+    uint32_t frame_;
+    VectorX<Scalar> state_;
+    SparseMatrix<Scalar> transport_;
+    Scalar d_;
+};
+
+}
+
+#endif

simucene/main.cpp

-#include "scalar.h"
 #include "matrix_builder.h"
+#include "glut_grapher.h"
 #include <iostream>
-#include <iomanip>
 
 using namespace simucene;
 
 //--------------------------------------------------------------------------------------------------
-template <typename T>
-void print_vector(VectorX<T> const& vec) {
-    std::cout << std::fixed << std::setprecision(3);
-    for (uint32_t i = 0; i < vec.size() - 1; ++i) {
-        std::cout << std::setw(5) << vec[i] << ", ";
-    }
-    std::cout << vec[vec.size() - 1] << '\n';
-}
+GlutGrapher grapher;
+void display() { grapher.display(); }
+void idle() { grapher.idle(); }
+void mouse(int, int, int, int) { exit(0); }
 
 //--------------------------------------------------------------------------------------------------
 SparseMatrix<Scalar> d_x_first_order(int32_t const dim) {
@@ -31,13 +26,22 @@ SparseMatrix<Scalar> diagonal_ones(int32_t const dim) {
 }
 
 //--------------------------------------------------------------------------------------------------
-int main(int const, char const**) {
+Scalar normal_pdf(Scalar x, Scalar mean, Scalar std_dev)
+{
+    // https://stackoverflow.com/questions/10847007/using-the-gaussian-probability-density-function-in-c
+    static Scalar const inv_sqrt_2pi = 0.3989422804014327;
+    Scalar a = (x - mean) / std_dev;
+    return inv_sqrt_2pi / std_dev * std::exp(static_cast<Scalar>(-0.5) * a * a);
+}
+
+//--------------------------------------------------------------------------------------------------
+int main(int, char**) {
     // primary parameters of the simulation
     Scalar const max_x = 1;
-    int32_t const dim = 20;
-    Scalar const h = 0.001;
-    Scalar const advection_coefficient = 0;
-    Scalar const diffusion_coefficient = 1;
+    int32_t const dim = 1000;
+    Scalar const h = 1e-8;
+    Scalar const advection_coefficient = -1e3;
+    Scalar const diffusion_coefficient = 1e1;
     VectorX<Scalar> reaction_coefficients = VectorX<Scalar>::Zero(dim);
 
     // derived quantities
@@ -49,11 +53,12 @@ int main(int const, char const**) {
     SparseMatrix<Scalar> const advection_matrix = alpha * d_x_first_order(dim);
     SparseMatrix<Scalar> const diffusion_matrix = delta * d_x2_first_order(dim);
     SparseMatrix<Scalar> const reaction_matrix = h * diagonal_ones(dim) * reaction_coefficients.asDiagonal();
-    SparseMatrix<Scalar> const transport_matrix = advection_matrix + diffusion_matrix + reaction_matrix;
+    SparseMatrix<Scalar> transport_matrix = advection_matrix + diffusion_matrix + reaction_matrix;
 
     // debug: print useful things
     std::cout << "alpha: " << alpha << '\n';
     std::cout << "delta: " << delta << '\n';
+    /*
     std::cout << "Advection\n";
     std::cout << advection_matrix << '\n';
     std::cout << "Diffusion\n";
@@ -62,14 +67,17 @@ int main(int const, char const**) {
     std::cout << reaction_matrix << '\n';
     std::cout << "Transport\n";
     std::cout << transport_matrix << '\n';
+    */
 
-    // define the initial state (not boundary conditions)
+    // define the initial state 
     VectorX<Scalar> state = VectorX<Scalar>::Zero(dim);
-    for (int32_t i = dim / 2 - 2; i < dim / 2 + 2; ++i) {
-        state[i] = 1;
+    for (int32_t i = 0; i < dim; ++i) {
+        Scalar const x = (i + static_cast<Scalar>(0.5)) * d;
+        state[i] = 0.1 * normal_pdf(x, 0.2 * max_x, 0.05 * max_x);
     }
 
     // iterate
+    /*
     for (uint32_t i = 0; i < 100; ++i) {
         std::cout << i << ": ";
         print_vector(state);
@@ -77,6 +85,18 @@ int main(int const, char const**) {
     }
     std::cout << 100 << ": ";
     print_vector(state);
+    */
+
+    grapher.initialize(
+        std::move(state),
+        std::move(transport_matrix),
+        BoundingBox<Vector2<Scalar>>({0, -0.1}, {max_x, 1 + 0.1})
+    );
+    grapher.init_glut();
+    glutDisplayFunc(display);
+    glutMouseFunc(mouse);
+    glutIdleFunc(idle);
+    grapher.run();
 
     return 0;
 }