Commit: 926a674fe804a2ac802650c9fe7b82618addc9a7
Author: Lukas Tönne
Date: Sat Nov 8 18:45:28 2014 +0100
Branches: master
https://developer.blender.org/rB926a674fe804a2ac802650c9fe7b82618addc9a7
Main solver step for generating a divergence-free hair velocity field
on the grid.
This uses the Eigen conjugate-gradient solver to solve the implicit
Poisson equation for the pressure Laplacian:
div(grad(p)) = div(v)
As described in "Detail Preserving Continuum Simulation of Straight Hair"
(McAdams, Selle, 2009).
Conflicts:
source/blender/physics/intern/BPH_mass_spring.cpp
===================================================================
M source/blender/physics/intern/BPH_mass_spring.cpp
M source/blender/physics/intern/eigen_utils.h
M source/blender/physics/intern/hair_volume.cpp
M source/blender/physics/intern/implicit.h
===================================================================
diff --git a/source/blender/physics/intern/BPH_mass_spring.cpp
b/source/blender/physics/intern/BPH_mass_spring.cpp
index dfc08b5..fbb1cb7 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 2112386..cd3bd5e 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 cda5694..cd4e650 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;
@@ Diff output truncated at 10240 characters. @@
_______________________________________________
Bf-blender-cvs mailing list
[email protected]
http://lists.blender.org/mailman/listinfo/bf-blender-cvs
