[Bf-blender-cvs] [90f9fa5] temp_hair_flow: Solver for the pressure Poisson equation on hair flow grids.

Lukas Tönne noreply at git.blender.org
Thu Jan 8 10:22:53 CET 2015


Commit: 90f9fa579a9bbbcdca9a5d1dea253a7e5a077130
Author: Lukas Tönne
Date:   Thu Jan 8 10:22:05 2015 +0100
Branches: temp_hair_flow
https://developer.blender.org/rB90f9fa579a9bbbcdca9a5d1dea253a7e5a077130

Solver for the pressure Poisson equation on hair flow grids.

Note: does not include boundary conditions so far.

===================================================================

M	source/blender/physics/intern/grid.cpp
M	source/blender/physics/intern/grid.h
M	source/blender/physics/intern/hair_flow.cpp

===================================================================

diff --git a/source/blender/physics/intern/grid.cpp b/source/blender/physics/intern/grid.cpp
index 51fddb8..ff188f0 100644
--- a/source/blender/physics/intern/grid.cpp
+++ b/source/blender/physics/intern/grid.cpp
@@ -224,6 +224,90 @@ void Grid::calc_divergence(GridHash<float> &divergence, const GridHash<bool> &so
 	}
 }
 
+/* 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
+ * 
+ * For a good overview of eulerian fluid sim methods, see
+ * http://www.proxyarch.com/util/techpapers/papers/Fluid%20flow%20for%20the%20rest%20of%20us.pdf
+ */
+void Grid::solve_pressure(GridHash<float> &pressure, const GridHash<float> &divergence)
+{
+	int stride[3] = { 1, res[0], res[0] * res[1] };
+	
+	lMatrix A(num_cells, num_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_cells, 7));
+	
+	for (int x = 0; x < res[0]; ++x) {
+		for (int y = 0; y < res[1]; ++y) {
+			for (int z = 0; z < res[2]; ++z) {
+				int u = x * stride[0] + y * stride[1] + z * stride[2];
+				bool is_margin = !(x > 0 && x < res[0]-1 && y > 0 && y < res[1]-1 && z > 0 && z < res[2]-1);
+				if (is_margin) {
+					A.insert(u, u) = 1.0f;
+					continue;
+				}
+				
+				/* check for upper bounds in advance
+				 * to get the correct number of neighbors,
+				 * needed for the diagonal element
+				 */
+				int neighbors_lo = 0;
+				int neighbors_hi = 0;
+				int non_solid_neighbors = 0;
+				int neighbor_lo_index[3];
+				int neighbor_hi_index[3];
+				if (z > 1)
+					neighbor_lo_index[neighbors_lo++] = u - stride[2];
+				if (y > 1)
+					neighbor_lo_index[neighbors_lo++] = u - stride[1];
+				if (x > 1)
+					neighbor_lo_index[neighbors_lo++] = u - stride[0];
+				if (x < res[0]-2)
+					neighbor_hi_index[neighbors_hi++] = u + stride[0];
+				if (y < res[1]-2)
+					neighbor_hi_index[neighbors_hi++] = u + stride[1];
+				if (z < res[2]-2)
+					neighbor_hi_index[neighbors_hi++] = u + stride[2];
+				
+				/*int liquid_neighbors = neighbors_lo + neighbors_hi;*/
+				non_solid_neighbors = 6;
+				
+				for (int n = 0; n < neighbors_lo; ++n)
+					A.insert(neighbor_lo_index[n], u) = -1.0f;
+				A.insert(u, u) = (float)non_solid_neighbors;
+				for (int n = 0; n < neighbors_hi; ++n)
+					A.insert(neighbor_hi_index[n], u) = -1.0f;
+			}
+		}
+	}
+	
+	ConjugateGradient cg;
+	cg.setMaxIterations(100);
+	cg.setTolerance(0.01f);
+	
+	cg.compute(A);
+	
+	lVector B(num_cells);
+	divergence.to_eigen(B);
+	
+	lVector p = cg.solve(B);
+	
+	if (cg.info() == Eigen::Success) {
+		pressure.from_eigen(p);
+	}
+	else
+		pressure.clear();
+}
+
 } /* namespace BPH */
 
 #include "implicit.h"
