Fix T43220, T47551: collider scaling or rotation causes smoke to explode.

The problem happens because smoke collides only with the surface of the
collider and uses incompressible fluid solver. This means that scaling
the collider tries to compress or decompress fluid within the volume of
the collider, which can't be handled by the simulation. Fast rotation
likely also causes transient scaling due to emulation of arcs by chords.

This can be fixed by finding compartments completely isolated by obstacles
from the rest of the domain, and forcing total divergence within each one
to be zero so that equations are solvable. Physical validity is somewhat
dubious, but without this the solver simply breaks down.

From the physics point of view, the effect of the correction should be
similar to opening a hole from every cell to another dimension that lets
an equal amount of air to pass through to balance the change in volume.

Reviewers: miikah, lukastoenne

Reviewed By: lukastoenne

Subscribers: dafassi, scorpion81, #physics

Maniphest Tasks: T43220, T47551

Differential Revision: https://developer.blender.org/D2112

intern/smoke/intern/FLUID_3D.cpp

@@ -993,6 +993,9 @@ void FLUID_3D::project()
 
 	copyBorderAll(_pressure, 0, _zRes);
 
+	// fix fluid compression caused in isolated components by obstacle movement
+	fixObstacleCompression(_divergence);
+
 	// solve Poisson equation
 	solvePressurePre(_pressure, _divergence, _obstacles);
 
@@ -1291,6 +1294,142 @@ void FLUID_3D::setObstacleBoundaries(float *_pressure, int zBegin, int zEnd)
 	}	// z-loop
 }
 
+void FLUID_3D::floodFillComponent(int *buffer, size_t *queue, size_t limit, size_t pos, int from, int to)
+{
+	/* Flood 'from' cells with 'to' in the grid. Rely on (from != 0 && from != to && edges == 0) to stop. */
+	int offsets[] = { -1, +1, -_xRes, +_xRes, -_slabSize, +_slabSize };
+	size_t qend = 0;
+
+	buffer[pos] = to;
+	queue[qend++] = pos;
+
+	for (size_t qidx = 0; qidx < qend; qidx++)
+	{
+		pos = queue[qidx];
+
+		for (int i = 0; i < 6; i++)
+		{
+			size_t next = pos + offsets[i];
+
+			if (next < limit && buffer[next] == from)
+			{
+				buffer[next] = to;
+				queue[qend++] = next;
+			}
+		}
+	}
+}
+
+void FLUID_3D::mergeComponents(int *buffer, size_t *queue, size_t cur, size_t other)
+{
+	/* Replace higher value with lower. */
+	if (buffer[other] < buffer[cur])
+	{
+		floodFillComponent(buffer, queue, cur, cur, buffer[cur], buffer[other]);
+	}
+	else if (buffer[cur] < buffer[other])
+	{
+		floodFillComponent(buffer, queue, cur, other, buffer[other], buffer[cur]);
+	}
+}
+
+void FLUID_3D::fixObstacleCompression(float *divergence)
+{
+	int x, y, z;
+	size_t index;
+
+	/* Find compartments completely separated by obstacles.
+	 * Edge of the domain is automatically component 0. */
+	int *component = new int[_totalCells];
+	size_t *queue = new size_t[_totalCells];
+
+	memset(component, 0, sizeof(int) * _totalCells);
+
+	int next_id = 1;
+
+	for (z = 1, index = _slabSize + _xRes + 1; z < _zRes - 1; z++, index += 2 * _xRes)
+	{
+		for (y = 1; y < _yRes - 1; y++, index += 2)
+		{
+			for (x = 1; x < _xRes - 1; x++, index++)
+			{
+				if(!_obstacles[index])
+				{
+					/* Check for connection to the domain edge at iteration end. */
+					if ((x == _xRes-2 && !_obstacles[index + 1]) ||
+					    (y == _yRes-2 && !_obstacles[index + _xRes]) ||
+					    (z == _zRes-2 && !_obstacles[index + _slabSize]))
+					{
+						component[index] = 0;
+					}
+					else {
+						component[index] = next_id;
+					}
+
+					if (!_obstacles[index - 1])
+						mergeComponents(component, queue, index, index - 1);
+					if (!_obstacles[index - _xRes])
+						mergeComponents(component, queue, index, index - _xRes);
+					if (!_obstacles[index - _slabSize])
+						mergeComponents(component, queue, index, index - _slabSize);
+
+					if (component[index] == next_id)
+						next_id++;
+				}
+			}
+		}
+	}
+
+	delete[] queue;
+
+	/* Compute average divergence within each component. */
+	float *total_divergence = new float[next_id];
+	int *component_size = new int[next_id];
+
+	memset(total_divergence, 0, sizeof(float) * next_id);
+	memset(component_size, 0, sizeof(int) * next_id);
+
+	for (z = 1, index = _slabSize + _xRes + 1; z < _zRes - 1; z++, index += 2 * _xRes)
+	{
+		for (y = 1; y < _yRes - 1; y++, index += 2)
+		{
+			for (x = 1; x < _xRes - 1; x++, index++)
+			{
+				if(!_obstacles[index])
+				{
+					int ci = component[index];
+
+					component_size[ci]++;
+					total_divergence[ci] += divergence[index];
+				}
+			}
+		}
+	}
+
+	/* Adjust divergence to make the average zero in each component except the edge. */
+	total_divergence[0] = 0.0f;
+
+	for (z = 1, index = _slabSize + _xRes + 1; z < _zRes - 1; z++, index += 2 * _xRes)
+	{
+		for (y = 1; y < _yRes - 1; y++, index += 2)
+		{
+			for (x = 1; x < _xRes - 1; x++, index++)
+			{
+				if(!_obstacles[index])
+				{
+					int ci = component[index];
+
+					divergence[index] -= total_divergence[ci] / component_size[ci];
+				}
+			}
+		}
+	}
+
+	delete[] component;
+	delete[] component_size;
+	delete[] total_divergence;
+}
+
 //////////////////////////////////////////////////////////////////////
 // add buoyancy forces
 //////////////////////////////////////////////////////////////////////
@@ -1650,4 +1789,4 @@ void FLUID_3D::updateFlame(float *react, float *flame, int total_cells)
 		else
 			flame[index] = 0.0f;
 	}
-}
+}

intern/smoke/intern/FLUID_3D.h

@@ -195,6 +195,8 @@ struct FLUID_3D
 		void setObstacleBoundaries(float *_pressure, int zBegin, int zEnd);
 		void setObstaclePressure(float *_pressure, int zBegin, int zEnd);
 
+		void fixObstacleCompression(float *divergence);
+
 	public:
 		// advection, accessed e.g. by WTURBULENCE class
 		//void advectMacCormack();
@@ -202,6 +204,9 @@ struct FLUID_3D
 		void advectMacCormackEnd1(int zBegin, int zEnd);
 		void advectMacCormackEnd2(int zBegin, int zEnd);
 
+		void floodFillComponent(int *components, size_t *queue, size_t limit, size_t start, int from, int to);
+		void mergeComponents(int *components, size_t *queue, size_t cur, size_t other);
+
 		/* burning */
 		float *_burning_rate; // RNA pointer
 		float *_flame_smoke; // RNA pointer