(* ::Package:: *)

Print["update3D.m is run"]
s=200*tol;n=0;NN=(ilast-iini)*(jlast-jini)*(klast-kini);
(*
txyzVn[0]=txyzV;
*)
While[s>tol,
{n=n+1;

For [i=iini+1,i<=ilast-1,i++,
For [j=jini+1,j<=jlast-1,j++,
For [k=kini+1,k<=klast-1,k++,
Vold[i,j,k]=V[i,j,k];
];(*end for k*)
];(*end for j*)
];(*end for i*)

For [i=iini+1,i<=ilast-1,i++,
For [j=jini+1,j<=jlast-1,j++,
For [k=kini+1,k<=klast-1,k++,

(***** for potential field ***************)
V[i,j,k]=(1/6.0)*
(V[i-1,j,k]+V[i+1,j,k]+V[i,j-1,k]+V[i,j+1,k]+V[i,j,k+1]+V[i,j,k-1]+rho[i,j,k]*Delta^3);
];(*end for k*)
];(*end for j*)
];(*end for i*)

(*BC*)
(* xy *)
k=kini;
For [i=iini,i<=ilast,i++,
For [j=jini,j<=jlast,j++,
V[i,j,k]=0;
];];(*end for j*)

(* xybar *)
k=klast;
For [i=iini,i<=ilast,i++,
For [j=jini,j<=jlast,j++,
V[i,j,k]=0;
];];(*end for j*)

(* yz *)
i=iini;
For [j=jini,j<=jlast,j++,
For [k=kini,k<=klast,k++,
V[i,j,k]=0;
];];(*end for j*)

(* yzbar *)
i=ilast;
For [j=jini,j<=jlast,j++,
For [k=kini,k<=klast,k++,
V[i,j,k]=0;
];];(*end for j*)


(* zx *)
j=jini;
For [i=iini,i<=ilast,i++,
For [k=kini,k<=klast,k++,
V[i,j,k]=0;
];];(*end for j*)

(* zxbar *)
j=jlast;
For [i=iini,i<=ilast,i++,
For [k=kini,k<=klast,k++,
V[i,j,k]=0;
];];(*end for j*)
(*BC*)

(*
txyzV=Table[{x[i],y[j],z[k],V[i,j,k]},{k,kini,klast},{j,jini,jlast},{i,iini,ilast}];
txyzV=Flatten[txyzV,2];
txyzVn[n]=txyzV;
*)
txyV=Table[{x[i],y[j],V[i,j,km]},{j,jini,jlast},{i,iini,ilast}];
txyV=Flatten[txyV,1];
txyVn[n]=txyV;

(*Print[ListContourPlot[txyVn[n]]];*)
(*Print[ListDensityPlot[txyVn[n]],Mesh->True];*)
(*Print[ListPlot3D[txyVn[n]]];*)
(*Print[
ListPlot3D[Flatten[Table[{x[i],y[j],Vold[i,j,km]},{i,iini+1,ilast-1},{j,jini+1,jlast-1}],1]]
];*)
(***** for potential field ***************)


(***** for vector field ***************)
For [i=iini+1,i<=ilast-1,i++,
For [j=jini+1,j<=jlast-1,j++,
For [k=kini+1,k<=klast-1,k++,
Ex[i,j,k]=-(V[i+1,j,k]-V[i-1,j,k])/(2*Delta);
Ey[i,j,k]=-(V[i,j+1,k]-V[i,j-1,k])/(2*Delta);
Ez[i,j,k]=-(V[i,j,k+1]-V[i,j,k-1])/(2*Delta);
];];];

txyzE=Table[{Ex[i,j,k],Ey[i,j,k],Ez[i,j,k]},{k,kini+1,klast-1},{j,jini+1,jlast-1},{i,iini+1,ilast-1}];
(*txyzE=Flatten[txyzE,2];*)
txyzEn[n]=txyzE;

txyE=Table[{Ex[i,j,km],Ey[i,j,km]},{j,jini+1,jlast-1},{i,iini+1,ilast-1}];
(*txyE=Flatten[txyE,1];*)
txyEn[n]=txyE;

(*Print[txyEn[n]];*)
(*Print[ListVectorPlot[txyEn[n]]];*)
(***** for vector field ***************)

(* check for convergernce *)
If[n==1,s=2*tol,
s=Sqrt[Sum[(V[i,j,k]-Vold[i,j,k])^2,{i,iini+1,ilast-1},{j,jini+1,jlast-1},{k,kini+1,klast-1}]/NN];
];(*end if *)

Print["n, s= ", {n,s}];
nlast=n;
}];(*end while*)


Print[Manipulate[
ListDensityPlot[txyVn[n],Mesh->True]
,{n,1,nlast,1}
]];

(*
Print[Manipulate[
ListVectorPlot[txyEn[n]]
,{n,1,nlast,1}
]];


Print[Manipulate[
ListVectorPlot3D[txyzEn[n],Mesh->True]
,{n,1,nlast,1}
]];
*)
