Derivation of the Discretization

• We begin by defining our mesh indexing.

• We use both global and local indexing.

• Local indices enable us to uniquely define certain quantities that are associated with a vertex or face center and a cell.

• The following figures illustrate the indexing.
  Figure 1: Global indexes for four vertices associated with cell (i, j, k) .
  
  
  Figure 2: Local and global indices for three of six face centers associated with cell (i, j, k) .
  
  
  Figure 3: Local and global indices for three of six face centers associated with cell (i, j, k) .
  
  
  Figure 4: Vertex shared by the Right, Top, and Up faces having local index RTU.
  
  

• The intensities (scalars) are defined to exist at both cell center: $\phi^{{C}}_{{i,j,k}}$ , and on the face centers: $\phi^{{L}}_{{i,j,k}}$ , $\phi^{{R}}_{{i,j,k}}$ , $\phi^{{B}}_{{i,j,k}}$ , $\phi^{{T}}_{{i,j,k}}$ , $\phi^{{D}}_{{i,j,k}}$ , $\phi^{{U}}_{{i,j,k}}$ .

• As previously noted, the use of local indices implies that a quantity is uniquely associated with a single cell. For instance, unless it is otherwise stated, one should assume that $\phi^{{R}}_{{i,j,k}}$ $\neq$ $\phi^{{L}}_{{i+1,j,k}}$ .

• Vectors are defined in terms of face-area components located at the face centers: f^L_{i, j, k} , f^R_{i, j, k} , f^B_{i, j, k} , f^T_{i, j, k} , f^D_{i, j, k} , f^U_{i, j, k} , where f^L_{i, j, k} denotes the dot product of $\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$}$ with the outward-directed area vector located at the center of the left face of cell i, j, k . The other face-area components are defined analogously.

• The area vector is defined as the integral of the outward-directed unit normal vector over the face, i.e.,
  $$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{A}$}$$ = $$\oint$$$$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{n}$}$$ dA   ,
  where $\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{n}$}$ is a unit vector that is normal to the face at each point on the face.

• The average outward-directed unit normal vector for the face is defined as follows:
  $$\left\langle\vphantom{ \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{n}$} }\right.$$$$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{n}$}$$$$\left.\vphantom{ \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{n}$} }\right\rangle$$ = $${\frac{{\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{A}$}}}{{\Vert \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{A}$} \Vert}}}$$   ,
  where |$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{A}$}$| denotes the magnitude (standard Euclidian norm) of $\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{A}$}$ .

• This definition can be used to convert face-area flux components to face-normal components if desired, e.g.
  

  $$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$} \left\langle\vphantom{ \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{n}$} }\right. \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{n}$} \left.\vphantom{ \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{n}$} }\right\rangle \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$} {\frac{{\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{A}$}}}{{\Vert \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{A}$} \Vert}}}$$ =

  $${\frac{{f}}{{\Vert \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{A}$}\Vert}}}$$ =

  
  

• Note that |$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{A}$}$| is equal to the face area only when the face is flat.

• Interestingly, the true face areas never arise in our discretization scheme.

• Since it takes three components to define a full vector, the full vectors are considered to be located at the cell vertices: $\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$}$^LBD_{i, j, k} , $\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$}$^RBD_{i, j, k} , $\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$}$^LTD_{i, j, k} , $\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$}$^RTD_{i, j, k} , $\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$}$^LBU_{i, j, k} , $\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$}$^RBU_{i, j, k} , $\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$}$^LTU_{i, j, k} , $\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$}$^RTU_{i, j, k} .

