(*global graphics parameters*)
pointsize=.06;
opacity=.5;

(*mesh parameters*)
domain={{0,1},{0,1},{0,1}};
α=1;
Δt=.0005; (*time step parameter*)
Δx=.2; (*both x and y direction mesh parameter*)
λ=((α*Δt)/((Δx)^2));
(*generate node records and adjacency matrix*)
xcoords=Table[domain[[1,1]]+i*Δx,{i,0,(domain[[1,2]]-domain[[1,1]])/Δx}];
ycoords=Table[domain[[2,1]]+i*Δx,{i,0,(domain[[2,2]]-domain[[2,1]])/Δx}];
zcoords=Table[domain[[3,1]]+i*Δx,{i,0,(domain[[3,2]]-domain[[3,1]])/Δx}];
max=Length[xcoords]Length[ycoords]Length[zcoords];
(*to determine the current position based on node: first, determine z-plane (how far back), then determine column (how far to the right), then determine row (how far up)*)
nodeRecord=Table[{i,xcoords[[Floor[Mod[i-1,Length[xcoords]Length[ycoords]]/Length[ycoords]]+1]],ycoords[[Mod[i,Length[ycoords],1]]],
zcoords[[Floor[(i-1)/(Length[xcoords]Length[ycoords])]+1]]
},
{i,1,Length[xcoords]Length[ycoords]Length[ycoords]}];
(*functions that return the node relations*)
above[i_]:=i+1;
below[i_]:=i-1;
left[i_]:=i-Length[ycoords];
right[i_]:=i+Length[ycoords];
front[i_]:=i-Length[xcoords]Length[ycoords];
back[i_]:=i+Length[xcoords]Length[ycoords];
(*nodes are arranged bottom left to upper right, front to back*)
(*load the matrix*)
A=(1-6λ)*IdentityMatrix[max];

