(* ::Package:: *)

(************************************************************************)
(* This file was generated automatically by the Mathematica front end.  *)
(* It contains Initialization cells from a Notebook file, which         *)
(* typically will have the same name as this file except ending in      *)
(* ".nb" instead of ".m".                                               *)
(*                                                                      *)
(* This file is intended to be loaded into the Mathematica kernel using *)
(* the package loading commands Get or Needs.  Doing so is equivalent   *)
(* to using the Evaluate Initialization Cells menu command in the front *)
(* end.                                                                 *)
(*                                                                      *)
(* DO NOT EDIT THIS FILE.  This entire file is regenerated              *)
(* automatically each time the parent Notebook file is saved in the     *)
(* Mathematica front end.  Any changes you make to this file will be    *)
(* overwritten.                                                         *)
(************************************************************************)



DiffMat1DPeriodic[nx_] := SparseArray[{
Band[{1,1}]->-2.0,
Band[{1,2}]->1.0,
Band[{2,1}]->1.0,
{1,nx+1}->1.0,
{nx+1,1}->1.0
},
{nx+1, nx+1}]
MatrixForm[DiffMat1DPeriodic[6]]


Clear[Relax3,\[Omega],\[Nu]x,\[Nu]y,\[Nu]z,U,F]
SetAttributes[Relax3,HoldRest]
Relax3 [{\[Omega]_,{\[Nu]x_,\[Nu]y_,\[Nu]z_}},U_,F_]:=Module[
{nx,ny,nz,i,j,k,\[Alpha],\[Beta]},
{nx,ny,nz}=Dimensions[F];
\[Alpha]=0.5\[Omega]/(\[Nu]x+\[Nu]y+\[Nu]z);
\[Beta] = 1-\[Omega];
(* red cells i+j+k even first *) 
Do[
U[[i,j,k]]=\[Alpha] (F[[i,j,k]] 
+ \[Nu]x (U[[i-1,j,k]]+U[[i+1,j,k]])
+ \[Nu]y (U[[i,j-1,k]]+U[[i,j+1,k]])
+ \[Nu]z (U[[i,j,k-1]]+U[[i,j,k+1]])
 ) + \[Beta] U[[i,j,k]],
{i,2,nx-1,2}, (* even *)
{j,2,ny-1,2}, (* even *)
{k,2,nz-1,2}     (* even *)
];
Do[
U[[i,j,k]]=\[Alpha] (F[[i,j,k]] 
+ \[Nu]x (U[[i-1,j,k]]+U[[i+1,j,k]])
+ \[Nu]y (U[[i,j-1,k]]+U[[i,j+1,k]])
+ \[Nu]z (U[[i,j,k-1]]+U[[i,j,k+1]])
 ) + \[Beta] U[[i,j,k]],
{i,2,nx-1,2},  (* even *)
{j,3,ny-1,2},  (* odd *)
{k,3,nz-1,2}     (* odd *)
];
Do[
U[[i,j,k]]=\[Alpha] (F[[i,j,k]] 
+ \[Nu]x (U[[i-1,j,k]]+U[[i+1,j,k]])
+ \[Nu]y (U[[i,j-1,k]]+U[[i,j+1,k]])
+ \[Nu]z (U[[i,j,k-1]]+U[[i,j,k+1]])
 ) + \[Beta] U[[i,j,k]],
{i,3,nx-1,2},  (* odd *)
{j,2,ny-1,2},  (* even *)
{k,3,nz-1,2}     (* odd *)
];
Do[
U[[i,j,k]]=\[Alpha] (F[[i,j,k]] 
+ \[Nu]x (U[[i-1,j,k]]+U[[i+1,j,k]])
+ \[Nu]y (U[[i,j-1,k]]+U[[i,j+1,k]])
+ \[Nu]z (U[[i,j,k-1]]+U[[i,j,k+1]])
 ) + \[Beta] U[[i,j,k]],
{i,3,nx-1,2},  (* odd *)
{j,3,ny-1,2},  (* odd *)
{k,2,nz-1,2}     (* even *)
];
(* Now black Cells i+j+k odd *)
Do[
U[[i,j,k]]=\[Alpha] (F[[i,j,k]] 
+ \[Nu]x (U[[i-1,j,k]]+U[[i+1,j,k]])
+ \[Nu]y (U[[i,j-1,k]]+U[[i,j+1,k]])
+ \[Nu]z (U[[i,j,k-1]]+U[[i,j,k+1]])
 ) + \[Beta] U[[i,j,k]],
{i,3,nx-1,2}, (* odd *)
{j,2,ny-1,2}, (* even *)
{k,2,nz-1,2}     (* even *)
];
Do[
U[[i,j,k]]=\[Alpha] (F[[i,j,k]] 
+ \[Nu]x (U[[i-1,j,k]]+U[[i+1,j,k]])
+ \[Nu]y (U[[i,j-1,k]]+U[[i,j+1,k]])
+ \[Nu]z (U[[i,j,k-1]]+U[[i,j,k+1]])
 ) + \[Beta] U[[i,j,k]],
{i,2,nx-1,2}, (* even *)
{j,3,ny-1,2}, (* odd *)
{k,2,nz-1,2}     (* even *)
];
Do[
U[[i,j,k]]=\[Alpha] (F[[i,j,k]] 
+ \[Nu]x (U[[i-1,j,k]]+U[[i+1,j,k]])
+ \[Nu]y (U[[i,j-1,k]]+U[[i,j+1,k]])
+ \[Nu]z (U[[i,j,k-1]]+U[[i,j,k+1]])
 ) + \[Beta] U[[i,j,k]],
{i,2,nx-1,2}, (* even *)
{j,2,ny-1,2}, (* even *)
{k,3,nz-1,2}     (* odd *)
];
Do[
U[[i,j,k]]=\[Alpha] (F[[i,j,k]] 
+ \[Nu]x (U[[i-1,j,k]]+U[[i+1,j,k]])
+ \[Nu]y (U[[i,j-1,k]]+U[[i,j+1,k]])
+ \[Nu]z (U[[i,j,k-1]]+U[[i,j,k+1]])
 ) + \[Beta] U[[i,j,k]],
{i,3,nx-1,2}, (* odd *)
{j,3,ny-1,2}, (* odd *)
{k,3,nz-1,2}     (* odd *)
];
]



