fix: save exact particles state

apps/benchmarks/CNT/Config.cpp

@@ -50,10 +50,11 @@ void saveConfig(const shared_ptr<data::ParticleStorage>& ps,
       sb << static_cast<double> (p.getPosition()[1]);
       sb << static_cast<double> (p.getPosition()[2]);
 
-      Vec3 orient = p.getRotation().getMatrix() * Vec3(1,0,0);
-      // Now orient has the shape (sin(th)cos(ph), sin(th)sin(ph), cos(th))
-      sb << static_cast<double> (std::acos(orient[2]));
-      sb << static_cast<double> (std::atan2(orient[1], orient[0]));
+      Quaternion q = p.getRotation().getQuaternion();
+      sb << static_cast<double> (q[0]);
+      sb << static_cast<double> (q[1]);
+      sb << static_cast<double> (q[2]);
+      sb << static_cast<double> (q[3]);
 
       sb << static_cast<double> (p.getLinearVelocity()[0]);
       sb << static_cast<double> (p.getLinearVelocity()[1]);
@@ -62,6 +63,15 @@ void saveConfig(const shared_ptr<data::ParticleStorage>& ps,
       sb << static_cast<double> (p.getAngularVelocity()[0]);
       sb << static_cast<double> (p.getAngularVelocity()[1]);
       sb << static_cast<double> (p.getAngularVelocity()[2]);
+
+      sb << static_cast<double> (p.getOldForce()[0]);
+      sb << static_cast<double> (p.getOldForce()[1]);
+      sb << static_cast<double> (p.getOldForce()[2]);
+
+      sb << static_cast<double> (p.getOldTorque()[0]);
+      sb << static_cast<double> (p.getOldTorque()[1]);
+      sb << static_cast<double> (p.getOldTorque()[2]);
+
    }
 
    switch (io_mode)
@@ -74,7 +84,7 @@ void saveConfig(const shared_ptr<data::ParticleStorage>& ps,
          WALBERLA_ROOT_SECTION()
          {
             uint_t dataSize = sizeof(mpi::RecvBuffer::ElementType) * rb.size();
-            size_t size     = rb.size() / 104;
+            size_t size     = rb.size() / 168;
             std::ofstream binfile;
             binfile.open(filename + ".sav", std::ios::out | std::ios::binary);
             binfile.write((char*) &size, sizeof(size_t));

apps/benchmarks/CNT/InitializeCNTs.cpp

@@ -164,27 +164,43 @@ int64_t loadCNTs(const std::string& filename,
          double x;
          double y;
          double z;
-         double theta;
-         double phi;
+         double q0;
+         double q1;
+         double q2;
+         double q3;
          double vx;
          double vy;
          double vz;
          double wx;
          double wy;
          double wz;
+         double fx;
+         double fy;
+         double fz;
+         double tx;
+         double ty;
+         double tz;
          binfile.read((char *) &sID, sizeof(int64_t));
          binfile.read((char *) &cID, sizeof(int64_t));
          binfile.read((char *) &x, sizeof(double));
          binfile.read((char *) &y, sizeof(double));
          binfile.read((char *) &z, sizeof(double));
-         binfile.read((char *) &theta, sizeof(double));
-         binfile.read((char *) &phi, sizeof(double));
+         binfile.read((char *) &q0, sizeof(double));
+         binfile.read((char *) &q1, sizeof(double));
+         binfile.read((char *) &q2, sizeof(double));
+         binfile.read((char *) &q3, sizeof(double));
          binfile.read((char *) &vx, sizeof(double));
          binfile.read((char *) &vy, sizeof(double));
          binfile.read((char *) &vz, sizeof(double));
          binfile.read((char *) &wx, sizeof(double));
          binfile.read((char *) &wy, sizeof(double));
          binfile.read((char *) &wz, sizeof(double));
+         binfile.read((char *) &fx, sizeof(double));
+         binfile.read((char *) &fy, sizeof(double));
+         binfile.read((char *) &fz, sizeof(double));
+         binfile.read((char *) &tx, sizeof(double));
+         binfile.read((char *) &ty, sizeof(double));
+         binfile.read((char *) &tz, sizeof(double));
          Vec3 pos;
          pos[0] = real_c(x);
          pos[1] = real_c(y);
@@ -198,10 +214,11 @@ int64_t loadCNTs(const std::string& filename,
             sp.setInteractionRadius(kernel::cnt::outer_radius);
             sp.setSegmentID(sID);
             sp.setClusterID(cID);
-            sp.getRotationRef().rotate(Vec3(0_r, 1_r, 0_r), -0.5_r * math::pi + real_c(theta));
-            sp.getRotationRef().rotate(Vec3(0_r, 0_r, 1_r), real_c(phi));
+            sp.setRotation(Rot3(Quaternion(q0, q1, q2, q3)));
             sp.setLinearVelocity(Vec3(real_c(vx), real_c(vy), real_c(vz)));
             sp.setAngularVelocity(Vec3(real_c(wx), real_c(wy), real_c(wz)));
+            sp.setOldForce(Vec3(real_c(fx), real_c(fy), real_c(fz)));
+            sp.setOldTorque(Vec3(real_c(tx), real_c(ty), real_c(tz)));
             numParticles++;
          }
       }