[Bf-blender-cvs] [e73df24] master: Added a margin to the number of cells used in the poisson grid solver, to ensure we always have one layer of empty cells around the fluid.

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


Commit: e73df249c7f0e5acfcb187e6eb86246b44ae06f9
Author: Lukas Tönne
Date:   Fri Nov 14 10:42:09 2014 +0100
Branches: master
https://developer.blender.org/rBe73df249c7f0e5acfcb187e6eb86246b44ae06f9

Added a margin to the number of cells used in the poisson grid solver,
to ensure we always have one layer of empty cells around the fluid.

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

M	source/blender/physics/intern/hair_volume.cpp

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

diff --git a/source/blender/physics/intern/hair_volume.cpp b/source/blender/physics/intern/hair_volume.cpp
index 4d10c04..3a0e702 100644
--- a/source/blender/physics/intern/hair_volume.cpp
+++ b/source/blender/physics/intern/hair_volume.cpp
@@ -550,43 +550,61 @@ bool BPH_hair_volume_solve_divergence(HairGrid *grid, float dt, float target_den
 	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 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 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];
+	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, u;
+	int i, j, k;
 	
-	BLI_assert(num_inner_cells >= 1);
+	#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_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;
+	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;
 				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];
+				float dx = NEIGHBOR_MARGIN_i1 ? 0.0f : vert_px->velocity[0] - v[0];
+				float dy = NEIGHBOR_MARGIN_j1 ? 0.0f : vert_py->velocity[1] - v[1];
+				float dz = NEIGHBOR_MARGIN_k1 ? 0.0f : vert_pz->velocity[2] - v[2];
 				
 				float divergence = (dx + dy + dz) * flowfac;
 				
@@ -610,20 +628,21 @@ bool BPH_hair_volume_solve_divergence(HairGrid *grid, float dt, float target_den
 	 * 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);
+	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_inner_cells, 7));
+	A.reserve(Eigen::VectorXi::Constant(num_cellsA, 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;
+	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 (vert->density > density_threshold) {
+				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;
@@ -635,30 +654,18 @@ bool BPH_hair_volume_solve_divergence(HairGrid *grid, float dt, float target_den
 					 * 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;
-					}
+					if (!NEIGHBOR_MARGIN_k0 && (vert - stride2)->density > density_threshold)
+						neighbor_lo_index[neighbors_lo++] = u - stride2;
+					if (!NEIGHBOR_MARGIN_j0 && (vert - stride1)->density > density_threshold)
+						neighbor_lo_index[neighbors_lo++] = u - stride1;
+					if (!NEIGHBOR_MARGIN_i0 && (vert - stride0)->density > density_threshold)
+						neighbor_lo_index[neighbors_lo++] = u - stride0;
+					if (!NEIGHBOR_MARGIN_i1 && (vert + stride0)->density > density_threshold)
+						neighbor_hi_index[neighbors_hi++] = u + stride0;
+					if (!NEIGHBOR_MARGIN_j1 && (vert + stride1)->density > density_threshold)
+						neighbor_hi_index[neighbors_hi++] = u + stride1;
+					if (!NEIGHBOR_MARGIN_k1 && (vert + stride2)->density > density_threshold)
+						neighbor_hi_index[neighbors_hi++] = u + stride2;
 					
 					/*int liquid_neighbors = neighbors_lo + neighbors_hi;*/
 					non_solid_neighbors = 6;
@@ -686,20 +693,23 @@ bool BPH_hair_volume_solve_divergence(HairGrid *grid, float dt, float target_den
 	
 	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;
+		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 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;
+						grad_p[0] = p0 - p[u - stride0];
+						grad_p[1] = p0 - p[u - stride1];
+						grad_p[2] = p0 - p[u - stride2];
 						
 						/* pressure gradient describes velocity delta */
 						madd_v3_v3v3fl(vert->velocity_smooth, vert->velocity, grad_p, inv_flowfac);
@@ -715,7 +725,7 @@ bool BPH_hair_volume_solve_divergence(HairGrid *grid, float dt, float target_den
 	}
 	else {
 		/* Clear result in case of error */
-		for (i = inner_cells_start, vert = grid->verts + inner_cells_start; i < num_inner_cells; ++i, ++vert) {
+		for (i = 0, vert = grid->verts; i < num_cells; ++i, ++vert) {
 			zero_v3(vert->velocity_smooth);
 		}




More information about the Bf-blender-cvs mailing list