Quadrupoles are now constructed

src/const.h

@@ -40,7 +40,7 @@
 #define const_max_u_change 0.1f
 
 /* Gravity stuff. */
-#define multipole_order 1
+#define multipole_order 2
 #define const_theta_max 0.57735f
 #define const_G 6.672e-8f     /* Gravitational constant. */
 #define const_epsilon 0.0014f /* Gravity blending distance. */

src/multipole.c

@@ -22,44 +22,140 @@
 #include "../config.h"
 
 /* Some standard headers. */
-#include <float.h>
-#include <limits.h>
-#include <math.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-
-/* MPI headers. */
-#ifdef WITH_MPI
-#include <mpi.h>
-#endif
+#include <strings.h>
 
 /* This object's header. */
 #include "multipole.h"
 
+/**
+* @brief Reset the data of a #multipole.
+*
+* @param m The #multipole.
+*/
+void multipole_reset(struct multipole *m) {
+
+  /* Just bzero the struct. */
+  bzero(m, sizeof(struct multipole));
+}
+
+/**
+* @brief Init a multipole from a set of particles.
+*
+* @param m The #multipole.
+* @param gparts The #gpart.
+* @param gcount The number of particles.
+*/
+void multipole_init(struct multipole *m, const struct gpart *gparts,
+                    int gcount) {
+
+#if multipole_order > 2
+#error "Multipoles of order >2 not yet implemented."
+#endif
+
+  /* Zero everything */
+  multipole_reset(m);
+
+  /* Temporary variables */
+  double mass = 0.0;
+  double com[3] = {0.0, 0.0, 0.0};
+
+#if multipole_order >= 2
+  double I_xx = 0.0, I_yy = 0.0, I_zz = 0.0;
+  double I_xy = 0.0, I_xz = 0.0, I_yz = 0.0;
+#endif
+
+  /* Collect the particle data. */
+  for (int k = 0; k < gcount; k++) {
+    const float w = gparts[k].mass;
+
+    mass += w;
+    com[0] += gparts[k].x[0] * w;
+    com[1] += gparts[k].x[1] * w;
+    com[2] += gparts[k].x[2] * w;
+
+#if multipole_order >= 2
+    I_xx += gparts[k].x[0] * gparts[k].x[0] * w;
+    I_yy += gparts[k].x[1] * gparts[k].x[1] * w;
+    I_zz += gparts[k].x[2] * gparts[k].x[2] * w;
+    I_xy += gparts[k].x[0] * gparts[k].x[1] * w;
+    I_xz += gparts[k].x[0] * gparts[k].x[2] * w;
+    I_yz += gparts[k].x[1] * gparts[k].x[2] * w;
+#endif
+  }
+
+  const double imass = 1.0 / mass;
+
+  /* Store the data on the multipole. */
+  m->mass = mass;
+  m->CoM[0] = com[0] * imass;
+  m->CoM[1] = com[1] * imass;
+  m->CoM[2] = com[2] * imass;
+
+#if multipole_order >= 2
+  m->I_xx = I_xx - imass * com[0] * com[0];
+  m->I_yy = I_yy - imass * com[1] * com[1];
+  m->I_zz = I_zz - imass * com[2] * com[2];
+  m->I_xy = I_xy - imass * com[0] * com[1];
+  m->I_xz = I_xz - imass * com[0] * com[2];
+  m->I_yz = I_yz - imass * com[1] * com[2];
+#endif
+}
+
 /**
  * @brief Add the second multipole to the first one.
  *
- * @param m_sum The #multipole which will contain the sum.
- * @param m_other The #multipole to add.
+ * @param ma The #multipole which will contain the sum.
+ * @param mb The #multipole to add.
  */
 