• Each vertex vector is constructed using the face-area components and area vectors associated with the three faces that share that vertex. For instance,
  $$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$}$$^LBD_{i, j, k} = $${\frac{{f^L \left( \mbox{$\stackrel{^{\mathstrut}\smash{\longrigh... ...\times \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{A}$}^D \right)}}}$$ +
  $${\frac{{f^B \left( \mbox{$\stackrel{^{\mathstrut}\smash{\longrigh... ...\times \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{A}$}^L \right)}}}$$ + $${\frac{{f^D \left( \mbox{$\stackrel{^{\mathstrut}\smash{\longrigh... ...\times \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{A}$}^B \right)}}}$$   .

• It is convenient for our purposes to define an algebraic flux vector, $\hat{{F}}$ , consisting of the three face-area components associated with the physical vector, $\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$}$ , e.g.,
  $$\hat{{F}}_{{LBD}}^{}$$ = $$\left(\vphantom{ f^{L}_{i,j,k}, \, f^{B}_{i,j,k}, \, f^{D}_{i,j,k} }\right.$$f^L_{i, j, k}, f^B_{i, j, k}, f^D_{i, j, k}$$\left.\vphantom{ f^{L}_{i,j,k}, \, f^{B}_{i,j,k}, \, f^{D}_{i,j,k} }\right)^{t}_{}$$   ,
  where a superscript ``t'' denotes ``transpose.'' The three face-area components associated with the Right-Top-Up vertex are illustrated in Fig.4. The algebraic flux vectors for the other vertices are defined analogously.

• Our next step is to discretize the integral identity over a single cell.

• We first discretize the surface integral:
  $$\oint_{{\partial V}}^{}$$$$\phi$$$$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{H}$}$$ ^. $$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{n}$}$$  dA $$\approx$$ $$\sum_{{f}}^{}$$$$\phi^{{f}}_{}$$h^f   ,
  where f is the face index and the sum is taken over all faces.

• Next we approximate the flux volumetric integral:
  $$\int_{V}^{}$$ - D^-1$$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{H}$}$$ ^. D$$\mbox{$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{\nabla}$}$}$$$$\phi$$  dV $$\approx$$
  $$\sum_{{v}}^{}$$D^-1$$\left(\vphantom{ \mbox{$\stackrel{^{\mathstrut}\smash{\longrighta... ... \cdot \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$}^{v} }\right.$$$$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{H}$}$$^{v .} $$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$}$$^v$$\left.\vphantom{ \mbox{$\stackrel{^{\mathstrut}\smash{\longrighta... ... \cdot \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$}^{v} }\right)$$V^v  ,
  where v is the vertex index, and V^v denotes the volumetric weight associated with vertex v .

• The vertex volumetric weight is a free parameter in the support-operators method. We have tried several different weight definitions.

• The best choice appears to be one-eighth the triple product of the three edge vectors that share the vertex. For instance, the volumetric weight for the Left-Bottom-Down vertex can be expressed as follows using the point numbering shown in Fig.2:
  V^LBD = $${\frac{{1}}{{8}}}$$$$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{R}$}$$_{1, 2} x $$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{R}$}$$_{1, 3} ^. $$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{R}$}$$_{1, 4}   ,
  where $\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{R}$}$_{i, j} denotes the vector from vertex i to vertex j .

• Since these weights do not necessarily sum to the total hexahedral volume, we normalize them to do so.

• Although the expression given previously for the physical flux vector can be used to evaluate the dot products in the flux volumetric integral, we find it preferable to evaluate them in terms of the algebraic face-area flux vectors.

• This is achieved by transforming the face-area vectors to Cartesian vectors and then taking the dot product.

• There is a separate transformation matrix for each vertex. Let us denote this matrix for the Left-Bottom-Down vertex by A^LBD . Although we cannot explicitly define this matrix, we can explicitly define its inverse.