(*corners first*)
i=1;(*bottom left, front face*)
A[[i,above[i]]]=λ;
A[[i,right[i]]]=λ;
A[[i,back[i]]]=λ;
i=Nest[above,1,Length[ycoords]-1];(*top left, front face*)
A[[i,below[i]]]=λ;
A[[i,right[i]]]=λ;
A[[i,back[i]]]=λ;
i=Nest[right,1,Length[xcoords]-1];(*bottom right, front face*)
A[[i,above[i]]]=λ;
A[[i,left[i]]]=λ;
A[[i,back[i]]]=λ;
i=Nest[right,Nest[above,1,Length[ycoords]-1],Length[xcoords]-1];(*top right, front face*)
A[[i,below[i]]]=λ;
A[[i,left[i]]]=λ;
A[[i,back[i]]]=λ;
i=Nest[back,1,Length[zcoords]-1];(*bottom left, back face*)
A[[i,above[i]]]=λ;
A[[i,right[i]]]=λ;
A[[i,front[i]]]=λ;
i=Nest[above,Nest[back,1,Length[zcoords]-1],Length[ycoords]-1];(*top left, back face*)
A[[i,below[i]]]=λ;
A[[i,right[i]]]=λ;
A[[i,front[i]]]=λ;
i=Nest[right,Nest[back,1,Length[zcoords]-1],Length[xcoords]-1];(*bottom right, back face*)
A[[i,above[i]]]=λ;
A[[i,left[i]]]=λ;
A[[i,front[i]]]=λ;
i=Nest[above,Nest[right,Nest[back,1,Length[zcoords]-1],Length[xcoords]-1],Length[ycoords]-1];(*upper right, back face*)
A[[i,left[i]]]=λ;
A[[i,below[i]]]=λ;
A[[i,front[i]]]=λ;
(*sides*)
(*left edge, front face*)
For[i=2,i<Length[ycoords],i++,
A[[i,above[i]]]=λ;
A[[i,below[i]]]=λ;
A[[i,right[i]]]=λ;
A[[i,back[i]]]=λ;
]
(*bottom edge, front face*)
For[i=2,i<Length[xcoords],i++,
k=Nest[right,1,i-1];(*move right starting from node 1*)
A[[k,left[k]]]=λ;
A[[k,right[k]]]=λ;
A[[k,above[k]]]=λ;
A[[k,back[k]]]=λ;
]
(*top edge, front face*)
For[i=2,i<Length[xcoords],i++,
k=Nest[right,Nest[above,1,Length[ycoords]-1],i-1];(*start at top left, move right*)
A[[k,left[k]]]=λ;
A[[k,right[k]]]=λ;
A[[k,below[k]]]=λ;
A[[k,back[k]]]=λ;
]
(*right edge, front face*)
For[i=2,i<Length[ycoords],i++,
k=Nest[above,Nest[right,1,Length[xcoords]-1],i-1];(*this is the form of every element on the right front edge*)
A[[k,left[k]]]=λ;
A[[k,below[k]]]=λ;
A[[k,above[k]]]=λ;
A[[k,back[k]]]=λ;
]
(*bottom left front-back edge*)
For[i=2,i<Length[zcoords],i++,
k=Nest[back,1,i-1];
A[[k,front[k]]]=λ;
A[[k,back[k]]]=λ;
A[[k,above[k]]]=λ;
A[[k,right[k]]]=λ;
]
(*top left front-back edge*)
For[i=2,i<Length[zcoords],i++,
k=Nest[back,Nest[above,1,Length[ycoords]-1],i-1];
A[[k,front[k]]]=λ;
A[[k,back[k]]]=λ;
A[[k,right[k]]]=λ;
A[[k,below[k]]]=λ;
]
(*bottom right front-back edge*)
For[i=2,i<Length[zcoords],i++,
k=Nest[back,Nest[right,1,Length[xcoords]-1],i-1];
A[[k,above[k]]]=λ;
A[[k,left[k]]]=λ;
A[[k,front[k]]]=λ;
A[[k,back[k]]]=λ;
]
(*top right front-back edge*)
For[i=2,i<Length[zcoords],i++,
k=Nest[back,Nest[above,Nest[right,1,Length[xcoords]-1],Length[ycoords]-1],i-1];
A[[k,left[k]]]=λ;
A[[k,below[k]]]=λ;
A[[k,front[k]]]=λ;
A[[k,back[k]]]=λ;
]
(*bottom edge, back face*)
For[i=2,i<Length[xcoords],i++,
k=Nest[right,Nest[back,1,Length[zcoords]-1],i-1];
A[[k,right[k]]]=λ;
A[[k,left[k]]]=λ;
A[[k,front[k]]]=λ;
A[[k,above[k]]]=λ;
]
(*top edge, back face*)
For[i=2,i<Length[xcoords],i++,
k=Nest[right,Nest[above,Nest[back,1,Length[zcoords]-1],Length[ycoords]-1],i-1];
A[[k,right[k]]]=λ;
A[[k,left[k]]]=λ;
A[[k,below[k]]]=λ;
A[[k,front[k]]]=λ;
]
(*left edge, back face*)
For[i=2,i<Length[ycoords],i++,
k=Nest[above,Nest[back,1,Length[zcoords]-1],i-1];
A[[k,above[k]]]=λ;
A[[k,below[k]]]=λ;
A[[k,front[k]]]=λ;
A[[k,right[k]]]=λ;
]
(*right edge, back face*)
For[i=2,i<Length[ycoords],i++,
k=Nest[above,Nest[right,Nest[back,1,Length[zcoords]-1],Length[xcoords]-1],i-1];
A[[k,above[k]]]=λ;
A[[k,below[k]]]=λ;
A[[k,left[k]]]=λ;
A[[k,front[k]]]=λ;
]
(*front face*)
For[i=2,i<Length[xcoords],i++,
For[j=2,j<Length[ycoords],j++,
k=Nest[above,Nest[right,1,i-1],j-1];
A[[k,above[k]]]=λ;
A[[k,below[k]]]=λ;
A[[k,right[k]]]=λ;
A[[k,left[k]]]=λ;
A[[k,back[k]]]=λ;
]
]
(*back face*)
For[i=2,i<Length[xcoords],i++,
For[j=2,j<Length[ycoords],j++,
k=Nest[above,Nest[right,Nest[back,1,Length[zcoords]-1],i-1],j-1];
A[[k,above[k]]]=λ;
A[[k,below[k]]]=λ;
A[[k,right[k]]]=λ;
A[[k,left[k]]]=λ;
A[[k,front[k]]]=λ;
]
]
(*left face*)
For[i=2,i<Length[ycoords],i++,
For[j=2,j<Length[zcoords],j++,
k=Nest[back,Nest[above,1,i-1],j-1];
A[[k,above[k]]]=λ;
A[[k,below[k]]]=λ;
A[[k,right[k]]]=λ;
A[[k,front[k]]]=λ;
A[[k,back[k]]]=λ;
]
]
(*right face*)
For[i=2,i<Length[ycoords],i++,
For[j=2,j<Length[zcoords],j++,
k=Nest[back,Nest[above,Nest[right,1,Length[xcoords]-1],i-1],j-1];
A[[k,above[k]]]=λ;
A[[k,below[k]]]=λ;
A[[k,front[k]]]=λ;
A[[k,back[k]]]=λ;
A[[k,left[k]]]=λ;
]
]
(*top face*)
For[i=2,i<Length[xcoords],i++,
For[j=2,j<Length[zcoords],j++,
k=Nest[back,Nest[right,Nest[above,1,Length[ycoords]-1],i-1],j-1];
A[[k,front[k]]]=λ;
A[[k,back[k]]]=λ;
A[[k,left[k]]]=λ;
A[[k,right[k]]]=λ;
A[[k,below[k]]]=λ;
]
]
(*bottom face*)
For[i=2,i<Length[xcoords],i++,
For[j=2,j<Length[zcoords],j++,
k=Nest[back,Nest[right,1,i-1],j-1];
A[[k,front[k]]]=λ;
A[[k,back[k]]]=λ;
A[[k,left[k]]]=λ;
A[[k,right[k]]]=λ;
A[[k,above[k]]]=λ;
]
]
(*middle nodes*)
For[i=2,i<Length[ycoords],i++,
For[j=2,j<Length[xcoords],j++,
For[k=2,k<Length[zcoords],k++,
n=(k-1)Length[xcoords]Length[ycoords]+(j-1)Length[ycoords]+i;
A[[n,above[n]]]=λ;
A[[n,below[n]]]=λ;
A[[n,front[n]]]=λ;
A[[n,back[n]]]=λ;
A[[n,left[n]]]=λ;
A[[n,right[n]]]=λ;
]
]
]

