         if (flag && nodes[offset + 6] == 0) nodes[offset + 6] = child + 6;
         if (flag && nodes[offset + 7] == 0) nodes[offset + 7] = child + 7;
         __syncthreads();
-#endif
 
         n_offset += n_total;    //Increase the offset in the array by the number of newly added nodes
 


        n_approx += n_total;

        while (n_approx >= DIM) 
        {
          n_approx -= DIM;
          const int address      = (approx[n_approx + tid] << 1) + approx[n_approx + tid];
#ifndef TEXTURES
          const float4 monopole  = multipole_data[address    ];
#if 0
          float4 octopole0 = multipole_data[address + 1];
          float4 octopole1 = multipole_data[address + 2];
#endif
#else
          const float4 monopole  = tex1Dfetch(texMultipole, address);
#if 0
          float4 octopole0 = tex1Dfetch(texMultipole, address + 1);
          float4 octopole1 = tex1Dfetch(texMultipole, address + 2);
#endif
#endif

          node_mon0[tid] = monopole.w;
          node_mon1[tid] = make_float3(monopole.x,  monopole.y,  monopole.z);
          __syncthreads();

#if 0
          const float f_dm   = 0.0f;
          const float f_star = 1.0f
            const float darkMatterMass = f_dm   * octopole1.w;
          /* eg: we need to be careful with the line below to avoid truncation error due to 
             subtraction of two large numbers, monopole.w and darkMatterMass both could be
             very large.
             Instead, we can use octopole1.w to be stellar mass, and DM mass to be 
             monopole.w, then we add the two together to get total mass, but this will
             require more changes to the kernel */
          const float    stellarMass = f_star * (monopole.w - darkMatterMass);
          const float hinv = 1.0f/hi;   /* eg: this can be precomputing to avoid division */
          density += interact(
              make_float3(pos_i.x, pos_i.y, pos_i.z), h, hinv,
              make_float3(monopole.x, monople.y, monopole.z), darkMatterMass + stellarMass);
          /* eg: the interact function still calls sqrtf(f), which invloves 1 div and 1 rsqrtf,
             so ideally we would like to take advantage of rsqrtf in add_acc, and then we only
             do 1 div */
#endif


#if 1
#pragma unroll 16

  __syncthreads();

  if (warpId > 0) val += shdata[warpId - 1];

  __syncthreads();
  if ((1 << BLOCKDIM2) - 1 == idx) shdata[0] = val;
  __syncthreads();