• This inverse matrix transforms Cartesian vectors to face-area vectors:
  $$\left[\vphantom{ {\bf A}^{LBD} }\right.$$$$\bf A^{{LBD}}_{}$$$$\left.\vphantom{ {\bf A}^{LBD} }\right]^{{-1}}_{}$$$$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{H}$}$$^LBD = $$\hat{{H}}^{{LBD}}_{}$$   ,
  where $\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{H}$}$ denotes a Left-Bottom-Down Cartesian flux vector,
  $$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{H}$}$$ = $$\left(\vphantom{ h^x, h^y, h^z }\right.$$h^x, h^y, h^z$$\left.\vphantom{ h^x, h^y, h^z }\right)^{t}_{}$$   ,
  $\hat{{H}}$ denotes a Left-Bottom-Down face-area flux vector,
  $$\hat{{H}}$$ = $$\left(\vphantom{ h^L, h^B, h^D }\right.$$h^L, h^B, h^D$$\left.\vphantom{ h^L, h^B, h^D }\right)^{t}_{}$$   ,
  and
  $$\left[\vphantom{ {\bf A}^{LBD} }\right.$$$$\bf A^{{LBD}}_{}$$$$\left.\vphantom{ {\bf A}^{LBD} }\right]^{{-1}}_{}$$ = $$\left[\vphantom{ \begin{array}{ccc} a^L_x & a^L_y & a^L_z \\ a^B_x & a^B_y & a^B_z \\ a^D_x & a^D_y & a^D_z \end{array} }\right.$$$$\begin{array}{ccc} a^L_x & a^L_y & a^L_z \\ a^B_x & a^B_y & a^B_z \\ a^D_x & a^D_y & a^D_z \end{array}$$$$\left.\vphantom{ \begin{array}{ccc} a^L_x & a^L_y & a^L_z \\ a^B_x & a^B_y & a^B_z \\ a^D_x & a^D_y & a^D_z \end{array} }\right]$$   ,
  where a^L_x denotes the x-component of the area vector associated with the left face. The remaining components of the matrix are defined analogously.

• To obtain A , we numerically invert A^-1 .

• We can now rewrite the approximation to the flux volumetric integral in terms of the face-area flux vectors as follows:
  $$\int_{V}^{}$$ - D^-1$$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{H}$}$$ ^. $$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$}$$  dV $$\approx$$
  $$\sum_{{v}}^{}$$D^-1$$\left(\vphantom{ \hat{H}^{v} \cdot \mathbf{S}^{v} \hat{F}^{v} }\right.$$$$\hat{{H}}^{{v}}_{}$$ ^. S^v$$\hat{{F}}^{{v}}_{}$$$$\left.\vphantom{ \hat{H}^{v} \cdot \mathbf{S}^{v} \hat{F}^{v} }\right)$$V^v  ,
  where
  S = A^tA   ,
  and the dot product is taken in the usual way.

• Finally, we approximate the divergence volumetric integral:
  $$\int_{V}^{}$$$$\phi$$$$\mbox{$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{\nabla}$} \cdot$}$$$$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{H}$}$$  dV $$\approx$$ $$\phi^{{C}}_{}$$$$\sum_{{f}}^{}$$h^f   ,
  where the sum is taken over all faces.

• Putting all of the pieces together, we obtain obtain the discretized integral identity:
  $$\sum_{{f}}^{}$$$$\phi^{{f}}_{}$$h^f + $$\sum_{{v}}^{}$$D^-1$$\left(\vphantom{ \hat{H}^{v} \cdot \mathbf{S}^{v} \hat{F}^{v} }\right.$$$$\hat{{H}}^{{v}}_{}$$ ^. S^v$$\hat{{F}}^{{v}}_{}$$$$\left.\vphantom{ \hat{H}^{v} \cdot \mathbf{S}^{v} \hat{F}^{v} }\right)$$V^v = $$\phi^{{C}}_{}$$$$\sum_{{f}}^{}$$h^f  .