(*boundary conditions go here*)
(*leave blank for full conductivity outwards*)
(*ex. insulate all sides*)
For[i=1,i≤Length[A],i++,
A[[i,i]]+=1-Sum[A[[i,k]],{k,1,Length[A]}];
]
(*end boundary conditions*)

(*initial conditions*)
(*values=Table[Append[Delete[nodeRecord[[i]],1],RandomReal[1]],{i,1,max}];*)
values=Table[Append[Delete[nodeRecord[[i]],1],0],{i,1,max}];
values[[max,4]]=500;
(*display*)
output=Point[Table[Delete[values[[i]],4],{i,1,Length[values]}],VertexColors->Table[ColorData["ThermometerColors"][values[[i,4]]],{i,1,Length[values]}]];
Graphics3D[{Opacity[opacity],PointSize[pointsize],output}]

(*visualization 1:*)
(*generate table of plots for successive iterations*)
showIterations={1,400};(*show iterations from {start,finish}*)
showEvery=100;(*show every 5 iterations*)
record={Point[Table[Delete[values[[i]],4],{i,1,Length[values]}],VertexColors->Table[ColorData["ThermometerColors"][values[[i,4]]],{i,1,Length[values]}]]};
iteratedvalues=values;
For[i=1,i≤showIterations[[2]],i++,
temp=Transpose[iteratedvalues];
u=temp[[4]];
u=A.u;
temp[[4]]=u;
iteratedvalues=Transpose[temp];
record=Append[record,Point[Table[Delete[iteratedvalues[[i]],4],{i,1,Length[iteratedvalues]}],VertexColors->Table[ColorData["ThermometerColors"][iteratedvalues[[i,4]]],{i,1,Length[iteratedvalues]}]]];
]
Table[Graphics3D[{Opacity[opacity],PointSize[pointsize],record[[i]]}],{i,showIterations[[1]],showIterations[[2]],showEvery}]

(*visualization 2:*)
(*reset*)
iteratedvalues=values;

(*animation: run cell, click on the plot to stop animation*)
record={values};
DynamicModule[{go=True},
Monitor[While[go,
temp=Transpose[iteratedvalues];
u=temp[[4]];
u=A.u;
temp[[4]]=u;
iteratedvalues=Transpose[temp];
record=Append[record,Point[Table[Delete[iteratedvalues[[i]],4],{i,1,Length[iteratedvalues]}],VertexColors->Table[ColorData["ThermometerColors"][iteratedvalues[[i,4]]],{i,1,Length[iteratedvalues]}]]];
p=Graphics3D[{Opacity[opacity],PointSize[pointsize],Point[Table[Delete[iteratedvalues[[i]],4],{i,1,Length[iteratedvalues]}],VertexColors->Table[ColorData["ThermometerColors"][iteratedvalues[[i,4]]],{i,1,Length[iteratedvalues]}]]}];
],EventHandler[
p,{"MouseClicked":>(go=False)}]]
];
p

ListAnimate[Table[Graphics3D[{Opacity[opacity],PointSize[pointsize],record[[i]]}],{i,1,Length[record]}]]

Export["3dheat.avi",%]