(* ::Package:: *)

Print["initilisation3D.m is run"];
L=20;tol=0.01;
iini=-L;ilast=L;
jini=-L;jlast=L;
kini=-L;klast=L;

xini=0;xlast=1.0;yini=0;ylast=1.0;zini=0;zlast=1.0;
Delta=(xlast-xini)/L;
For[i=iini,i<=ilast,i++,
For[j=jini,j<=jlast,j++,
For[k=kini,k<=klast,k++,
V[i,j,k]=0.0;
rho[i,j,k]=0;
];
];(*end for i*)
];(*end for i*)

im=IntegerPart[(1/2)*(iini+ilast)];jm=IntegerPart[(1/2)*(jini+jlast)];km=IntegerPart[(1/2)*(kini+klast)];
rho[im,jm,km]=1.0/(Delta^3);

For[i=iini,i<=ilast,i++,
x[i]=xini+Delta*i;
];(*end for i*)
For[j=jini,j<=jlast,j++,
y[j]=yini+Delta*j;
];(*end for i*)

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

V[im,jm,km]=1;

(*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,1];*)
