Simulation: Rename `physics` directory to `simulation`

Function names will be updated in a separate commit.

This will be the place for the new particle system and other
code related to the Simulation data block. We don't want
to have all that code in blenkernel.

Approved by brecht.

1 files changed, 0 insertions, 1274 deletions

diff --git a/source/blender/physics/intern/hair_volume.cpp b/source/blender/physics/intern/hair_volume.cpp
deleted file mode 100644
index 1764d0a910c..00000000000
--- a/source/blender/physics/intern/hair_volume.cpp
+++ /dev/null
@@ -1,1274 +0,0 @@
-/*
- * This program is free software; you can redistribute it and/or
- * modify it under the terms of the GNU General Public License
- * as published by the Free Software Foundation; either version 2
- * of the License, or (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program; if not, write to the Free Software Foundation,
- * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
- *
- * The Original Code is Copyright (C) Blender Foundation
- * All rights reserved.
- */
-
-/** \file
- * \ingroup bph
- */
-
-#include "MEM_guardedalloc.h"
-
-#include "BLI_math.h"
-#include "BLI_utildefines.h"
-
-#include "DNA_texture_types.h"
-
-#include "BKE_effect.h"
-
-#include "eigen_utils.h"
-#include "implicit.h"
-
-/* ================ Volumetric Hair Interaction ================
- * adapted from
- *
- * Volumetric Methods for Simulation and Rendering of Hair
- *     (Petrovic, Henne, Anderson, Pixar Technical Memo #06-08, Pixar Animation Studios)
- *
- * as well as
- *
- * "Detail Preserving Continuum Simulation of Straight Hair"
- *     (McAdams, Selle 2009)
- */
-
-/* Note about array indexing:
- * Generally the arrays here are one-dimensional.
- * The relation between 3D indices and the array offset is
- *   offset = x + res_x * y + res_x * res_y * z
- */
-
-static float I[3][3] = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
-
-BLI_INLINE int floor_int(float value)
-{
-  return value > 0.0f ? (int)value : ((int)value) - 1;
-}
-
-BLI_INLINE float floor_mod(float value)
-{
-  return value - floorf(value);
-}
-
-BLI_INLINE int hair_grid_size(const int res[3])
-{
-  return res[0] * res[1] * res[2];
-}
-
-typedef struct HairGridVert {
-  int samples;
-  float velocity[3];
-  float density;
-
-  float velocity_smooth[3];
-} HairGridVert;
-
-typedef struct HairGrid {
-  HairGridVert *verts;
-  int res[3];
-  float gmin[3], gmax[3];
-  float cellsize, inv_cellsize;
-} HairGrid;
-
-#define HAIR_GRID_INDEX_AXIS(vec, res, gmin, scale, axis) \
-  (min_ii(max_ii((int)((vec[axis] - gmin[axis]) * scale), 0), res[axis] - 2))
-
-BLI_INLINE int hair_grid_offset(const float vec[3],
-                                const int res[3],
-                                const float gmin[3],
-                                float scale)
-{
-  int i, j, k;
-  i = HAIR_GRID_INDEX_AXIS(vec, res, gmin, scale, 0);
-  j = HAIR_GRID_INDEX_AXIS(vec, res, gmin, scale, 1);
-  k = HAIR_GRID_INDEX_AXIS(vec, res, gmin, scale, 2);
-  return i + (j + k * res[1]) * res[0];
-}
-
-BLI_INLINE int hair_grid_interp_weights(
-    const int res[3], const float gmin[3], float scale, const float vec[3], float uvw[3])
-{
-  int i, j, k, offset;
-
-  i = HAIR_GRID_INDEX_AXIS(vec, res, gmin, scale, 0);
-  j = HAIR_GRID_INDEX_AXIS(vec, res, gmin, scale, 1);
-  k = HAIR_GRID_INDEX_AXIS(vec, res, gmin, scale, 2);
-  offset = i + (j + k * res[1]) * res[0];
-
-  uvw[0] = (vec[0] - gmin[0]) * scale - (float)i;
-  uvw[1] = (vec[1] - gmin[1]) * scale - (float)j;
-  uvw[2] = (vec[2] - gmin[2]) * scale - (float)k;
-
-  //  BLI_assert(0.0f <= uvw[0] && uvw[0] <= 1.0001f);
-  //  BLI_assert(0.0f <= uvw[1] && uvw[1] <= 1.0001f);
-  //  BLI_assert(0.0f <= uvw[2] && uvw[2] <= 1.0001f);
-
-  return offset;
-}
-
-BLI_INLINE void hair_grid_interpolate(const HairGridVert *grid,
-                                      const int res[3],
-                                      const float gmin[3],
-                                      float scale,
-                                      const float vec[3],
-                                      float *density,
-                                      float velocity[3],
-                                      float vel_smooth[3],
-                                      float density_gradient[3],
-                                      float velocity_gradient[3][3])
-{
-  HairGridVert data[8];
-  float uvw[3], muvw[3];
-  int res2 = res[1] * res[0];
-  int offset;
-
-  offset = hair_grid_interp_weights(res, gmin, scale, vec, uvw);
-  muvw[0] = 1.0f - uvw[0];
-  muvw[1] = 1.0f - uvw[1];
-  muvw[2] = 1.0f - uvw[2];
-
-  data[0] = grid[offset];
-  data[1] = grid[offset + 1];
-  data[2] = grid[offset + res[0]];
-  data[3] = grid[offset + res[0] + 1];
-  data[4] = grid[offset + res2];
-  data[5] = grid[offset + res2 + 1];
-  data[6] = grid[offset + res2 + res[0]];
-  data[7] = grid[offset + res2 + res[0] + 1];
-
-  if (density) {
-    *density = muvw[2] * (muvw[1] * (muvw[0] * data[0].density + uvw[0] * data[1].density) +
-                          uvw[1] * (muvw[0] * data[2].density + uvw[0] * data[3].density)) +
-               uvw[2] * (muvw[1] * (muvw[0] * data[4].density + uvw[0] * data[5].density) +
-                         uvw[1] * (muvw[0] * data[6].density + uvw[0] * data[7].density));
-  }
-
-  if (velocity) {
-    int k;
-    for (k = 0; k < 3; k++) {
-      velocity[k] = muvw[2] *
-                        (muvw[1] * (muvw[0] * data[0].velocity[k] + uvw[0] * data[1].velocity[k]) +
-                         uvw[1] * (muvw[0] * data[2].velocity[k] + uvw[0] * data[3].velocity[k])) +
-                    uvw[2] *
-                        (muvw[1] * (muvw[0] * data[4].velocity[k] + uvw[0] * data[5].velocity[k]) +
-                         uvw[1] * (muvw[0] * data[6].velocity[k] + uvw[0] * data[7].velocity[k]));
-    }
-  }
-
-  if (vel_smooth) {
-    int k;
-    for (k = 0; k < 3; k++) {
-      vel_smooth[k] = muvw[2] * (muvw[1] * (muvw[0] * data[0].velocity_smooth[k] +
-                                            uvw[0] * data[1].velocity_smooth[k]) +
-                                 uvw[1] * (muvw[0] * data[2].velocity_smooth[k] +
-                                           uvw[0] * data[3].velocity_smooth[k])) +
-                      uvw[2] * (muvw[1] * (muvw[0] * data[4].velocity_smooth[k] +
-                                           uvw[0] * data[5].velocity_smooth[k]) +
-                                uvw[1] * (muvw[0] * data[6].velocity_smooth[k] +
-                                          uvw[0] * data[7].velocity_smooth[k]));
-    }
-  }
-
-  if (density_gradient) {
-    density_gradient[0] = muvw[1] * muvw[2] * (data[0].density - data[1].density) +
-                          uvw[1] * muvw[2] * (data[2].density - data[3].density) +
-                          muvw[1] * uvw[2] * (data[4].density - data[5].density) +
-                          uvw[1] * uvw[2] * (data[6].density - data[7].density);
-
-    density_gradient[1] = muvw[2] * muvw[0] * (data[0].density - data[2].density) +
-                          uvw[2] * muvw[0] * (data[4].density - data[6].density) +
-                          muvw[2] * uvw[0] * (data[1].density - data[3].density) +
-                          uvw[2] * uvw[0] * (data[5].density - data[7].density);
-
-    density_gradient[2] = muvw[2] * muvw[0] * (data[0].density - data[4].density) +
-                          uvw[2] * muvw[0] * (data[1].density - data[5].density) +
-                          muvw[2] * uvw[0] * (data[2].density - data[6].density) +
-                          uvw[2] * uvw[0] * (data[3].density - data[7].density);
-  }
-
-  if (velocity_gradient) {
-    /* XXX TODO */
-    zero_m3(velocity_gradient);
-  }
-}
-
-void BPH_hair_volume_vertex_grid_forces(HairGrid *grid,
-                                        const float x[3],
-                                        const float v[3],
-                                        float smoothfac,
-                                        float pressurefac,
-                                        float minpressure,
-                                        float f[3],
-                                        float dfdx[3][3],
-                                        float dfdv[3][3])
-{
-  float gdensity, gvelocity[3], ggrad[3], gvelgrad[3][3], gradlen;
-
-  hair_grid_interpolate(grid->verts,
-                        grid->res,
-                        grid->gmin,
-                        grid->inv_cellsize,
-                        x,
-                        &gdensity,
-                        gvelocity,
-                        NULL,
-                        ggrad,
-                        gvelgrad);
-
-  zero_v3(f);
-  sub_v3_v3(gvelocity, v);
-  mul_v3_v3fl(f, gvelocity, smoothfac);
-
-  gradlen = normalize_v3(ggrad) - minpressure;
-  if (gradlen > 0.0f) {
-    mul_v3_fl(ggrad, gradlen);
-    madd_v3_v3fl(f, ggrad, pressurefac);
-  }
-
-  zero_m3(dfdx);
-
-  sub_m3_m3m3(dfdv, gvelgrad, I);
-  mul_m3_fl(dfdv, smoothfac);
-}
-
-void BPH_hair_volume_grid_interpolate(HairGrid *grid,
-                                      const float x[3],
-                                      float *density,
-                                      float velocity[3],
-                                      float velocity_smooth[3],
-                                      float density_gradient[3],
-                                      float velocity_gradient[3][3])
-{
-  hair_grid_interpolate(grid->verts,
-                        grid->res,
-                        grid->gmin,
-                        grid->inv_cellsize,
-                        x,
-                        density,
-                        velocity,
-                        velocity_smooth,
-                        density_gradient,
-                        velocity_gradient);
-}
-
-void BPH_hair_volume_grid_velocity(
-    HairGrid *grid, const float x[3], const float v[3], float fluid_factor, float r_v[3])
-{
-  float gdensity, gvelocity[3], gvel_smooth[3], ggrad[3], gvelgrad[3][3];
-  float v_pic[3], v_flip[3];
-
-  hair_grid_interpolate(grid->verts,
-                        grid->res,
-                        grid->gmin,
-                        grid->inv_cellsize,
-                        x,
-                        &gdensity,
-                        gvelocity,
-                        gvel_smooth,
-                        ggrad,
-                        gvelgrad);
-
-  /* velocity according to PIC method (Particle-in-Cell) */
-  copy_v3_v3(v_pic, gvel_smooth);
-
-  /* velocity according to FLIP method (Fluid-Implicit-Particle) */
-  sub_v3_v3v3(v_flip, gvel_smooth, gvelocity);
-  add_v3_v3(v_flip, v);
-
-  interp_v3_v3v3(r_v, v_pic, v_flip, fluid_factor);
-}
-
-void BPH_hair_volume_grid_clear(HairGrid *grid)
-{
-  const int size = hair_grid_size(grid->res);
-  int i;
-  for (i = 0; i < size; i++) {
-    zero_v3(grid->verts[i].velocity);
-    zero_v3(grid->verts[i].velocity_smooth);
-    grid->verts[i].density = 0.0f;
-    grid->verts[i].samples = 0;
-  }
-}
-
-BLI_INLINE bool hair_grid_point_valid(const float vec[3], const float gmin[3], const float gmax[3])
-{
-  return !(vec[0] < gmin[0] || vec[1] < gmin[1] || vec[2] < gmin[2] || vec[0] > gmax[0] ||
-           vec[1] > gmax[1] || vec[2] > gmax[2]);
-}
-
-BLI_INLINE float dist_tent_v3f3(const float a[3], float x, float y, float z)
-{
-  float w = (1.0f - fabsf(a[0] - x)) * (1.0f - fabsf(a[1] - y)) * (1.0f - fabsf(a[2] - z));
-  return w;
-}
-
-BLI_INLINE float weights_sum(const float weights[8])
-{
-  float totweight = 0.0f;
-  int i;
-  for (i = 0; i < 8; i++) {
-    totweight += weights[i];
-  }
-  return totweight;
-}
-
-/* returns the grid array offset as well to avoid redundant calculation */
-BLI_INLINE int hair_grid_weights(
-    const int res[3], const float gmin[3], float scale, const float vec[3], float weights[8])
-{
-  int i, j, k, offset;
-  float uvw[3];
-
-  i = HAIR_GRID_INDEX_AXIS(vec, res, gmin, scale, 0);
-  j = HAIR_GRID_INDEX_AXIS(vec, res, gmin, scale, 1);
-  k = HAIR_GRID_INDEX_AXIS(vec, res, gmin, scale, 2);
-  offset = i + (j + k * res[1]) * res[0];
-
-  uvw[0] = (vec[0] - gmin[0]) * scale;
-  uvw[1] = (vec[1] - gmin[1]) * scale;
-  uvw[2] = (vec[2] - gmin[2]) * scale;
-
-  weights[0] = dist_tent_v3f3(uvw, (float)i, (float)j, (float)k);
-  weights[1] = dist_tent_v3f3(uvw, (float)(i + 1), (float)j, (float)k);
-  weights[2] = dist_tent_v3f3(uvw, (float)i, (float)(j + 1), (float)k);
-  weights[3] = dist_tent_v3f3(uvw, (float)(i + 1), (float)(j + 1), (float)k);
-  weights[4] = dist_tent_v3f3(uvw, (float)i, (float)j, (float)(k + 1));
-  weights[5] = dist_tent_v3f3(uvw, (float)(i + 1), (float)j, (float)(k + 1));
-  weights[6] = dist_tent_v3f3(uvw, (float)i, (float)(j + 1), (float)(k + 1));
-  weights[7] = dist_tent_v3f3(uvw, (float)(i + 1), (float)(j + 1), (float)(k + 1));
-
-  //  BLI_assert(fabsf(weights_sum(weights) - 1.0f) < 0.0001f);
-
-  return offset;
-}
-
-BLI_INLINE void grid_to_world(HairGrid *grid, float vecw[3], const float vec[3])
-{
-  copy_v3_v3(vecw, vec);
-  mul_v3_fl(vecw, grid->cellsize);
-  add_v3_v3(vecw, grid->gmin);
-}
-
-void BPH_hair_volume_add_vertex(HairGrid *grid, const float x[3], const float v[3])
-{
-  const int res[3] = {grid->res[0], grid->res[1], grid->res[2]};
-  float weights[8];
-  int di, dj, dk;
-  int offset;
-
-  if (!hair_grid_point_valid(x, grid->gmin, grid->gmax)) {
-    return;
-  }
-
-  offset = hair_grid_weights(res, grid->gmin, grid->inv_cellsize, x, weights);
-
-  for (di = 0; di < 2; di++) {
-    for (dj = 0; dj < 2; dj++) {
-      for (dk = 0; dk < 2; dk++) {
-        int voffset = offset + di + (dj + dk * res[1]) * res[0];
-        int iw = di + dj * 2 + dk * 4;
-
-        grid->verts[voffset].density += weights[iw];
-        madd_v3_v3fl(grid->verts[voffset].velocity, v, weights[iw]);
-      }
-    }
-  }
-}
-
-#if 0
-BLI_INLINE void hair_volume_eval_grid_vertex(HairGridVert *vert,
-                                             const float loc[3],
-                                             float radius,
-                                             float dist_scale,
-                                             const float x2[3],
-                                             const float v2[3],
-                                             const float x3[3],
-                                             const float v3[3])
-{
-  float closest[3], lambda, dist, weight;
-
-  lambda = closest_to_line_v3(closest, loc, x2, x3);
-  dist = len_v3v3(closest, loc);
-
-  weight = (radius - dist) * dist_scale;
-
-  if (weight > 0.0f) {
-    float vel[3];
-
-    interp_v3_v3v3(vel, v2, v3, lambda);
-    madd_v3_v3fl(vert->velocity, vel, weight);
-    vert->density += weight;
-    vert->samples += 1;
-  }
-}
-
-BLI_INLINE int major_axis_v3(const float v[3])
-{
-  const float a = fabsf(v[0]);
-  const float b = fabsf(v[1]);
-  const float c = fabsf(v[2]);
-  return a > b ? (a > c ? 0 : 2) : (b > c ? 1 : 2);
-}
-
-BLI_INLINE void hair_volume_add_segment_2D(HairGrid *grid,
-                                           const float UNUSED(x1[3]),
-                                           const float UNUSED(v1[3]),
-                                           const float x2[3],
-                                           const float v2[3],
-                                           const float x3[3],
-                                           const float v3[3],
-                                           const float UNUSED(x4[3]),
-                                           const float UNUSED(v4[3]),
-                                           const float UNUSED(dir1[3]),
-                                           const float dir2[3],
-                                           const float UNUSED(dir3[3]),
-                                           int resj,
-                                           int resk,
-                                           int jmin,
-                                           int jmax,
-                                           int kmin,
-                                           int kmax,
-                                           HairGridVert *vert,
-                                           int stride_j,
-                                           int stride_k,
-                                           const float loc[3],
-                                           int axis_j,
-                                           int axis_k,
-                                           int debug_i)
-{
-  const float radius = 1.5f;
-  const float dist_scale = grid->inv_cellsize;
-
-  int j, k;
-
-  /* boundary checks to be safe */
-  CLAMP_MIN(jmin, 0);
-  CLAMP_MAX(jmax, resj - 1);
-  CLAMP_MIN(kmin, 0);
-  CLAMP_MAX(kmax, resk - 1);
-
-  HairGridVert *vert_j = vert + jmin * stride_j;
-  float loc_j[3] = {loc[0], loc[1], loc[2]};
-  loc_j[axis_j] += (float)jmin;
-  for (j = jmin; j <= jmax; j++, vert_j += stride_j, loc_j[axis_j] += 1.0f) {
-
-    HairGridVert *vert_k = vert_j + kmin * stride_k;
-    float loc_k[3] = {loc_j[0], loc_j[1], loc_j[2]};
-    loc_k[axis_k] += (float)kmin;
-    for (k = kmin; k <= kmax; k++, vert_k += stride_k, loc_k[axis_k] += 1.0f) {
-
-      hair_volume_eval_grid_vertex(vert_k, loc_k, radius, dist_scale, x2, v2, x3, v3);
-
-#  if 0
-      {
-        float wloc[3], x2w[3], x3w[3];
-        grid_to_world(grid, wloc, loc_k);
-        grid_to_world(grid, x2w, x2);
-        grid_to_world(grid, x3w, x3);
-
-        if (vert_k->samples > 0) {
-          BKE_sim_debug_data_add_circle(wloc, 0.01f, 1.0, 1.0, 0.3, "grid", 2525, debug_i, j, k);
-        }
-
-        if (grid->debug_value) {
-          BKE_sim_debug_data_add_dot(wloc, 1, 0, 0, "grid", 93, debug_i, j, k);
-          BKE_sim_debug_data_add_dot(x2w, 0.1, 0.1, 0.7, "grid", 649, debug_i, j, k);
-          BKE_sim_debug_data_add_line(wloc, x2w, 0.3, 0.8, 0.3, "grid", 253, debug_i, j, k);
-          BKE_sim_debug_data_add_line(wloc, x3w, 0.8, 0.3, 0.3, "grid", 254, debug_i, j, k);
-          // BKE_sim_debug_data_add_circle(
-          //     x2w, len_v3v3(wloc, x2w), 0.2, 0.7, 0.2,
-          //     "grid", 255, i, j, k);
-        }
-      }
-#  endif
-    }
-  }
-}
-
-/* Uses a variation of Bresenham's algorithm for rasterizing a 3D grid with a line segment.
- *
- * The radius of influence around a segment is assumed to be at most 2*cellsize,
- * i.e. only cells containing the segment and their direct neighbors are examined.
- */
-void BPH_hair_volume_add_segment(HairGrid *grid,
-                                 const float x1[3],
-                                 const float v1[3],
-                                 const float x2[3],
-                                 const float v2[3],
-                                 const float x3[3],
-                                 const float v3[3],
-                                 const float x4[3],
-                                 const float v4[3],
-                                 const float dir1[3],
-                                 const float dir2[3],
-                                 const float dir3[3])
-{
-  const int res[3] = {grid->res[0], grid->res[1], grid->res[2]};
-
-  /* find the primary direction from the major axis of the direction vector */
-  const int axis0 = major_axis_v3(dir2);
-  const int axis1 = (axis0 + 1) % 3;
-  const int axis2 = (axis0 + 2) % 3;
-
-  /* vertex buffer offset factors along cardinal axes */
-  const int strides[3] = {1, res[0], res[0] * res[1]};
-  const int stride0 = strides[axis0];
-  const int stride1 = strides[axis1];
-  const int stride2 = strides[axis2];
-
-  /* increment of secondary directions per step in the primary direction
-   * note: we always go in the positive direction along axis0, so the sign can be inverted
-   */
-  const float inc1 = dir2[axis1] / dir2[axis0];
-  const float inc2 = dir2[axis2] / dir2[axis0];
-
-  /* start/end points, so increment along axis0 is always positive */
-  const float *start = x2[axis0] < x3[axis0] ? x2 : x3;
-  const float *end = x2[axis0] < x3[axis0] ? x3 : x2;
-  const float start0 = start[axis0], start1 = start[axis1], start2 = start[axis2];
-  const float end0 = end[axis0];
-
-  /* range along primary direction */
-  const int imin = max_ii(floor_int(start[axis0]) - 1, 0);
-  const int imax = min_ii(floor_int(end[axis0]) + 2, res[axis0] - 1);
-
-  float h = 0.0f;
-  HairGridVert *vert0;
-  float loc0[3];
-  int j0, k0, j0_prev, k0_prev;
-  int i;
-
-  for (i = imin; i <= imax; i++) {
-    float shift1, shift2; /* fraction of a full cell shift [0.0, 1.0) */
-    int jmin, jmax, kmin, kmax;
-
-    h = CLAMPIS((float)i, start0, end0);
-
-    shift1 = start1 + (h - start0) * inc1;
-    shift2 = start2 + (h - start0) * inc2;
-
-    j0_prev = j0;
-    j0 = floor_int(shift1);
-
-    k0_prev = k0;
-    k0 = floor_int(shift2);
-
-    if (i > imin) {
-      jmin = min_ii(j0, j0_prev);
-      jmax = max_ii(j0, j0_prev);
-      kmin = min_ii(k0, k0_prev);
-      kmax = max_ii(k0, k0_prev);
-    }
-    else {
-      jmin = jmax = j0;
-      kmin = kmax = k0;
-    }
-
-    vert0 = grid->verts + i * stride0;
-    loc0[axis0] = (float)i;
-    loc0[axis1] = 0.0f;
-    loc0[axis2] = 0.0f;
-
-    hair_volume_add_segment_2D(grid,
-                               x1,
-                               v1,
-                               x2,
-                               v2,
-                               x3,
-                               v3,
-                               x4,
-                               v4,
-                               dir1,
-                               dir2,
-                               dir3,
-                               res[axis1],
-                               res[axis2],
-                               jmin - 1,
-                               jmax + 2,
-                               kmin - 1,
-                               kmax + 2,
-                               vert0,
-                               stride1,
-                               stride2,
-                               loc0,
-                               axis1,
-                               axis2,
-                               i);
-  }
-}
-#else
-BLI_INLINE void hair_volume_eval_grid_vertex_sample(HairGridVert *vert,
-                                                    const float loc[3],
-                                                    float radius,
-                                                    float dist_scale,
-                                                    const float x[3],
-                                                    const float v[3])
-{
-  float dist, weight;
-
-  dist = len_v3v3(x, loc);
-
-  weight = (radius - dist) * dist_scale;
-
-  if (weight > 0.0f) {
-    madd_v3_v3fl(vert->velocity, v, weight);
-    vert->density += weight;
-    vert->samples += 1;
-  }
-}
-
-/* XXX simplified test implementation using a series of discrete sample along the segment,
- * instead of finding the closest point for all affected grid vertices.
- */
-void BPH_hair_volume_add_segment(HairGrid *grid,
-                                 const float UNUSED(x1[3]),
-                                 const float UNUSED(v1[3]),
-                                 const float x2[3],
-                                 const float v2[3],
-                                 const float x3[3],
-                                 const float v3[3],
-                                 const float UNUSED(x4[3]),
-                                 const float UNUSED(v4[3]),
-                                 const float UNUSED(dir1[3]),
-                                 const float UNUSED(dir2[3]),
-                                 const float UNUSED(dir3[3]))
-{
-  const float radius = 1.5f;
-  const float dist_scale = grid->inv_cellsize;
-
-  const int res[3] = {grid->res[0], grid->res[1], grid->res[2]};
-  const int stride[3] = {1, res[0], res[0] * res[1]};
-  const int num_samples = 10;
-
-  int s;
-
-  for (s = 0; s < num_samples; s++) {
-    float x[3], v[3];
-    int i, j, k;
-
-    float f = (float)s / (float)(num_samples - 1);
-    interp_v3_v3v3(x, x2, x3, f);
-    interp_v3_v3v3(v, v2, v3, f);
-
-    int imin = max_ii(floor_int(x[0]) - 2, 0);
-    int imax = min_ii(floor_int(x[0]) + 2, res[0] - 1);
-    int jmin = max_ii(floor_int(x[1]) - 2, 0);
-    int jmax = min_ii(floor_int(x[1]) + 2, res[1] - 1);
-    int kmin = max_ii(floor_int(x[2]) - 2, 0);
-    int kmax = min_ii(floor_int(x[2]) + 2, res[2] - 1);
-
-    for (k = kmin; k <= kmax; k++) {
-      for (j = jmin; j <= jmax; j++) {
-        for (i = imin; i <= imax; i++) {
-          float loc[3] = {(float)i, (float)j, (float)k};
-          HairGridVert *vert = grid->verts + i * stride[0] + j * stride[1] + k * stride[2];
-
-          hair_volume_eval_grid_vertex_sample(vert, loc, radius, dist_scale, x, v);
-        }
-      }
-    }
-  }
-}
-#endif
-
-void BPH_hair_volume_normalize_vertex_grid(HairGrid *grid)
-{
-  int i, size = hair_grid_size(grid->res);
-  /* divide velocity with density */
-  for (i = 0; i < size; i++) {
-    float density = grid->verts[i].density;
-    if (density > 0.0f) {
-      mul_v3_fl(grid->verts[i].velocity, 1.0f / density);
-    }
-  }
-}
-
-/* Cells with density below this are considered empty. */
-static const float density_threshold = 0.001f;
-
-/* Contribution of target density pressure to the laplacian in the pressure poisson equation.
- * This is based on the model found in
- * "Two-way Coupled SPH and Particle Level Set Fluid Simulation" (Losasso et al., 2008)
- */
-BLI_INLINE float hair_volume_density_divergence(float density,
-                                                float target_density,
-                                                float strength)
-{
-  if (density > density_threshold && density > target_density) {
-    return strength * logf(target_density / density);
-  }
-  else {
-    return 0.0f;
-  }
-}
-
-bool BPH_hair_volume_solve_divergence(HairGrid *grid,
-                                      float /*dt*/,
-                                      float target_density,
-                                      float target_strength)
-{
-  const float flowfac = grid->cellsize;
-  const float inv_flowfac = 1.0f / grid->cellsize;
-
-  /*const int num_cells = hair_grid_size(grid->res);*/
-  const int res[3] = {grid->res[0], grid->res[1], grid->res[2]};
-  const int resA[3] = {grid->res[0] + 2, grid->res[1] + 2, grid->res[2] + 2};
-
-  const int stride0 = 1;
-  const int stride1 = grid->res[0];
-  const int stride2 = grid->res[1] * grid->res[0];
-  const int strideA0 = 1;
-  const int strideA1 = grid->res[0] + 2;
-  const int strideA2 = (grid->res[1] + 2) * (grid->res[0] + 2);
-
-  const int num_cells = res[0] * res[1] * res[2];
-  const int num_cellsA = (res[0] + 2) * (res[1] + 2) * (res[2] + 2);
-
-  HairGridVert *vert_start = grid->verts - (stride0 + stride1 + stride2);
-  HairGridVert *vert;
-  int i, j, k;
-
-#define MARGIN_i0 (i < 1)
-#define MARGIN_j0 (j < 1)
-#define MARGIN_k0 (k < 1)
-#define MARGIN_i1 (i >= resA[0] - 1)
-#define MARGIN_j1 (j >= resA[1] - 1)
-#define MARGIN_k1 (k >= resA[2] - 1)
-
-#define NEIGHBOR_MARGIN_i0 (i < 2)
-#define NEIGHBOR_MARGIN_j0 (j < 2)
-#define NEIGHBOR_MARGIN_k0 (k < 2)
-#define NEIGHBOR_MARGIN_i1 (i >= resA[0] - 2)
-#define NEIGHBOR_MARGIN_j1 (j >= resA[1] - 2)
-#define NEIGHBOR_MARGIN_k1 (k >= resA[2] - 2)
-
-  BLI_assert(num_cells >= 1);
-
-  /* Calculate divergence */
-  lVector B(num_cellsA);
-  for (k = 0; k < resA[2]; k++) {
-    for (j = 0; j < resA[1]; j++) {
-      for (i = 0; i < resA[0]; i++) {
-        int u = i * strideA0 + j * strideA1 + k * strideA2;
-        bool is_margin = MARGIN_i0 || MARGIN_i1 || MARGIN_j0 || MARGIN_j1 || MARGIN_k0 ||
-                         MARGIN_k1;
-
-        if (is_margin) {
-          B[u] = 0.0f;
-          continue;
-        }
-
-        vert = vert_start + i * stride0 + j * stride1 + k * stride2;
-
-        const float *v0 = vert->velocity;
-        float dx = 0.0f, dy = 0.0f, dz = 0.0f;
-        if (!NEIGHBOR_MARGIN_i0) {
-          dx += v0[0] - (vert - stride0)->velocity[0];
-        }
-        if (!NEIGHBOR_MARGIN_i1) {
-          dx += (vert + stride0)->velocity[0] - v0[0];
-        }
-        if (!NEIGHBOR_MARGIN_j0) {
-          dy += v0[1] - (vert - stride1)->velocity[1];
-        }
-        if (!NEIGHBOR_MARGIN_j1) {
-          dy += (vert + stride1)->velocity[1] - v0[1];
-        }
-        if (!NEIGHBOR_MARGIN_k0) {
-          dz += v0[2] - (vert - stride2)->velocity[2];
-        }
-        if (!NEIGHBOR_MARGIN_k1) {
-          dz += (vert + stride2)->velocity[2] - v0[2];
-        }
-
-        float divergence = -0.5f * flowfac * (dx + dy + dz);
-
-        /* adjustment term for target density */
-        float target = hair_volume_density_divergence(
-            vert->density, target_density, target_strength);
-
-        /* B vector contains the finite difference approximation of the velocity divergence.
-         * Note: according to the discretized Navier-Stokes equation the rhs vector
-         * and resulting pressure gradient should be multiplied by the (inverse) density;
-         * however, this is already included in the weighting of hair velocities on the grid!
-         */
-        B[u] = divergence - target;
-
-#if 0
-        {
-          float wloc[3], loc[3];
-          float col0[3] = {0.0, 0.0, 0.0};
-          float colp[3] = {0.0, 1.0, 1.0};
-          float coln[3] = {1.0, 0.0, 1.0};
-          float col[3];
-          float fac;
-
-          loc[0] = (float)(i - 1);
-          loc[1] = (float)(j - 1);
-          loc[2] = (float)(k - 1);
-          grid_to_world(grid, wloc, loc);
-
-          if (divergence > 0.0f) {
-            fac = CLAMPIS(divergence * target_strength, 0.0, 1.0);
-            interp_v3_v3v3(col, col0, colp, fac);
-          }
-          else {
-            fac = CLAMPIS(-divergence * target_strength, 0.0, 1.0);
-            interp_v3_v3v3(col, col0, coln, fac);
-          }
-          if (fac > 0.05f) {
-            BKE_sim_debug_data_add_circle(
-                grid->debug_data, wloc, 0.01f, col[0], col[1], col[2], "grid", 5522, i, j, k);
-          }
-        }
-#endif
-      }
-    }
-  }
-
-  /* Main Poisson equation system:
-   * This is derived from the discretezation of the Poisson equation
-   *   div(grad(p)) = div(v)
-   *
-   * The finite difference approximation yields the linear equation system described here:
-   * https://en.wikipedia.org/wiki/Discrete_Poisson_equation
-   */
-  lMatrix A(num_cellsA, num_cellsA);
-  /* Reserve space for the base equation system (without boundary conditions).
-   * Each column contains a factor 6 on the diagonal
-   * and up to 6 factors -1 on other places.
-   */
-  A.reserve(Eigen::VectorXi::Constant(num_cellsA, 7));
-
-  for (k = 0; k < resA[2]; k++) {
-    for (j = 0; j < resA[1]; j++) {
-      for (i = 0; i < resA[0]; i++) {
-        int u = i * strideA0 + j * strideA1 + k * strideA2;
-        bool is_margin = MARGIN_i0 || MARGIN_i1 || MARGIN_j0 || MARGIN_j1 || MARGIN_k0 ||
-                         MARGIN_k1;
-
-        vert = vert_start + i * stride0 + j * stride1 + k * stride2;
-        if (!is_margin && vert->density > density_threshold) {
-          int neighbors_lo = 0;
-          int neighbors_hi = 0;
-          int non_solid_neighbors = 0;
-          int neighbor_lo_index[3];
-          int neighbor_hi_index[3];
-          int n;
-
-          /* check for upper bounds in advance
-           * to get the correct number of neighbors,
-           * needed for the diagonal element
-           */
-          if (!NEIGHBOR_MARGIN_k0 && (vert - stride2)->density > density_threshold) {
-            neighbor_lo_index[neighbors_lo++] = u - strideA2;
-          }
-          if (!NEIGHBOR_MARGIN_j0 && (vert - stride1)->density > density_threshold) {
-            neighbor_lo_index[neighbors_lo++] = u - strideA1;
-          }
-          if (!NEIGHBOR_MARGIN_i0 && (vert - stride0)->density > density_threshold) {
-            neighbor_lo_index[neighbors_lo++] = u - strideA0;
-          }
-          if (!NEIGHBOR_MARGIN_i1 && (vert + stride0)->density > density_threshold) {
-            neighbor_hi_index[neighbors_hi++] = u + strideA0;
-          }
-          if (!NEIGHBOR_MARGIN_j1 && (vert + stride1)->density > density_threshold) {
-            neighbor_hi_index[neighbors_hi++] = u + strideA1;
-          }
-          if (!NEIGHBOR_MARGIN_k1 && (vert + stride2)->density > density_threshold) {
-            neighbor_hi_index[neighbors_hi++] = u + strideA2;
-          }
-
-          /*int liquid_neighbors = neighbors_lo + neighbors_hi;*/
-          non_solid_neighbors = 6;
-
-          for (n = 0; n < neighbors_lo; n++) {
-            A.insert(neighbor_lo_index[n], u) = -1.0f;
-          }
-          A.insert(u, u) = (float)non_solid_neighbors;
-          for (n = 0; n < neighbors_hi; n++) {
-            A.insert(neighbor_hi_index[n], u) = -1.0f;
-          }
-        }
-        else {
-          A.insert(u, u) = 1.0f;
-        }
-      }
-    }
-  }
-
-  ConjugateGradient cg;
-  cg.setMaxIterations(100);
-  cg.setTolerance(0.01f);
-
-  cg.compute(A);
-
-  lVector p = cg.solve(B);
-
-  if (cg.info() == Eigen::Success) {
-    /* Calculate velocity = grad(p) */
-    for (k = 0; k < resA[2]; k++) {
-      for (j = 0; j < resA[1]; j++) {
-        for (i = 0; i < resA[0]; i++) {
-          int u = i * strideA0 + j * strideA1 + k * strideA2;
-          bool is_margin = MARGIN_i0 || MARGIN_i1 || MARGIN_j0 || MARGIN_j1 || MARGIN_k0 ||
-                           MARGIN_k1;
-          if (is_margin) {
-            continue;
-          }
-
-          vert = vert_start + i * stride0 + j * stride1 + k * stride2;
-          if (vert->density > density_threshold) {
-            float p_left = p[u - strideA0];
-            float p_right = p[u + strideA0];
-            float p_down = p[u - strideA1];
-            float p_up = p[u + strideA1];
-            float p_bottom = p[u - strideA2];
-            float p_top = p[u + strideA2];
-
-            /* finite difference estimate of pressure gradient */
-            float dvel[3];
-            dvel[0] = p_right - p_left;
-            dvel[1] = p_up - p_down;
-            dvel[2] = p_top - p_bottom;
-            mul_v3_fl(dvel, -0.5f * inv_flowfac);
-
-            /* pressure gradient describes velocity delta */
-            add_v3_v3v3(vert->velocity_smooth, vert->velocity, dvel);
-          }
-          else {
-            zero_v3(vert->velocity_smooth);
-          }
-        }
-      }
-    }
-
-#if 0
-    {
-      int axis = 0;
-      float offset = 0.0f;
-
-      int slice = (offset - grid->gmin[axis]) / grid->cellsize;
-
-      for (k = 0; k < resA[2]; k++) {
-        for (j = 0; j < resA[1]; j++) {
-          for (i = 0; i < resA[0]; i++) {
-            int u = i * strideA0 + j * strideA1 + k * strideA2;
-            bool is_margin = MARGIN_i0 || MARGIN_i1 || MARGIN_j0 || MARGIN_j1 || MARGIN_k0 ||
-                             MARGIN_k1;
-            if (i != slice) {
-              continue;
-            }
-
-            vert = vert_start + i * stride0 + j * stride1 + k * stride2;
-
-            float wloc[3], loc[3];
-            float col0[3] = {0.0, 0.0, 0.0};
-            float colp[3] = {0.0, 1.0, 1.0};
-            float coln[3] = {1.0, 0.0, 1.0};
-            float col[3];
-            float fac;
-
-            loc[0] = (float)(i - 1);
-            loc[1] = (float)(j - 1);
-            loc[2] = (float)(k - 1);
-            grid_to_world(grid, wloc, loc);
-
-            float pressure = p[u];
-            if (pressure > 0.0f) {
-              fac = CLAMPIS(pressure * grid->debug1, 0.0, 1.0);
-              interp_v3_v3v3(col, col0, colp, fac);
-            }
-            else {
-              fac = CLAMPIS(-pressure * grid->debug1, 0.0, 1.0);
-              interp_v3_v3v3(col, col0, coln, fac);
-            }
-            if (fac > 0.05f) {
-              BKE_sim_debug_data_add_circle(
-                  grid->debug_data, wloc, 0.01f, col[0], col[1], col[2], "grid", 5533, i, j, k);
-            }
-
-            if (!is_margin) {
-              float dvel[3];
-              sub_v3_v3v3(dvel, vert->velocity_smooth, vert->velocity);
-              // BKE_sim_debug_data_add_vector(
-              //     grid->debug_data, wloc, dvel, 1, 1, 1,
-              //     "grid", 5566, i, j, k);
-            }
-
-            if (!is_margin) {
-              float d = CLAMPIS(vert->density * grid->debug2, 0.0f, 1.0f);
-              float col0[3] = {0.3, 0.3, 0.3};
-              float colp[3] = {0.0, 0.0, 1.0};
-              float col[3];
-
-              interp_v3_v3v3(col, col0, colp, d);
-              // if (d > 0.05f) {
-              // BKE_sim_debug_data_add_dot(
-              //     grid->debug_data, wloc, col[0], col[1], col[2],
-              //     "grid", 5544, i, j, k);
-              // }
-            }
-          }
-        }
-      }
-    }
-#endif
-
-    return true;
-  }
-  else {
-    /* Clear result in case of error */
-    for (i = 0, vert = grid->verts; i < num_cells; i++, vert++) {
-      zero_v3(vert->velocity_smooth);
-    }
-
-    return false;
-  }
-}
-
-#if 0 /* XXX weighting is incorrect, disabled for now */
-/* Velocity filter kernel
- * See https://en.wikipedia.org/wiki/Filter_%28large_eddy_simulation%29
- */
-
-BLI_INLINE void hair_volume_filter_box_convolute(
-    HairVertexGrid *grid, float invD, const int kernel_size[3], int i, int j, int k)
-{
-  int res = grid->res;
-  int p, q, r;
-  int minp = max_ii(i - kernel_size[0], 0), maxp = min_ii(i + kernel_size[0], res - 1);
-  int minq = max_ii(j - kernel_size[1], 0), maxq = min_ii(j + kernel_size[1], res - 1);
-  int minr = max_ii(k - kernel_size[2], 0), maxr = min_ii(k + kernel_size[2], res - 1);
-  int offset, kernel_offset, kernel_dq, kernel_dr;
-  HairGridVert *verts;
-  float *vel_smooth;
-
-  offset = i + (j + k * res) * res;
-  verts = grid->verts;
-  vel_smooth = verts[offset].velocity_smooth;
-
-  kernel_offset = minp + (minq + minr * res) * res;
-  kernel_dq = res;
-  kernel_dr = res * res;
-  for (r = minr; r <= maxr; r++) {
-    for (q = minq; q <= maxq; q++) {
-      for (p = minp; p <= maxp; p++) {
-
-        madd_v3_v3fl(vel_smooth, verts[kernel_offset].velocity, invD);
-
-        kernel_offset += 1;
-      }
-      kernel_offset += kernel_dq;
-    }
-    kernel_offset += kernel_dr;
-  }
-}
-
-void BPH_hair_volume_vertex_grid_filter_box(HairVertexGrid *grid, int kernel_size)
-{
-  int size = hair_grid_size(grid->res);
-  int kernel_sizev[3] = {kernel_size, kernel_size, kernel_size};
-  int tot;
-  float invD;
-  int i, j, k;
-
-  if (kernel_size <= 0) {
-    return;
-  }
-
-  tot = kernel_size * 2 + 1;
-  invD = 1.0f / (float)(tot * tot * tot);
-
-  /* clear values for convolution */
-  for (i = 0; i < size; i++) {
-    zero_v3(grid->verts[i].velocity_smooth);
-  }
-
-  for (i = 0; i < grid->res; i++) {
-    for (j = 0; j < grid->res; j++) {
-      for (k = 0; k < grid->res; k++) {
-        hair_volume_filter_box_convolute(grid, invD, kernel_sizev, i, j, k);
-      }
-    }
-  }
-
-  /* apply as new velocity */
-  for (i = 0; i < size; i++) {
-    copy_v3_v3(grid->verts[i].velocity, grid->verts[i].velocity_smooth);
-  }
-}
-#endif
-
-HairGrid *BPH_hair_volume_create_vertex_grid(float cellsize,
-                                             const float gmin[3],
-                                             const float gmax[3])
-{
-  float scale;
-  float extent[3];
-  int resmin[3], resmax[3], res[3];
-  float gmin_margin[3], gmax_margin[3];
-  int size;
-  HairGrid *grid;
-  int i;
-
-  /* sanity check */
-  if (cellsize <= 0.0f) {
-    cellsize = 1.0f;
-  }
-  scale = 1.0f / cellsize;
-
-  sub_v3_v3v3(extent, gmax, gmin);
-  for (i = 0; i < 3; i++) {
-    resmin[i] = floor_int(gmin[i] * scale);
-    resmax[i] = floor_int(gmax[i] * scale) + 1;
-
-    /* add margin of 1 cell */
-    resmin[i] -= 1;
-    resmax[i] += 1;
-
-    res[i] = resmax[i] - resmin[i] + 1;
-    /* sanity check: avoid null-sized grid */
-    if (res[i] < 4) {
-      res[i] = 4;
-      resmax[i] = resmin[i] + 4;
-    }
-    /* sanity check: avoid too large grid size */
-    if (res[i] > MAX_HAIR_GRID_RES) {
-      res[i] = MAX_HAIR_GRID_RES;
-      resmax[i] = resmin[i] + MAX_HAIR_GRID_RES;
-    }
-
-    gmin_margin[i] = (float)resmin[i] * cellsize;
-    gmax_margin[i] = (float)resmax[i] * cellsize;
-  }
-  size = hair_grid_size(res);
-
-  grid = (HairGrid *)MEM_callocN(sizeof(HairGrid), "hair grid");
-  grid->res[0] = res[0];
-  grid->res[1] = res[1];
-  grid->res[2] = res[2];
-  copy_v3_v3(grid->gmin, gmin_margin);
-  copy_v3_v3(grid->gmax, gmax_margin);
-  grid->cellsize = cellsize;
-  grid->inv_cellsize = scale;
-  grid->verts = (HairGridVert *)MEM_callocN(sizeof(HairGridVert) * size, "hair voxel data");
-
-  return grid;
-}
-
-void BPH_hair_volume_free_vertex_grid(HairGrid *grid)
-{
-  if (grid) {
-    if (grid->verts) {
-      MEM_freeN(grid->verts);
-    }
-    MEM_freeN(grid);
-  }
-}
-
-void BPH_hair_volume_grid_geometry(
-    HairGrid *grid, float *cellsize, int res[3], float gmin[3], float gmax[3])
-{
-  if (cellsize) {
-    *cellsize = grid->cellsize;
-  }
-  if (res) {
-    copy_v3_v3_int(res, grid->res);
-  }
-  if (gmin) {
-    copy_v3_v3(gmin, grid->gmin);
-  }
-  if (gmax) {
-    copy_v3_v3(gmax, grid->gmax);
-  }
-}
-
-#if 0
-static HairGridVert *hair_volume_create_collision_grid(ClothModifierData *clmd,
-                                                       lfVector *lX,
-                                                       unsigned int numverts)
-{
-  int res = hair_grid_res;
-  int size = hair_grid_size(res);
-  HairGridVert *collgrid;
-  ListBase *colliders;
-  ColliderCache *col = NULL;
-  float gmin[3], gmax[3], scale[3];
-  /* 2.0f is an experimental value that seems to give good results */
-  float collfac = 2.0f * clmd->sim_parms->collider_friction;
-  unsigned int v = 0;
-  int i = 0;
-
-  hair_volume_get_boundbox(lX, numverts, gmin, gmax);
-  hair_grid_get_scale(res, gmin, gmax, scale);
-
-  collgrid = MEM_mallocN(sizeof(HairGridVert) * size, "hair collider voxel data");
-
-  /* initialize grid */
-  for (i = 0; i < size; i++) {
-    zero_v3(collgrid[i].velocity);
-    collgrid[i].density = 0.0f;
-  }
-
-  /* gather colliders */
-  colliders = BKE_collider_cache_create(depsgraph, NULL, NULL);
-  if (colliders && collfac > 0.0f) {
-    for (col = colliders->first; col; col = col->next) {
-      MVert *loc0 = col->collmd->x;
-      MVert *loc1 = col->collmd->xnew;
-      float vel[3];
-      float weights[8];
-      int di, dj, dk;
-
-      for (v = 0; v < col->collmd->numverts; v++, loc0++, loc1++) {
-        int offset;
-
-        if (!hair_grid_point_valid(loc1->co, gmin, gmax)) {
-          continue;
-        }
-
-        offset = hair_grid_weights(res, gmin, scale, lX[v], weights);
-
-        sub_v3_v3v3(vel, loc1->co, loc0->co);
-
-        for (di = 0; di < 2; di++) {
-          for (dj = 0; dj < 2; dj++) {
-            for (dk = 0; dk < 2; dk++) {
-              int voffset = offset + di + (dj + dk * res) * res;
-              int iw = di + dj * 2 + dk * 4;
-
-              collgrid[voffset].density += weights[iw];
-              madd_v3_v3fl(collgrid[voffset].velocity, vel, weights[iw]);
-            }
-          }
-        }
-      }
-    }
-  }
-  BKE_collider_cache_free(&colliders);
-
-  /* divide velocity with density */
-  for (i = 0; i < size; i++) {
-    float density = collgrid[i].density;
-    if (density > 0.0f) {
-      mul_v3_fl(collgrid[i].velocity, 1.0f / density);
-    }
-  }
-
-  return collgrid;
-}
-#endif