• This identity must hold for all $\hat{{H}}$ . This requirement uniquely determines the six face-area components of the flux operator in terms of the cell-center and face-center intensities. In particular the equation for a given flux component is obtained by setting that component to unity while setting the remaining five components to zero. A matrix equation of the following form is obtained:
  W^-1$$\hat{{\mathcal{F}}}$$ = $$\Delta$$$$\hat{{\Phi}}$$   ,
  where
  $$\hat{{\mathcal{F}}}$$ = $$\left(\vphantom{ f^L, f^R, f^B, f^T, f^D, f^U }\right.$$f^L, f^R, f^B, f^T, f^D, f^U$$\left.\vphantom{ f^L, f^R, f^B, f^T, f^D, f^U }\right)^{{t}}_{}$$   ,
  and
  $$\Delta$$$$\hat{{\Phi}}$$ = $$\left(\vphantom{ \phi^C -\phi^L, \phi^C - \phi^R, \ldots , \phi^C - \phi^U }\right.$$$$\phi^{C}_{}$$ - $$\phi^{L}_{}$$,$$\phi^{C}_{}$$ - $$\phi^{R}_{}$$,...,$$\phi^{C}_{}$$ - $$\phi^{U}_{}$$$$\left.\vphantom{ \phi^C -\phi^L, \phi^C - \phi^R, \ldots , \phi^C - \phi^U }\right)^{{t}}_{}$$  .

• Numerically inverting the 6 x 6 matrix, W^-1 , we obtain the desired expressions for the six face-area flux components associated with the cell:
  $$\hat{{\mathcal{F}}}$$ = W$$\Delta$$$$\hat{{\Phi}}$$  .

• Having enforced the discrete identity over each single cell, we next enforce it over the entire mesh simply by requiring continuity of the intensity and the flow at each cell face on the mesh interior. For instance, at the interior face (i + ${\frac{{1}}{{2}}}$, j, k) , requiring continuity of the intensity yields:
  $$\phi^{R}_{{i,j,k}}$$ = $$\phi^{L}_{{i+1,j,k}}$$ = $$\phi_{{i+\frac{1}{2},j,k}}^{}$$  ,
  while requiring continuity of the flow yields:
  - f^R_{i, j, k} - f^L_{i+1, j, k} = 0  ,

• Eliminating the fluxes via the W -matrices and eliminating half of the face-center intensities via the continuity requirement, we are left with one intensity unknown at each cell center, and one intensity unknown with a unique global index at each face center.

• The equation for the intensity at the center of cell (i, j, k) is the balance equation for cell (i, j, k) :
  $${\frac{{\partial \phi^C}}{{\partial t}}}$$V + f^L + f^R + f^B + f^T + f^D + f^U = Q^CV   ,
  where V denotes the total cell volume, Q_C denotes the cell-center inhomogeneous source, and the cell index (i, j, k) has been supressed for simplicity.

• The equation for each face-center intensity on the mesh interior is a continuity-of-flow equation. For instance, the equation for $\phi_{{i+\frac{1}{2},j,k}}^{}$ is:
  - f^R_{i, j, k} - f^L_{i+1, j, k} = 0  .

• The equation for each face-center intensities on the outer mesh boundary is also a continuity-of-flux equation, but an extrapolated boundary condition is used to define the face-area flux component in the "ghost cell" adjacent to each boundary cell. For instance, the equation for $\phi_{{\frac{1}{2},j,k}}^{}$ takes the following form with a Marshak-type extrapolated boundary condition:
  - f^R_{0, j, k} - f^L_{1, j, k} = 0  ,
  where
  - f^R_{0, j, k} = $${\frac{{1}}{{2}}}$$|$$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{A}$}$$|^L_{1, j, k}$$\left(\vphantom{ \phi_{\frac{1}{2},j,k} - \phi^e_{\frac{1}{2},j,k} }\right.$$$$\phi_{{\frac{1}{2},j,k}}^{}$$ - $$\phi^{e}_{{\frac{1}{2},j,k}}$$$$\left.\vphantom{ \phi_{\frac{1}{2},j,k} - \phi^e_{\frac{1}{2},j,k} }\right)$$  ,
  where $\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{A}$}$^L_{1, j, k} is the Left area vector associated with cell (1, j, k) , and $\phi^{e}_{{\frac{1}{2},j,k}}$ is the extrapolated intensity associated with the boundary face.

• This completes the derivation of our discretization scheme.