-void multipole_add(struct multipole *m_sum, const struct multipole *m_other) {
+void multipole_add(struct multipole *ma, const struct multipole *mb) {
 
-#if multipole_order != 1
-#error "Multipoles of order >1 not yet implemented."
+#if multipole_order > 2
+#error "Multipoles of order >2 not yet implemented."
 #endif
 
   /* Correct the position. */
-  const float M_tot = m_sum->mass + m_other->mass;
-  const float M_tot_inv = 1.f / M_tot;
-  for (int k = 0; k < 3; k++)
-    m_sum->CoM[k] =
-        (m_sum->CoM[k] * m_sum->mass + m_other->CoM[k] * m_other->mass) *
-        M_tot_inv;
-
-  /* Add the particle to the moments. */
-  m_sum->mass = M_tot;
+  const double ma_mass = ma->mass;
+  const double mb_mass = mb->mass;
+  const double M_tot = ma_mass + mb_mass;
+  const double M_tot_inv = 1.0 / M_tot;
+
+  const double ma_CoM[3] = {ma->CoM[0], ma->CoM[1], ma->CoM[2]};
+  const double mb_CoM[3] = {mb->CoM[0], mb->CoM[1], mb->CoM[2]};
+
+#if multipole_order >= 2
+  const double ma_I_xx = (double)ma->I_xx + ma_mass * ma_CoM[0] * ma_CoM[0];
+  const double ma_I_yy = (double)ma->I_yy + ma_mass * ma_CoM[1] * ma_CoM[1];
+  const double ma_I_zz = (double)ma->I_zz + ma_mass * ma_CoM[2] * ma_CoM[2];
+  const double ma_I_xy = (double)ma->I_xy + ma_mass * ma_CoM[0] * ma_CoM[1];
+  const double ma_I_xz = (double)ma->I_xz + ma_mass * ma_CoM[0] * ma_CoM[2];
+  const double ma_I_yz = (double)ma->I_yz + ma_mass * ma_CoM[1] * ma_CoM[2];
+
+  const double mb_I_xx = (double)mb->I_xx + mb_mass * mb_CoM[0] * mb_CoM[0];
+  const double mb_I_yy = (double)mb->I_yy + mb_mass * mb_CoM[1] * mb_CoM[1];
+  const double mb_I_zz = (double)mb->I_zz + mb_mass * mb_CoM[2] * mb_CoM[2];
+  const double mb_I_xy = (double)mb->I_xy + mb_mass * mb_CoM[0] * mb_CoM[1];
+  const double mb_I_xz = (double)mb->I_xz + mb_mass * mb_CoM[0] * mb_CoM[2];
+  const double mb_I_yz = (double)mb->I_yz + mb_mass * mb_CoM[1] * mb_CoM[2];
+#endif
+
+  /* New mass */
+  ma->mass = M_tot;
+
+  /* New CoM */
+  ma->CoM[0] = (ma_CoM[0] * ma_mass + mb_CoM[0] * mb_mass) * M_tot_inv;
+  ma->CoM[1] = (ma_CoM[1] * ma_mass + mb_CoM[1] * mb_mass) * M_tot_inv;
+  ma->CoM[2] = (ma_CoM[2] * ma_mass + mb_CoM[2] * mb_mass) * M_tot_inv;
+
+/* New quadrupole */
+#if multipole_order >= 2
+  ma->I_xx = (ma_I_xx + mb_I_xx) - M_tot * ma->CoM[0] * ma->CoM[0];
+  ma->I_yy = (ma_I_yy + mb_I_yy) - M_tot * ma->CoM[1] * ma->CoM[1];
+  ma->I_zz = (ma_I_zz + mb_I_zz) - M_tot * ma->CoM[2] * ma->CoM[2];
+  ma->I_xy = (ma_I_xy + mb_I_xy) - M_tot * ma->CoM[0] * ma->CoM[1];
+  ma->I_xz = (ma_I_xz + mb_I_xz) - M_tot * ma->CoM[0] * ma->CoM[2];
+  ma->I_yz = (ma_I_yz + mb_I_yz) - M_tot * ma->CoM[1] * ma->CoM[2];
+#endif
 }
 
 /**
@@ -124,50 +220,3 @@ void multipole_addparts(struct multipole *m, struct gpart *p, int N) {
    */
   /* #endif */
 }
-
-/**
-* @brief Reset the data of a #multipole.
-*
-* @param m The #multipole.
-*/
-void multipole_reset(struct multipole *m) {
-
-  /* Just bzero the struct. */
-  bzero(m, sizeof(struct multipole));
-}
-
-/**
-* @brief Init a multipole from a set of particles.
-*
-* @param m The #multipole.
-* @param parts The #gpart.
-* @param N The number of particles.
-*/
-void multipole_init(struct multipole *m, const struct gpart *gparts,
-                    int gcount) {
-
-#if multipole_order != 1
-#error "Multipoles of order >1 not yet implemented."
-#endif
-
-  /* Zero everything */
-  multipole_reset(m);
-
-  float mass = 0.0f;
-  double x[3] = {0.0, 0.0, 0.0};
-
-  /* Collect the particle data. */
-  for (int k = 0; k < gcount; k++) {
-    const float w = gparts[k].mass;
-    mass += w;
-    x[0] += gparts[k].x[0] * w;
-    x[1] += gparts[k].x[1] * w;
-    x[2] += gparts[k].x[2] * w;
-  }
-
-  /* Store the data on the multipole. */
-  m->mass = mass;
-  m->CoM[0] = x[0] / mass;
-  m->CoM[1] = x[1] / mass;
-  m->CoM[2] = x[2] / mass;
-}

src/multipole.h

@@ -38,11 +38,11 @@ struct multipole {
   /* Multipole mass */
   float mass;
 
-  /* /\* Acceleration on this multipole. *\/ */
-  /* float a[3]; */
-
-  /* /\* Multipole coefficients. *\/ */
-  /* float coeffs[multipole_order * multipole_order]; */
+#if multipole_order >= 2
+  /* Quadrupole terms */
+  float I_xx, I_yy, I_zz;
+  float I_xy, I_xz, I_yz;
+#endif
 };
 
 /* Multipole function prototypes. */