From e73df249c7f0e5acfcb187e6eb86246b44ae06f9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lukas=20T=C3=B6nne?= Date: Fri, 14 Nov 2014 10:42:09 +0100 Subject: 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. --- source/blender/physics/intern/hair_volume.cpp | 132 ++++++++++++++------------ 1 file changed, 71 insertions(+), 61 deletions(-) (limited to 'source') diff --git a/source/blender/physics/intern/hair_volume.cpp b/source/blender/physics/intern/hair_volume.cpp index 4d10c042f46..3a0e70269fc 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); } -- cgit v1.2.3