Quadrupole construction is now entirely recursive.

examples/test_fmm.c

@@ -216,12 +216,12 @@ void comp_com(struct cell *c) {
       com[2] += cp->new.com[2] * cp_mass;
       mass += cp_mass;
 
-      I_xx += cp->new.I_xx;
-      I_yy += cp->new.I_yy;
-      I_zz += cp->new.I_zz;
-      I_xy += cp->new.I_xy;
-      I_xz += cp->new.I_xz;
-      I_yz += cp->new.I_yz;
+      I_xx += cp->new.I_xx + cp_mass * cp->new.com[0] * cp->new.com[0];
+      I_yy += cp->new.I_yy + cp_mass * cp->new.com[1] * cp->new.com[1];
+      I_zz += cp->new.I_zz + cp_mass * cp->new.com[2] * cp->new.com[2];
+      I_xy += cp->new.I_xy + cp_mass * cp->new.com[0] * cp->new.com[1];
+      I_xz += cp->new.I_xz + cp_mass * cp->new.com[0] * cp->new.com[2];
+      I_yz += cp->new.I_yz + cp_mass * cp->new.com[1] * cp->new.com[2];
     }
 
     /* Otherwise, collect the multipole from the particles. */
@@ -249,12 +249,12 @@ void comp_com(struct cell *c) {
   c->new.com[0] = com[0] * imass;
   c->new.com[1] = com[1] * imass;
   c->new.com[2] = com[2] * imass;
-  c->new.I_xx = I_xx;
-  c->new.I_yy = I_yy;
-  c->new.I_zz = I_zz;
-  c->new.I_xy = I_xy;
-  c->new.I_xz = I_xz;
-  c->new.I_yz = I_yz;
+  c->new.I_xx = I_xx - imass * com[0] * com[0];
+  c->new.I_yy = I_yy - imass * com[1] * com[1];
+  c->new.I_zz = I_zz - imass * com[2] * com[2];
+  c->new.I_xy = I_xy - imass * com[0] * com[1];
+  c->new.I_xz = I_xz - imass * com[0] * com[2];
+  c->new.I_yz = I_yz - imass * com[1] * com[2];
 
   /* Reset COM acceleration */
   c->a[0] = 0.;
@@ -1460,6 +1460,30 @@ void test_bh(int N, int nr_threads, int runs, char *fileName) {
   message("[check] root CoMx= %f %f", root->legacy.com[0], root->new.com[0]);
   message("[check] root CoMy= %f %f", root->legacy.com[1], root->new.com[1]);
   message("[check] root CoMz= %f %f", root->legacy.com[2], root->new.com[2]);
+  
+  /* message("[check]     | %f %f %f |", root->new.I_xx, root->new.I_xy, root->new.I_xz); */
+  /* message("[check] I = | %f %f %f |", root->new.I_xy, root->new.I_yy, root->new.I_yz); */
+  /* message("[check]     | %f %f %f |", root->new.I_xz, root->new.I_yz, root->new.I_zz); */
+
+  /* double I_xx = 0., I_yy = 0., I_zz = 0.,  I_xy = 0., I_xz = 0., I_yz = 0.; */
+  /* for (i = 0; i < N; ++i) { */
+  /*   I_xx += root->parts[i].x[0] * root->parts[i].x[0] * root->parts[i].mass; */
+  /*   I_yy += root->parts[i].x[1] * root->parts[i].x[1] * root->parts[i].mass; */
+  /*   I_zz += root->parts[i].x[2] * root->parts[i].x[2] * root->parts[i].mass; */
+  /*   I_xy += root->parts[i].x[0] * root->parts[i].x[1] * root->parts[i].mass; */
+  /*   I_xz += root->parts[i].x[0] * root->parts[i].x[2] * root->parts[i].mass; */
+  /*   I_yz += root->parts[i].x[1] * root->parts[i].x[2] * root->parts[i].mass; */
+  /* } */
+  /* I_xx = I_xx - root->new.mass * root->new.com[0] * root->new.com[0]; */
+  /* I_yy = I_yy - root->new.mass * root->new.com[1] * root->new.com[1]; */
+  /* I_zz = I_zz - root->new.mass * root->new.com[2] * root->new.com[2]; */
+  /* I_xy = I_xy - root->new.mass * root->new.com[0] * root->new.com[1]; */
+  /* I_xz = I_xz - root->new.mass * root->new.com[0] * root->new.com[2]; */
+  /* I_yz = I_yz - root->new.mass * root->new.com[1] * root->new.com[2]; */
+  /* message("[check]     | %f %f %f |", I_xx, I_xy, I_xz); */
+  /* message("[check] I = | %f %f %f |", I_xy, I_yy, I_yz); */
+  /* message("[check]     | %f %f %f |", I_xz, I_yz, I_zz); */
+
 #if ICHECK >= 0
   for (i = 0; i < N; ++i)
     if (root->parts[i].id == ICHECK)