Merge pull request #2382 from RudolfWeeber/rotation_cleanup

Rotation cleanup

src/core/communication.cpp

@@ -644,16 +644,17 @@ void mpi_send_rotational_inertia_slave(int pnode, int part) {
 #endif
 }
 
-void mpi_rotate_particle(int pnode, int part, double axis[3], double angle) {
+void mpi_rotate_particle(int pnode, int part, const Vector3d &axis,
+ double angle) {
 #ifdef ROTATION
  mpi_call(mpi_rotate_particle_slave, pnode, part);
 
  if (pnode == this_node) {
  Particle *p = local_particles[part];
- local_rotate_particle(p, axis, angle);
+ local_rotate_particle(*p, axis, angle);
  } else {
- MPI_Send(axis, 3, MPI_DOUBLE, pnode, SOME_TAG, comm_cart);
- MPI_Send(&angle, 1, MPI_DOUBLE, pnode, SOME_TAG, comm_cart);
+ comm_cart.send(pnode, SOME_TAG, axis);
+ comm_cart.send(pnode, SOME_TAG, angle);
  }
 
  on_particle_change();
@@ -664,10 +665,11 @@ void mpi_rotate_particle_slave(int pnode, int part) {
 #ifdef ROTATION
  if (pnode == this_node) {
  Particle *p = local_particles[part];
- double axis[3], angle;
- MPI_Recv(axis, 3, MPI_DOUBLE, 0, SOME_TAG, comm_cart, MPI_STATUS_IGNORE);
- MPI_Recv(&angle, 1, MPI_DOUBLE, 0, SOME_TAG, comm_cart, MPI_STATUS_IGNORE);
- local_rotate_particle(p, axis, angle);
+ Vector3d axis;
+ double angle;
+ comm_cart.recv(0, SOME_TAG, axis);
+ comm_cart.recv(0, SOME_TAG, angle);
+ local_rotate_particle(*p, axis, angle);
  }
 
  on_particle_change();
@@ -815,18 +817,15 @@ void mpi_send_quat_slave(int pnode, int part) {
 
 /********************* REQ_SET_OMEGA ********/
 
-void mpi_send_omega(int pnode, int part, double omega[3]) {
+void mpi_send_omega(int pnode, int part, const Vector3d &omega) {
 #ifdef ROTATION
  mpi_call(mpi_send_omega_slave, pnode, part);
 
  if (pnode == this_node) {
  Particle *p = local_particles[part];
- /* memmove(p->omega, omega, 3*sizeof(double));*/
- p->m.omega[0] = omega[0];
- p->m.omega[1] = omega[1];
- p->m.omega[2] = omega[2];
+ p->m.omega = omega;
  } else {
- MPI_Send(omega, 3, MPI_DOUBLE, pnode, SOME_TAG, comm_cart);
+ comm_cart.send(pnode, SOME_TAG, omega);
  }
 
  on_particle_change();
@@ -837,27 +836,23 @@ void mpi_send_omega_slave(int pnode, int part) {
 #ifdef ROTATION
  if (pnode == this_node) {
  Particle *p = local_particles[part];
- MPI_Recv(p->m.omega.data(), 3, MPI_DOUBLE, 0, SOME_TAG, comm_cart,
- MPI_STATUS_IGNORE);
+ comm_cart.recv(0, SOME_TAG, p->m.omega);
  }
-
  on_particle_change();
 #endif
 }
 
 /********************* REQ_SET_TORQUE ********/
 
-void mpi_send_torque(int pnode, int part, double torque[3]) {
+void mpi_send_torque(int pnode, int part, const Vector3d &torque) {
 #ifdef ROTATION
  mpi_call(mpi_send_torque_slave, pnode, part);
 
  if (pnode == this_node) {
  Particle *p = local_particles[part];
- p->f.torque[0] = torque[0];
- p->f.torque[1] = torque[1];
- p->f.torque[2] = torque[2];
+ p->f.torque = torque;
  } else {
- MPI_Send(torque, 3, MPI_DOUBLE, pnode, SOME_TAG, comm_cart);
+ comm_cart.send(pnode, SOME_TAG, torque);
  }
 
  on_particle_change();
@@ -868,8 +863,7 @@ void mpi_send_torque_slave(int pnode, int part) {
 #ifdef ROTATION
  if (pnode == this_node) {
  Particle *p = local_particles[part];
- MPI_Recv(p->f.torque.data(), 3, MPI_DOUBLE, 0, SOME_TAG, comm_cart,
- MPI_STATUS_IGNORE);
+ comm_cart.recv(0, SOME_TAG, p->f.torque);
  }
 
  on_particle_change();

src/core/communication.hpp

@@ -232,7 +232,8 @@ void mpi_send_rotational_inertia(int node, int part, double rinertia[3]);
  \param axis rotation axis
  \param angle rotation angle
 */
-void mpi_rotate_particle(int node, int part, double axis[3], double angle);
+void mpi_rotate_particle(int pnode, int part, const Vector3d &axis,
+ double angle);
 #endif
 
 #ifdef AFFINITY
@@ -278,15 +279,15 @@ void mpi_send_rotation(int pnode, int part, short int rot);
  \param node the node it is attached to.
  \param omega its new angular velocity.
 */
-void mpi_send_omega(int node, int part, double omega[3]);
+void mpi_send_omega(int node, int part, const Vector3d &omega);
 
 /** Issue REQ_SET_TORQUE: send particle torque.
  Also calls \ref on_particle_change.
  \param part the particle.
  \param node the node it is attached to.
  \param torque its new torque.
 */
-void mpi_send_torque(int node, int part, double torque[3]);
+void mpi_send_torque(int node, int part, const Vector3d &torque);
 #endif
 
 #ifdef DIPOLES

src/core/minimize_energy.cpp

@@ -95,18 +95,11 @@ bool steepest_descent_step(void) {
  }
 #ifdef ROTATION
  // Rotational increment
- double dq[3]; // Vector parallel to torque
+ Vector3d dq = params->gamma * p.f.torque;
+ t += p.f.torque.norm2();
 
- for (int j = 0; j < 3; j++) {
- dq[j] = 0;
- // Square of torque
- t += Utils::sqr(p.f.torque[j]);
-
- // Rotational increment
- dq[j] = params->gamma * p.f.torque[j];
- }
  // Normalize rotation axis and compute amount of rotation
- double l = normr(dq);
+ double l = dq.norm();
  if (l > 0.0) {
  for (int j = 0; j < 3; j++)
  dq[j] /= l;
@@ -116,7 +109,7 @@ bool steepest_descent_step(void) {
  l = sgn(l) * params->max_displacement;
 
  // Rotate the particle around axis dq by amount l
- local_rotate_particle(&(p), dq, l);
+ local_rotate_particle(p, dq, l);
  }
 #endif
 

src/core/observables/ParticleAngularVelocities.hpp

@@ -33,20 +33,8 @@ class ParticleAngularVelocities : public PidObservable {
  std::vector<double> res(n_values());
  for (int i = 0; i < ids().size(); i++) {
 #ifdef ROTATION
-
- double RMat[9];
- double omega[3];
- define_rotation_matrix(partCfg[ids()[i]], RMat);
- omega[0] = RMat[0 + 3 * 0] * partCfg[ids()[i]].m.omega[0] +
- RMat[1 + 3 * 0] * partCfg[ids()[i]].m.omega[1] +
- RMat[2 + 3 * 0] * partCfg[ids()[i]].m.omega[2];
- omega[1] = RMat[0 + 3 * 1] * partCfg[ids()[i]].m.omega[0] +
- RMat[1 + 3 * 1] * partCfg[ids()[i]].m.omega[1] +
- RMat[2 + 3 * 1] * partCfg[ids()[i]].m.omega[2];
- omega[2] = RMat[0 + 3 * 2] * partCfg[ids()[i]].m.omega[0] +
- RMat[1 + 3 * 2] * partCfg[ids()[i]].m.omega[1] +
- RMat[2 + 3 * 2] * partCfg[ids()[i]].m.omega[2];
-
+ const Vector3d omega =
+ convert_vector_body_to_space(partCfg[i], partCfg[i].m.omega);
  res[3 * i + 0] = omega[0];
  res[3 * i + 1] = omega[1];
  res[3 * i + 2] = omega[2];

src/core/observables/ParticleBodyVelocities.hpp

@@ -35,20 +35,8 @@ class ParticleBodyVelocities : public PidObservable {
 #ifdef ROTATION
 
  double RMat[9];
- double vel_lab[3];
- double vel_body[3];
-
- vel_lab[0] = partCfg[ids()[i]].m.v[0];
- vel_lab[1] = partCfg[ids()[i]].m.v[1];
- vel_lab[2] = partCfg[ids()[i]].m.v[2];
- define_rotation_matrix(partCfg[ids()[i]], RMat);
-
- vel_body[0] = RMat[0 + 3 * 0] * vel_lab[0] +
- RMat[0 + 3 * 1] * vel_lab[1] + RMat[0 + 3 * 2] * vel_lab[2];
- vel_body[1] = RMat[1 + 3 * 0] * vel_lab[0] +
- RMat[1 + 3 * 1] * vel_lab[1] + RMat[1 + 3 * 2] * vel_lab[2];
- vel_body[2] = RMat[2 + 3 * 0] * vel_lab[0] +
- RMat[2 + 3 * 1] * vel_lab[1] + RMat[2 + 3 * 2] * vel_lab[2];
+ const Vector3d vel_body =
+ convert_vector_space_to_body(partCfg[i], partCfg[i].m.v);
 
  res[3 * i + 0] = vel_body[0];
  res[3 * i + 1] = vel_body[1];

src/core/particle_data.cpp

@@ -530,7 +530,7 @@ int set_particle_rotation(int part, int rot) {
 }
 #endif
 #ifdef ROTATION
-int rotate_particle(int part, double axis[3], double angle) {
+int rotate_particle(int part, const Vector3d &axis, double angle) {
  auto const pnode = get_particle_node(part);
 
  mpi_rotate_particle(pnode, part, axis, angle);
@@ -664,59 +664,31 @@ int set_particle_quat(int part, double quat[4]) {
  return ES_OK;
 }
 
-int set_particle_omega_lab(int part, double omega_lab[3]) {
+int set_particle_omega_lab(int part, const Vector3d &omega_lab) {
  auto const &particle = get_particle_data(part);
 
- /* Internal functions require the body coordinates
- so we need to convert to these from the lab frame */
-
- double A[9];
- double omega[3];
-
- define_rotation_matrix(particle, A);
-
- omega[0] = A[0 + 3 * 0] * omega_lab[0] + A[0 + 3 * 1] * omega_lab[1] +
- A[0 + 3 * 2] * omega_lab[2];
- omega[1] = A[1 + 3 * 0] * omega_lab[0] + A[1 + 3 * 1] * omega_lab[1] +
- A[1 + 3 * 2] * omega_lab[2];
- omega[2] = A[2 + 3 * 0] * omega_lab[0] + A[2 + 3 * 1] * omega_lab[1] +
- A[2 + 3 * 2] * omega_lab[2];
-
  auto const pnode = get_particle_node(part);
- mpi_send_omega(pnode, part, omega);
+ mpi_send_omega(pnode, part,
+ convert_vector_space_to_body(particle, omega_lab));
  return ES_OK;
 }
 
-int set_particle_omega_body(int part, double omega[3]) {
+int set_particle_omega_body(int part, const Vector3d &omega) {
  auto const pnode = get_particle_node(part);
  mpi_send_omega(pnode, part, omega);
  return ES_OK;
 }
 
-int set_particle_torque_lab(int part, double torque_lab[3]) {
+int set_particle_torque_lab(int part, const Vector3d &torque_lab) {
  auto const &particle = get_particle_data(part);
 
- /* Internal functions require the body coordinates
- so we need to convert to these from the lab frame */
-
- double A[9];
- double torque[3];
-
- define_rotation_matrix(particle, A);
-
- torque[0] = A[0 + 3 * 0] * torque_lab[0] + A[0 + 3 * 1] * torque_lab[1] +
- A[0 + 3 * 2] * torque_lab[2];
- torque[1] = A[1 + 3 * 0] * torque_lab[0] + A[1 + 3 * 1] * torque_lab[1] +
- A[1 + 3 * 2] * torque_lab[2];
- torque[2] = A[2 + 3 * 0] * torque_lab[0] + A[2 + 3 * 1] * torque_lab[1] +
- A[2 + 3 * 2] * torque_lab[2];
-
  auto const pnode = get_particle_node(part);
- mpi_send_torque(pnode, part, torque);
+ mpi_send_torque(pnode, part,
+ convert_vector_space_to_body(particle, torque_lab));
  return ES_OK;
 }
 
-int set_particle_torque_body(int part, double torque[3]) {
+int set_particle_torque_body(int part, const Vector3d &torque) {
  auto const pnode = get_particle_node(part);
 
  /* Nothing to be done but pass, since the coordinates
@@ -1243,10 +1215,6 @@ void pointer_to_omega_body(Particle const *p, double const *&res) {
  res = p->m.omega.data();
 }
 
-void pointer_to_torque_lab(Particle const *p, double const *&res) {
- res = p->f.torque.data();
-}
-
 void pointer_to_quat(Particle const *p, double const *&res) {
  res = p->r.quat.data();
 }

src/core/particle_data.hpp

@@ -672,7 +672,7 @@ int set_particle_rotation(int part, int rot);
  @param axis rotation axis
  @param angle rotation angle
 */
-int rotate_particle(int part, double axis[3], double angle);
+int rotate_particle(int part, const Vector3d &axis, double angle);
 
 #ifdef AFFINITY
 /** Call only on the master node: set particle affinity.
@@ -736,28 +736,28 @@ int set_particle_quat(int part, double quat[4]);
  @param omega its new angular velocity.
  @return ES_OK if particle existed
 */
-int set_particle_omega_lab(int part, double omega[3]);
+int set_particle_omega_lab(int part, const Vector3d &omega);
 
 /** Call only on the master node: set particle angular velocity in body frame.
  @param part the particle.
  @param omega its new angular velocity.
  @return ES_OK if particle existed
 */
-int set_particle_omega_body(int part, double omega[3]);
+int set_particle_omega_body(int part, const Vector3d &omega);
 
 /** Call only on the master node: set particle torque from lab frame.
  @param part the particle.
  @param torque its new torque.
  @return ES_OK if particle existed
 */
-int set_particle_torque_lab(int part, double torque[3]);
+int set_particle_torque_lab(int part, const Vector3d &torque);
 
 /** Call only on the master node: set particle torque in body frame.
  @param part the particle.
  @param torque its new torque.
  @return ES_OK if particle existed
 */
-int set_particle_torque_body(int part, double torque[3]);
+int set_particle_torque_body(int part, const Vector3d &torque);
 #endif
 
 #ifdef DIPOLES
@@ -995,7 +995,7 @@ int number_of_particles_with_type(int type);
 #ifdef ROTATION
 void pointer_to_omega_body(Particle const *p, double const *&res);
 
-void pointer_to_torque_lab(Particle const *p, double const *&res);
+inline Vector3d get_torque_body(const Particle &p) { return p.f.torque; }
 
 void pointer_to_quat(Particle const *p, double const *&res);
 

src/core/rotate_system.cpp

@@ -48,7 +48,7 @@ void local_rotate_system(double phi, double theta, double alpha) {
  mpi::all_reduce(comm_cart, local_com, std::plus<Vector3d>()) / total_mass;
 
  // Rotation axis in Cartesian coordinates
- double axis[3];
+ Vector3d axis;
  axis[0] = sin(theta) * cos(phi);
  axis[1] = sin(theta) * sin(phi);
  axis[2] = cos(theta);
@@ -67,7 +67,7 @@ void local_rotate_system(double phi, double theta, double alpha) {
  p.r.p[j] = com[j] + res[j];
  }
 #ifdef ROTATION
- local_rotate_particle(&p, axis, alpha);
+ local_rotate_particle(p, axis, alpha);
 #endif
  }