[Bf-blender-cvs] [926a674f] master: Main solver step for generating a divergence-free hair velocity field on the grid.

Lukas Tönne noreply at git.blender.org
Tue Jan 20 09:52:34 CET 2015


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. @@




More information about the Bf-blender-cvs mailing list