diff --git a/source/blender/physics/intern/grid.h b/source/blender/physics/intern/grid.h
index fce336b..206225f 100644
--- a/source/blender/physics/intern/grid.h
+++ b/source/blender/physics/intern/grid.h
@@ -57,6 +57,9 @@ struct GridHash {
 			MEM_freeN(m_data);
 	}
 	
+	const int* resolution() const { return m_res; }
+	int size() const { return m_size; }
+	
 	void resize(const int res[3])
 	{
 		if (m_data)
@@ -104,6 +107,22 @@ struct GridHash {
 		return *(m_data + x + y * m_stride[1] + z * m_stride[2]);
 	}
 	
+	void from_eigen(lVector &r) const
+	{
+		BLI_assert(r.rows() == m_size);
+		for (int i = 0; i < m_size; ++i) {
+			m_data[i] = r.coeff(i);
+		}
+	}
+	
+	void to_eigen(lVector &r) const
+	{
+		BLI_assert(r.rows() == m_size);
+		for (int i = 0; i < m_size; ++i) {
+			r.coeffRef(i) = m_data[i];
+		}
+	}
+	
 private:
 	T *m_data; /* XXX TODO stub array for now, actually use a hash table here! */
 	int m_res[3];
diff --git a/source/blender/physics/intern/hair_flow.cpp b/source/blender/physics/intern/hair_flow.cpp
index 06833c0..e0c0696 100644
--- a/source/blender/physics/intern/hair_flow.cpp
+++ b/source/blender/physics/intern/hair_flow.cpp
@@ -144,11 +144,15 @@ HairFlowData *BPH_strands_solve_hair_flow(Scene *scene, Object *ob, float max_le
 	divergence.resize(data->grid.res);
 	data->grid.calc_divergence(divergence, source, source_normal);
 	
+	GridHash<float> pressure;
+	pressure.resize(data->grid.res);
+	data->grid.solve_pressure(pressure, divergence);
+	
 	{
 		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 divfac = 10.0f;
+		float colfac = 10.0f;
 		
 		BKE_sim_debug_data_clear_category(debug_data, "hair flow");
 		
@@ -161,20 +165,22 @@ HairFlowData *BPH_strands_solve_hair_flow(Scene *scene, Object *ob, float max_le
 					
 					bool is_source = *source.get(x, y, z);
 					float div = *divergence.get(x, y, z);
+					float prs = *pressure.get(x, y, z);
 					
 //					if (is_source)
 //						BKE_sim_debug_data_add_circle(debug_data, vec, 0.02f, 1,0,0, "hair flow", 111, x, y, z);
 //					else
 //						BKE_sim_debug_data_add_circle(debug_data, vec, 0.02f, 0,1,0, "hair flow", 111, x, y, z);
 					
+					float input = prs;
 					float fac;
 					float col[3];
 					if (div > 0.0f) {
-						fac = CLAMPIS(div * divfac, 0.0, 1.0);
+						fac = CLAMPIS(input * colfac, 0.0, 1.0);
 						interp_v3_v3v3(col, col0, colp, fac);
 					}
 					else {
-						fac = CLAMPIS(-div * divfac, 0.0, 1.0);
+						fac = CLAMPIS(-input * colfac, 0.0, 1.0);
 						interp_v3_v3v3(col, col0, coln, fac);
 					}
 					if (fac > 0.05f)
@@ -182,32 +188,6 @@ HairFlowData *BPH_strands_solve_hair_flow(Scene *scene, Object *ob, float max_le
 				}
 			}
 		}
-#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
 	}
 	
 	return data;




More information about the Bf-blender-cvs mailing list