diff options
-rw-r--r-- | source/blender/physics/intern/BPH_mass_spring.cpp | 15 | ||||
-rw-r--r-- | source/blender/physics/intern/eigen_utils.h | 34 | ||||
-rw-r--r-- | source/blender/physics/intern/hair_volume.cpp | 220 | ||||
-rw-r--r-- | source/blender/physics/intern/implicit.h | 6 |
4 files changed, 239 insertions, 36 deletions
diff --git a/source/blender/physics/intern/BPH_mass_spring.cpp b/source/blender/physics/intern/BPH_mass_spring.cpp index dfc08b5574b..fbb1cb7cf2c 100644 --- a/source/blender/physics/intern/BPH_mass_spring.cpp +++ b/source/blender/physics/intern/BPH_mass_spring.cpp @@ -745,7 +745,7 @@ static void cloth_continuum_fill_grid(HairGrid *grid, Cloth *cloth) BPH_hair_volume_normalize_vertex_grid(grid);
}
-static void cloth_continuum_step(ClothModifierData *clmd)
+static void cloth_continuum_step(ClothModifierData *clmd, float dt)
{
ClothSimSettings *parms = clmd->sim_parms;
Cloth *cloth = clmd->clothObject;
@@ -774,6 +774,9 @@ static void cloth_continuum_step(ClothModifierData *clmd) cloth_continuum_fill_grid(grid, cloth);
+ /* main hair continuum solver */
+ BPH_hair_volume_solve_divergence(grid, dt);
+
for (i = 0, vert = cloth->verts; i < numverts; i++, vert++) {
float x[3], v[3], nv[3];
@@ -811,17 +814,19 @@ static void cloth_continuum_step(ClothModifierData *clmd) BKE_sim_debug_data_clear_category(clmd->debug_data, "grid velocity");
for (j = 0; j < size; ++j) {
for (i = 0; i < size; ++i) {
- float x[3], v[3], gvel[3], gdensity;
+ float x[3], v[3], gvel[3], gvel_smooth[3], gdensity;
madd_v3_v3v3fl(x, offset, a, (float)i / (float)(size-1));
madd_v3_v3fl(x, b, (float)j / (float)(size-1));
zero_v3(v);
- BPH_hair_volume_grid_interpolate(grid, x, &gdensity, gvel, NULL, NULL);
+ BPH_hair_volume_grid_interpolate(grid, x, &gdensity, gvel, gvel_smooth, NULL, NULL);
// BKE_sim_debug_data_add_circle(clmd->debug_data, x, gdensity, 0.7, 0.3, 1, "grid density", hash_int_2d(hash_int_2d(i, j), 3111));
- if (!is_zero_v3(gvel))
+ if (!is_zero_v3(gvel) || !is_zero_v3(gvel_smooth)) {
BKE_sim_debug_data_add_vector(clmd->debug_data, x, gvel, 0.4, 0, 1, "grid velocity", hash_int_2d(hash_int_2d(i, j), 3112));
+ BKE_sim_debug_data_add_vector(clmd->debug_data, x, gvel_smooth, 0.6, 4, 1, "grid velocity", hash_int_2d(hash_int_2d(i, j), 3113));
+ }
}
}
}
@@ -950,7 +955,7 @@ int BPH_cloth_solve(Object *ob, float frame, ClothModifierData *clmd, ListBase * BPH_mass_spring_solve_velocities(id, dt, &result);
cloth_record_result(clmd, &result, clmd->sim_parms->stepsPerFrame);
- cloth_continuum_step(clmd);
+ cloth_continuum_step(clmd, dt);
BPH_mass_spring_solve_positions(id, dt);
diff --git a/source/blender/physics/intern/eigen_utils.h b/source/blender/physics/intern/eigen_utils.h index 2112386781f..cd3bd5ee45d 100644 --- a/source/blender/physics/intern/eigen_utils.h +++ b/source/blender/physics/intern/eigen_utils.h @@ -53,21 +53,21 @@ typedef float Scalar; /* slightly extended Eigen vector class * with conversion to/from plain C float array */ -class fVector : public Eigen::Vector3f { +class Vector3 : public Eigen::Vector3f { public: typedef float *ctype; - fVector() + Vector3() { } - fVector(const ctype &v) + Vector3(const ctype &v) { for (int k = 0; k < 3; ++k) coeffRef(k) = v[k]; } - fVector& operator = (const ctype &v) + Vector3& operator = (const ctype &v) { for (int k = 0; k < 3; ++k) coeffRef(k) = v[k]; @@ -83,22 +83,22 @@ public: /* slightly extended Eigen matrix class * with conversion to/from plain C float array */ -class fMatrix : public Eigen::Matrix3f { +class Matrix3 : public Eigen::Matrix3f { public: typedef float (*ctype)[3]; - fMatrix() + Matrix3() { } - fMatrix(const ctype &v) + Matrix3(const ctype &v) { for (int k = 0; k < 3; ++k) for (int l = 0; l < 3; ++l) coeffRef(l, k) = v[k][l]; } - fMatrix& operator = (const ctype &v) + Matrix3& operator = (const ctype &v) { for (int k = 0; k < 3; ++k) for (int l = 0; l < 3; ++l) @@ -112,19 +112,21 @@ public: } }; +typedef Eigen::VectorXf lVector; + /* Extension of dense Eigen vectors, * providing 3-float block access for blenlib math functions */ -class lVector : public Eigen::VectorXf { +class lVector3f : public Eigen::VectorXf { public: typedef Eigen::VectorXf base_t; - lVector() + lVector3f() { } template <typename T> - lVector& operator = (T rhs) + lVector3f& operator = (T rhs) { base_t::operator=(rhs); return *this; @@ -151,8 +153,8 @@ typedef Eigen::SparseMatrix<Scalar> lMatrix; * This should be used for building lMatrix instead of writing to such lMatrix directly (which is very inefficient). * After all elements have been defined using the set() method, the actual matrix can be filled using construct(). */ -struct lMatrixCtor { - lMatrixCtor() +struct lMatrix3fCtor { + lMatrix3fCtor() { } @@ -167,7 +169,7 @@ struct lMatrixCtor { m_trips.reserve(numverts * 9); } - void add(int i, int j, const fMatrix &m) + void add(int i, int j, const Matrix3 &m) { i *= 3; j *= 3; @@ -176,7 +178,7 @@ struct lMatrixCtor { m_trips.push_back(Triplet(i + k, j + l, m.coeff(l, k))); } - void sub(int i, int j, const fMatrix &m) + void sub(int i, int j, const Matrix3 &m) { i *= 3; j *= 3; @@ -199,7 +201,7 @@ typedef Eigen::ConjugateGradient<lMatrix, Eigen::Lower, Eigen::DiagonalPrecondit using Eigen::ComputationInfo; -BLI_INLINE void print_lvector(const lVector &v) +BLI_INLINE void print_lvector(const lVector3f &v) { for (int i = 0; i < v.rows(); ++i) { if (i > 0 && i % 3 == 0) diff --git a/source/blender/physics/intern/hair_volume.cpp b/source/blender/physics/intern/hair_volume.cpp index cda56946345..cd4e6503965 100644 --- a/source/blender/physics/intern/hair_volume.cpp +++ b/source/blender/physics/intern/hair_volume.cpp @@ -39,6 +39,7 @@ #include "BKE_effect.h" #include "implicit.h" +#include "eigen_utils.h" /* ================ Volumetric Hair Interaction ================ * adapted from @@ -113,7 +114,7 @@ BLI_INLINE int hair_grid_interp_weights(const int res[3], const float gmin[3], f } 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 density_gradient[3], float velocity_gradient[3][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]; @@ -151,6 +152,16 @@ BLI_INLINE void hair_grid_interpolate(const HairGridVert *grid, const int res[3] } } + 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 ) + @@ -180,7 +191,7 @@ void BPH_hair_volume_vertex_grid_forces(HairGrid *grid, const float x[3], const { 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, ggrad, gvelgrad); + 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); @@ -199,21 +210,32 @@ void BPH_hair_volume_vertex_grid_forces(HairGrid *grid, const float x[3], const } void BPH_hair_volume_grid_interpolate(HairGrid *grid, const float x[3], - float *density, float velocity[3], float density_gradient[3], float velocity_gradient[3][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, density_gradient, velocity_gradient); + 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], ggrad[3], gvelgrad[3][3]; + float gdensity, gvelocity[3], gvel_smooth[3], ggrad[3], gvelgrad[3][3]; - hair_grid_interpolate(grid->verts, grid->res, grid->gmin, grid->inv_cellsize, x, &gdensity, gvelocity, ggrad, gvelgrad); + hair_grid_interpolate(grid->verts, grid->res, grid->gmin, grid->inv_cellsize, x, &gdensity, gvelocity, gvel_smooth, ggrad, gvelgrad); /* XXX TODO implement FLIP method and use fluid_factor to blend between FLIP and PIC */ - copy_v3_v3(r_v, gvelocity); + copy_v3_v3(r_v, gvel_smooth); +} + +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; + } } BLI_INLINE bool hair_grid_point_valid(const float vec[3], float gmin[3], float gmax[3]) @@ -464,6 +486,182 @@ void BPH_hair_volume_normalize_vertex_grid(HairGrid *grid) } } +bool BPH_hair_volume_solve_divergence(HairGrid *grid, float dt) +{ + const float density_threshold = 0.001f; /* cells with density below this are considered empty */ + + const float flowfac = grid->cellsize / dt; + const float inv_flowfac = dt / grid->cellsize; + + /*const int num_cells = hair_grid_size(grid->res);*/ + const int inner_res[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); + + /* NB: to avoid many boundary checks, we only solve the system + * for the inner vertices, excluding a 1-cell margin. + */ + const int inner_cells_start = stride0 + stride1 + stride2; + const int num_inner_cells = inner_res[0] * inner_res[1] * inner_res[2]; + + HairGridVert *vert; + int i, j, k, u; + + BLI_assert(num_inner_cells >= 1); + + /* Calculate divergence */ + lVector B(num_inner_cells); + for (k = 0; k < inner_res[2]; ++k) { + for (j = 0; j < inner_res[1]; ++j) { + for (i = 0; i < inner_res[0]; ++i) { + u = i * strideA0 + j * strideA1 + k * strideA2; + vert = grid->verts + inner_cells_start + i * stride0 + j * stride1 + k * stride2; + + HairGridVert *vert_px = vert + stride0; + HairGridVert *vert_py = vert + stride1; + HairGridVert *vert_pz = vert + stride2; + + const float *v = vert->velocity; + float dx = vert_px->velocity[0] - v[0]; + float dy = vert_py->velocity[1] - v[1]; + float dz = vert_pz->velocity[2] - v[2]; + + /* 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] = (dx + dy + dz) * flowfac; + } + } + } + + /* 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: + * http://en.wikipedia.org/wiki/Discrete_Poisson_equation + */ + lMatrix A(num_inner_cells, num_inner_cells); + /* 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_inner_cells, 7)); + + for (k = 0; k < inner_res[2]; ++k) { + for (j = 0; j < inner_res[1]; ++j) { + for (i = 0; i < inner_res[0]; ++i) { + u = i * strideA0 + j * strideA1 + k * strideA2; + vert = grid->verts + inner_cells_start + i * stride0 + j * stride1 + k * stride2; + + if (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 (k >= 1) { + if ((vert - stride2)->density > density_threshold) + neighbor_lo_index[neighbors_lo++] = u - strideA2; + } + if (j >= 1) { + if ((vert - stride1)->density > density_threshold) + neighbor_lo_index[neighbors_lo++] = u - strideA1; + } + if (i >= 1) { + if ((vert - stride0)->density > density_threshold) + neighbor_lo_index[neighbors_lo++] = u - strideA0; + } + if (i < inner_res[0] - 1) { + if ((vert + stride0)->density > density_threshold) + neighbor_hi_index[neighbors_hi++] = u + strideA0; + } + if (j < inner_res[1] - 1) { + if ((vert + stride1)->density > density_threshold) + neighbor_hi_index[neighbors_hi++] = u + strideA1; + } + if (k < inner_res[2] - 1) { + if ((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 < inner_res[2]; ++k) { + for (j = 0; j < inner_res[1]; ++j) { + for (i = 0; i < inner_res[0]; ++i) { + u = i * strideA0 + j * strideA1 + k * strideA2; + vert = grid->verts + inner_cells_start + i * stride0 + j * stride1 + k * stride2; + + if (vert->density > density_threshold) { + float p0 = p[u]; + + /* finite difference estimate of pressure gradient */ + float grad_p[3]; + grad_p[0] = i >= 1 ? p0 - p[u - strideA0] : 0.0f; + grad_p[1] = j >= 1 ? p0 - p[u - strideA1] : 0.0f; + grad_p[2] = k >= 1 ? p0 - p[u - strideA2] : 0.0f; + + /* pressure gradient describes velocity delta */ + madd_v3_v3v3fl(vert->velocity_smooth, vert->velocity, grad_p, inv_flowfac); + } + else { + zero_v3(vert->velocity_smooth); + } + } + } + } + + return true; + } + else { + /* Clear result in case of error */ + for (i = inner_cells_start, vert = grid->verts + inner_cells_start; i < num_inner_cells; ++i, ++vert) { + zero_v3(vert->velocity_smooth); + } + + return false; + } +} + #if 0 /* XXX weighting is incorrect, disabled for now */ /* Velocity filter kernel * See http://en.wikipedia.org/wiki/Filter_%28large_eddy_simulation%29 @@ -584,14 +782,8 @@ HairGrid *BPH_hair_volume_create_vertex_grid(float cellsize, const float gmin[3] copy_v3_v3(grid->gmax, gmax_margin); grid->cellsize = cellsize; grid->inv_cellsize = scale; - grid->verts = (HairGridVert *)MEM_mallocN(sizeof(HairGridVert) * size, "hair voxel data"); + grid->verts = (HairGridVert *)MEM_callocN(sizeof(HairGridVert) * size, "hair voxel data"); - /* initialize grid */ - for (i = 0; i < size; ++i) { - zero_v3(grid->verts[i].velocity); - grid->verts[i].density = 0.0f; - } - return grid; } diff --git a/source/blender/physics/intern/implicit.h b/source/blender/physics/intern/implicit.h index 373a73cf53f..58883a26c85 100644 --- a/source/blender/physics/intern/implicit.h +++ b/source/blender/physics/intern/implicit.h @@ -177,6 +177,7 @@ void BPH_hair_volume_free_vertex_grid(struct HairGrid *grid); void BPH_hair_volume_set_debug_data(struct HairGrid *grid, struct SimDebugData *debug_data); void BPH_hair_volume_grid_geometry(struct HairGrid *grid, float *cellsize, int res[3], float gmin[3], float gmax[3]); +void BPH_hair_volume_grid_clear(struct HairGrid *grid); void BPH_hair_volume_add_vertex(struct HairGrid *grid, const float x[3], const float v[3]); void BPH_hair_volume_add_segment(struct HairGrid *grid, const float x1[3], const float v1[3], const float x2[3], const float v2[3], @@ -184,12 +185,15 @@ void BPH_hair_volume_add_segment(struct HairGrid *grid, const float dir1[3], const float dir2[3], const float dir3[3]); void BPH_hair_volume_normalize_vertex_grid(struct HairGrid *grid); + +bool BPH_hair_volume_solve_divergence(struct HairGrid *grid, float dt); #if 0 /* XXX weighting is incorrect, disabled for now */ void BPH_hair_volume_vertex_grid_filter_box(struct HairVertexGrid *grid, int kernel_size); #endif void BPH_hair_volume_grid_interpolate(struct HairGrid *grid, const float x[3], - float *density, float velocity[3], float density_gradient[3], float velocity_gradient[3][3]); + float *density, float velocity[3], float velocity_smooth[3], + float density_gradient[3], float velocity_gradient[3][3]); /* Effect of fluid simulation grid on velocities. * fluid_factor controls blending between PIC (Particle-in-Cell) |