Updated main.cpp

CMakeLists.txt

@@ -25,7 +25,7 @@ include_directories(SYSTEM deps/libtimestep/include)
 include_directories(SYSTEM deps/eigen)
 include_directories(include)
 
-add_executable(libgran main.cpp writer.cpp)
+add_executable(libgran main.cpp writer.cpp test/compute_energy.cpp)
 # A particle and a lattice colliding
 add_executable(contact_test test/contact.cpp test/compute_energy.cpp)
 # A particle and a lattice with VdW interactions colliding

main.cpp

@@ -12,7 +12,8 @@
 #include <libgran/hamaker_force/hamaker_force.h>
 #include <libgran/granular_system/granular_system.h>
 
-#include "../writer.h"
+#include "writer.h"
+#include "test/compute_energy.h"
 
 // The driver program is responsible for:
 // 1) Instantiating the force models
@@ -108,11 +109,23 @@ int main() {
                                                              v0, theta0, omega0, 0.0, Eigen::Vector3d::Zero(), 0.0,
                                                              step_handler_instance, contact_force_model, hamaker_model);
 
+
+    // Output stream for data
+    std::ofstream ofs("energies.dat");
+
+    ofs << "t\tE_trs\tE_rot\tE_tot\tP\tL" << std::endl;
+
     auto start_time = std::chrono::high_resolution_clock::now();
     for (size_t n = 0; n < n_steps; n ++) {
         if (n % dump_period == 0) {
             std::cout << "Dump #" << n / dump_period << std::endl;
             write_particles("run", system.get_x(), system.get_theta(), r_part);
+            std::cout << double(n) * dt << "\t"
+                << compute_translational_kinetic_energy(system.get_v(), mass) << "\t"
+                << compute_rotational_kinetic_energy(system.get_v(), mass) << "\t"
+                << compute_total_kinetic_energy(system.get_v(), system.get_omega(), mass, inertia) << "\t"
+                << compute_linear_momentum(system.get_v(), mass) << "\t"
+                << compute_angular_momentum(system.get_x(), system.get_v(), mass, system.get_omega(), inertia) << std::endl;
         }
         system.do_step(dt);
     }

test/compute_energy.cpp

@@ -4,9 +4,9 @@
 
 #include "compute_energy.h"
 
-double compute_kinetic_energy(std::vector<Eigen::Vector3d> const & v,
-                              std::vector<Eigen::Vector3d> const & omega,
-                              double mass, double inertia) {
+double compute_total_kinetic_energy(std::vector<Eigen::Vector3d> const & v,
+                                    std::vector<Eigen::Vector3d> const & omega,
+                                    double mass, double inertia) {
 
     double energy = 0.0;
 
@@ -27,3 +27,32 @@ double compute_linear_momentum(std::vector<Eigen::Vector3d> const & v, double ma
 
     return mass * momentum.norm();
 }
+
+double compute_rotational_kinetic_energy(std::vector<Eigen::Vector3d> const & omegas, double inertia) {
+    double val = 0.0;
+    for (auto const & omega : omegas) {
+        val += inertia * omega.dot(omega);
+    }
+    return val / 2.0;
+}
+
+double compute_translational_kinetic_energy(std::vector<Eigen::Vector3d> const & vs, double m) {
+    double val = 0.0;
+    for (auto const & v : vs) {
+        val += m * v.dot(v);
+    }
+    return val / 2.0;
+}
+
+double compute_angular_momentum(std::vector<Eigen::Vector3d> const & xs, std::vector<Eigen::Vector3d> const & vs, double m, std::vector<Eigen::Vector3d> const & omegas, double inertia) {
+    Eigen::Vector3d val = Eigen::Vector3d::Zero();
+    for (auto const & omega : omegas) {
+        val += inertia * omega;
+    }
+    for (size_t i = 0; i < xs.size(); i ++) {
+        auto const & x = xs[i];
+        auto const & v = vs[i];
+        val += x.cross(m*v);
+    }
+    return val.norm();
+}

test/compute_energy.h

@@ -9,10 +9,18 @@
 
 #include <Eigen/Eigen>
 
-double compute_kinetic_energy(std::vector<Eigen::Vector3d> const & v,
-                              std::vector<Eigen::Vector3d> const & omega,
-                              double mass, double inertia);
+double compute_total_kinetic_energy(std::vector<Eigen::Vector3d> const & v,
+                                    std::vector<Eigen::Vector3d> const & omega,
+                                    double mass, double inertia);
+
+double compute_linear_momentum(std::vector<Eigen::Vector3d> const & vs, double mass);
+
+double compute_rotational_kinetic_energy(std::vector<Eigen::Vector3d> const & omegas, double inertia);
+
+double compute_translational_kinetic_energy(std::vector<Eigen::Vector3d> const & vs, double m);
+
+double compute_angular_momentum(std::vector<Eigen::Vector3d> const & xs, std::vector<Eigen::Vector3d> const & vs, double m,
+                                std::vector<Eigen::Vector3d> const & omegas, double inertia);
 
-double compute_linear_momentum(std::vector<Eigen::Vector3d> const & v, double mass);
 
 #endif //LIBGRAN_COMPUTE_ENERGY_H

test/contact.cpp

@@ -95,7 +95,7 @@ int main() {
     double ke_tolerance = 1.0; // Percent
 
     double linear_momentum = compute_linear_momentum(system.get_v(), mass);
-    double ke = compute_kinetic_energy(system.get_v(), system.get_omega(), mass, inertia);
+    double ke = compute_total_kinetic_energy(system.get_v(), system.get_omega(), mass, inertia);
 
     if (abs(ke - target_ke) / target_ke > ke_tolerance / 100.0 ||
         abs(linear_momentum - target_linear_momentum) / target_linear_momentum > linear_momentum_tolerance / 100.0)

test/hamaker.cpp

@@ -109,7 +109,7 @@ int main() {
     double rg_tolerance = 1.0; // Percent
 
     double linear_momentum = compute_linear_momentum(system.get_v(), mass);
-    double ke = compute_kinetic_energy(system.get_v(), system.get_omega(), mass, inertia);
+    double ke = compute_total_kinetic_energy(system.get_v(), system.get_omega(), mass, inertia);
     double rg = radius_of_gyration(system.get_x());
 
     if (abs(ke - target_ke) / target_ke > ke_tolerance / 100.0 ||