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;
_______________________________________________
Bf-blender-cvs mailing list
[email protected]
http://lists.blender.org/mailman/listinfo/bf-blender-cvs
