// This is case2, the actual container with 12 drums inside it
function D = 0.062;
function C0 = x*y*z;
function al = 0.1;

double dt = 0.1;
double i;

vector n = (15,40,15); /* This works for sphere on left top of cyl 1 */
//vector n = (10,10,20); /* Amazing!! :) This worked!! for the lil cyl in front bottom of cyl 1 */
//vector n = (20,20,60); /* Yo!! This works too, just takes much longer */

vector a = (0,0,0);
vector b = (2,6,2);
mesh M = structured(n,a,b);

scene S = pov("case2conv.pov");
domain O = domain(S, outside(<1,0,0>) and outside(<0,1,1>) and outside(<0.9,1,1>) and outside(<0.8,0.9,1>) and outside(<0.7,0.8,1>) and outside(<0.6,0.7,1>) and outside(<0.5,0.6,1>) and outside(<0.4,0.5,1>) and outside(<0.3,0.4,1>) and outside(<0.2,0.3,1>) and outside(<0.1,0.2,1>) and outside(<0,0.1,1>) and outside(<0,0,1>));

// SDP: "optimization here" 
mesh Omesh = surface(O,M);
mesh Sphere = extract(Omesh,1); // 1 because <1,0,0> is the first color for domain definition

for (i=0; i<5; i++)
{
	cout << "\nIteration no. " << i+1 << "\n";
	solve(C) in O by M
	{
		pde(C)
			C - div(D*dt*grad(C)) = C0; /* Equation: dC/dt = D * d/dx(dC/dx) */
		dnu(C) = 0 on M;
	       	al*C + dnu(C) = 0.5 on Sphere;
	};
	C0 = C;
}
save(medit, "case2test-5", M);
save(medit, "case2test-5", C, M);