The Support-Operators Method

In this section we describe the support-operators method. It is convenient at this point to define a flux operator given by - D$\mbox{$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{\nabla}$}$}$ . The diffusion operator of interest is given by the product of the divergence operator and the flux operator: - $\mbox{$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{\nabla}$} \cdot$}$D$\mbox{$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{\nabla}$}$}$ . The support-operators method is based upon the following three facts:

• Given appropriately defined scalar and vector inner products, the divergence and flux operators are adjoint to one another.
• The adjoint of an operator varies with the definition of its associated inner products, but is unique for fixed inner products.
• The product of an operator and its adjoint is a self-adjoint positive-definite operator.

The mathematical details relating to these facts are given in [8]. As explained in [8], the adjoint relationship between the flux and divergence operators is embodied in the following integral identity:

$$\oint_{{\partial V}}^{} \phi \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{H}$} \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{n}$} \int_{V}^{} \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{H}$} \mbox{$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{\nabla}$}$} \phi \int_{V}^{} \phi \mbox{$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{\nabla}$} \cdot$} \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{H}$}$$ (4)

where $\phi$ is an arbitrary scalar function, $\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{H}$}$ is an arbitrary vector function, V denotes a volume, $\partial$V denotes its surface, and $\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{n}$}$ denotes the outward-directed unit normal associated with that surface. Our support-operators method can be conceptually described in the simplest terms as follows:

1. Define discrete scalar and vector spaces to be used in a discretization of Eq. (4).
2. Fully discretize all but the flux operator in Eq. (4) over a single arbitrary cell. The flux operator is left in the general form of a discrete vector as defined in Step 1.
3. Solve for the discrete flux operator (i.e., for its vector components) on a single arbitrary cell by requiring that the discrete version of Eq. (4) hold for all elements of the discrete scalar and vector spaces defined in Step 1.
4. Obtain the interior-mesh discretization of Eq. (4) by connecting adjacent mesh cells in such a way as to ensure that Eq. (4) is satisfied over the whole grid. This simply amounts to enforcing continuity of intensity and flux at the cell interfaces.
5. Change the flux operator at those cell faces on the exterior mesh boundary so as to satisfy the appropriate boundary conditions.
6. Combine the global divergence matrix and the global flux matrix to obtain the global diffusion matrix.

The actual method is somewhat more complicated because of the presence of both cell-center and cell-face intensities, but this description nonetheless conveys the central theme of the method.

To make this process concrete, we next generate the diffusion matrix for a hexahedral mesh in Cartesian geometry. To simplify the presentation, we assume a logically-rectangular mesh. However, our discretization scheme can be used with unstructured meshes as well. The assumption of a logically-rectangular mesh merely simplifies our notation and mesh indexing. Our first step is to define that indexing. For reasons explained later, both global and local indices are used. Let us first consider the global indices. The cell centers carry integral global indices, e.g., (i, j, k) ; cell vertices carry half-integral global indices, e.g., (i + ${\frac{{1}}{{2}}}$, j + ${\frac{{1}}{{2}}}$, k + ${\frac{{1}}{{2}}}$) ; and face centers carry mixed global indices composed of both integral and half-integral indices, e.g., (i + ${\frac{{1}}{{2}}}$, j, k) . The global indices for four of the vertices associated with cell (i, j, k) are illustrated in Fig. 1.

Figure 1: Global indices for four vertices associated with cell (i, j, k) .

Local indices allow us to uniquely define certain quantities that are associated with a vertex or face center and a cell. For instance, the local indices for the six faces associated with each cell are given by L, R, B, T, D, and U, which denote Left, Right, Bottom, Top, Down, and Up respectively. This local face indexing is illustrated for cell (i, j, k) in Fig. 2 and Fig. 3 together with a mapping between the local indices and the corresponding global indices.

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) .

Note that the index i increases when moving from Left to Right, the index j increases when moving from Bottom to Top, and the index k increases when moving from the Down to Up. The local indices for the vertices follow directly from the face indices in that each vertex is uniquely shared by three faces of the cell. Thus the vertex shared by the Right, Top, and Up faces is denoted by the index RTU. This vertex is illustrated in Fig. 4.

Figure 4: Vertex shared by the Right, Top, and Up faces having local index RTU.

The vector and matrix notation used from this point forward in this paper is as follows. Each vector is denoted by an upper-case symbol and the components of that vector are denoted by the corresponding lower-case symbol. An arrow is placed over the upper-case symbol if the vector is physical, while a chevron is placed above the upper-case symbol if the vector is algebraic. Each matrix is denoted by a bold-face upper-case symbol and the elements of that matrix are denoted by the corresponding lower-case symbol.

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}$}$$ (5)

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}}}$$ (6)

where |$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{A}$}$| denotes the magnitude (standard Euclidean norm) of $\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{A}$}$ . Equation (6) 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}}}$$ = (7)

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}$} {\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)}}}$$ (8)

It is convenient for our purposes to define an algebraic 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. \left.\vphantom{ f^{L}_{i,j,k}, \, f^{B}_{i,j,k}, \, f^{D}_{i,j,k} }\right)^{t}_{}$$ (9)

where a superscript ``t'' denotes ``transpose.'' The three face-area components associated with the Right-Top-Up vertex are illustrated in Fig. 5.

Figure 5: Three face-center face-area components defining the flux vector at vertex RTU.

The other vertex vectors are defined in analogy with Eqs. (8) and (9).

As explained in Reference [8], the adjoint relationship between the gradient and divergence operators is embodied in the following integral identity:

$$\oint_{{\partial V}}^{} \phi \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{H}$} \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{n}$} \int_{V}^{} \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{H}$} \mbox{$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{\nabla}$}$} \phi \int_{V}^{} \phi \mbox{$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{\nabla}$} \cdot$} \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{H}$}$$ (10)

where $\phi$ is an arbitrary scalar function, $\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{H}$}$ is an arbitrary vector function, V denotes a volume, $\partial$V denotes its surface, and $\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{n}$}$ denotes the outward-directed unit normal associated with that surface. The vector $\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{H}$}$ has the same mesh locations as the flux vector $\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$}$ , but is not necessarily equal to - D$\mbox{$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{\nabla}$}$}$$\phi$ . We stress that the function $\phi$ at this point represents an arbitrary scalar function, and not necessarily the solution of the diffusion equation. The next step in our support-operators method is to discretize Eq. (10) over a single arbitrary cell in a special manner. Specifically, we explicitly discretize all but the flux operator, which is expressed in an implicit form consistent with our choice of discrete vector unknowns. We assume indices of i, j, k for the arbitrary cell, but suppress these indices whenever possible in the discrete approximation to Eq. (10) that follows. We first discretize the surface integral:

$$\oint_{{\partial V}}^{} \phi \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{H}$} \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{n}$} \approx \phi^{{L}}_{} \phi^{{R}}_{} \phi^{{B}}_{} \phi^{{T}}_{} \phi^{{D}}_{} \phi^{{U}}_{}$$ (11)

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$$

$$\left(\vphantom{ \mbox{$\stackrel{^{\mathstrut}\smash{\longrighta... ...cdot \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$}^{LBD} }\right. \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{H}$} \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$} \left.\vphantom{ \mbox{$\stackrel{^{\mathstrut}\smash{\longrighta... ...cdot \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$}^{LBD} }\right) \left(\vphantom{ \mbox{$\stackrel{^{\mathstrut}\smash{\longrighta... ...cdot \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$}^{RBD} }\right. \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{H}$} \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$} \left.\vphantom{ \mbox{$\stackrel{^{\mathstrut}\smash{\longrighta... ...cdot \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$}^{RBD} }\right)$$ +

$$\left(\vphantom{ \mbox{$\stackrel{^{\mathstrut}\smash{\longrighta... ...cdot \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$}^{LTD} }\right. \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{H}$} \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$} \left.\vphantom{ \mbox{$\stackrel{^{\mathstrut}\smash{\longrighta... ...cdot \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$}^{LTD} }\right) \left(\vphantom{ \mbox{$\stackrel{^{\mathstrut}\smash{\longrighta... ...cdot \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$}^{RTD} }\right. \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{H}$} \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$} \left.\vphantom{ \mbox{$\stackrel{^{\mathstrut}\smash{\longrighta... ...cdot \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$}^{RTD} }\right)$$ +

$$\left(\vphantom{ \mbox{$\stackrel{^{\mathstrut}\smash{\longrighta... ...cdot \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$}^{LBU} }\right. \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{H}$} \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$} \left.\vphantom{ \mbox{$\stackrel{^{\mathstrut}\smash{\longrighta... ...cdot \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$}^{LBU} }\right) \left(\vphantom{ \mbox{$\stackrel{^{\mathstrut}\smash{\longrighta... ...cdot \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$}^{RBU} }\right. \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{H}$} \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$} \left.\vphantom{ \mbox{$\stackrel{^{\mathstrut}\smash{\longrighta... ...cdot \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$}^{RBU} }\right)$$ +

$$\left(\vphantom{ \mbox{$\stackrel{^{\mathstrut}\smash{\longrighta... ...cdot \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$}^{LTU} }\right. \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{H}$} \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$} \left.\vphantom{ \mbox{$\stackrel{^{\mathstrut}\smash{\longrighta... ...cdot \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$}^{LTU} }\right) \left(\vphantom{ \mbox{$\stackrel{^{\mathstrut}\smash{\longrighta... ...cdot \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$}^{RTU} }\right. \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{H}$} \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$} \left.\vphantom{ \mbox{$\stackrel{^{\mathstrut}\smash{\longrighta... ...cdot \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$}^{RTU} }\right)$$ + (12)

where $\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$}$^LBD denotes - D$\mbox{$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{\nabla}$}$}$$\phi$ at the Left-Bottom-Down vertex, and V^LBD denotes the volumetric weight associated with the Left-Bottom-Down vertex. The remaining flux vectors and vertex volumetric weights are analogously indexed. The choice of weights is one of the many free parameters in the support-operators method. We have investigated several different choices. Specifically:

1. Each vertex weight can be given by one-eighth the triple product associated with the vertex. For instance, using the local vertex indexing shown in Fig. 2, the volumetric weight for the Left-Bottom-Down vertex is given by

   $${\frac{{1}}{{8}}} \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{R}$} \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{R}$} \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{R}$}$$ (13)

   where $\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{R}$}$_{i, j} denotes the vector from vertex i to vertex j . Note that these vertex weights do not sum to the total volume of the hexahedron unless the hexahedron is a parallelepiped. We refer to these weights as the triple-product weights.

2. The weights given in Eq. (13) can be normalized, i.e., multiplied by a single constant, so that they sum to the exact cell volume. We refer to these weights as the normalized triple product weights.

3. Each vertex weight can be set equal to the volume of an associated sub-hexahedron. The sub-hexahedra are obtained by using four straight lines to connect each face center with the four edge centers adjacent to it, and by using six straight lines to connect the cell-center with the six face centers. A sub-hexahedron is illustrated in Fig. 6.

   Figure 6: Sub-hexahedron associated with vertex.

   
   Although it may not be obvious, each outer face of each sub-hexahedron coincides with a face of the hexahedron. Thus the volumes of the sub-hexahedra always sum to the total hexahedron volume. This is perhaps the most natural choice for the volumetric weights. We refer to these weights as the sub-hexahedron weights.

4. Each vertex weight can be set to one-eighth of the total hexahedron volume. We refer to these weights as the one-eighth weights.

Computational testing indicates that the sub-hexahedron and one-eighth weights are decidedly inferior to the triple-product and normalized triple-product weights. In particular, the triple-product and normalized triple-product weights both yield a second-order-accurate diffusion discretization, whereas the sub-hexahedron and one-eighth weights yield a first-order accurate diffusion discretization. Although they both give second-order accuracy, the normalized triple-product weights seem to be slightly more accurate than the triple product weights. Thus we use the normalized triple-product weights.

One can evaluate the dot products in Eq. (12) using Eq. (8), but we find it better for our purposes to evaluate them with the algebraic face-area flux vectors defined by Eq. (9). This is achieved by first transforming the face-area vectors to Cartesian vectors and then taking the dot product. Rather than explicitly define the matrix that transforms face-area vectors to Cartesian vectors, we explicitly define its inverse. The desired transformation matrix can then be obtained by either algebraic or numerical inversion. For instance, let us consider the Left-Bottom-Down vertex vectors. We denote the matrix that transforms face-area vectors to Cartesian vectors as A^LBD . Its inverse is the matrix that transforms Cartesian vectors to face-area vectors:

$$\hat{{H}}^{{LBD}}_{} \left[\vphantom{ {\bf A}^{LBD} }\right. \bf A^{{LBD}}_{} \left.\vphantom{ {\bf A}^{LBD} }\right]^{{-1}}_{} \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{H}$}$$ (14)

where $\hat{{H}}$ denotes a Left-Bottom-Down face-area flux vector,

$$\hat{{H}} \left(\vphantom{ h^L, h^B, h^D }\right. \left.\vphantom{ h^L, h^B, h^D }\right)^{t}_{}$$ (15)

and $\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. \left.\vphantom{ h^x, h^y, h^z }\right)^{t}_{}$$ (16)

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]$$ (17)

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. Transforming the face-area vector for the Left-Bottom-Down vertex, we obtain:

$$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{H}$} \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$} \hat{{H}}^{{LBD}}_{} \hat{{F}}^{{LBD}}_{}$$ =

$$\hat{{H}}^{{LBD}}_{} \hat{{F}}^{{LBD}}_{}$$ = (18)

where

$$\bf S^{{LBD}}_{} \left[\vphantom{ {\bf A}^{LBD} }\right. \bf A^{{LBD}}_{} \left.\vphantom{ {\bf A}^{LBD} }\right]^{t}_{} \bf A^{{LBD}}_{}$$ (19)

Following Eq. (19), We now rewrite Eq. (12) in terms of face-area vectors as follows:

$$\int_{V}^{}$$ - D^-1$$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{H}$}$$ ^. D$$\mbox{$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{\nabla}$}$}$$$$\phi$$  dV $$\approx$$

$$\left(\vphantom{ \hat{H}^{LBD} \cdot \mathbf{S}^{LBD} \hat{F}^{LBD} }\right. \hat{{H}}^{{LBD}}_{} \hat{{F}}^{{LBD}}_{} \left.\vphantom{ \hat{H}^{LBD} \cdot \mathbf{S}^{LBD} \hat{F}^{LBD} }\right) \left(\vphantom{ \hat{H}^{RBD} \cdot \mathbf{S}^{RBD} \hat{F}^{RBD} }\right. \hat{{H}}^{{RBD}}_{} \hat{{F}}^{{RBD}}_{} \left.\vphantom{ \hat{H}^{RBD} \cdot \mathbf{S}^{RBD} \hat{F}^{RBD} }\right)$$ +

$$\left(\vphantom{ \hat{H}^{LTD} \cdot \mathbf{S}^{LTD} \hat{F}^{LTD} }\right. \hat{{H}}^{{LTD}}_{} \hat{{F}}^{{LTD}}_{} \left.\vphantom{ \hat{H}^{LTD} \cdot \mathbf{S}^{LTD} \hat{F}^{LTD} }\right) \left(\vphantom{ \hat{H}^{RTD} \cdot \mathbf{S}^{RTD} \hat{F}^{RTD} }\right. \hat{{H}}^{{RTD}}_{} \hat{{F}}^{{RTD}}_{} \left.\vphantom{ \hat{H}^{RTD} \cdot \mathbf{S}^{RTD} \hat{F}^{RTD} }\right)$$ +

$$\left(\vphantom{ \hat{H}^{LBU} \cdot \mathbf{S}^{LBU} \hat{F}^{LBU} }\right. \hat{{H}}^{{LBU}}_{} \hat{{F}}^{{LBU}}_{} \left.\vphantom{ \hat{H}^{LBU} \cdot \mathbf{S}^{LBU} \hat{F}^{LBU} }\right) \left(\vphantom{ \hat{H}^{RBU} \cdot \mathbf{S}^{RBU} \hat{F}^{RBU} }\right. \hat{{H}}^{{RBU}}_{} \hat{{F}}^{{RBU}}_{} \left.\vphantom{ \hat{H}^{RBU} \cdot \mathbf{S}^{RBU} \hat{F}^{RBU} }\right)$$ +

$$\left(\vphantom{ \hat{H}^{LTU} \cdot \mathbf{S}^{LTU} \hat{F}^{LTU} }\right. \hat{{H}}^{{LTU}}_{} \hat{{F}}^{{LTU}}_{} \left.\vphantom{ \hat{H}^{LTU} \cdot \mathbf{S}^{LTU} \hat{F}^{LTU} }\right) \left(\vphantom{ \hat{H}^{RTU} \cdot \mathbf{S}^{RTU} \hat{F}^{RTU} }\right. \hat{{H}}^{{RTU}}_{} \hat{{F}}^{{RTU}}_{} \left.\vphantom{ \hat{H}^{RTU} \cdot \mathbf{S}^{RTU} \hat{F}^{RTU} }\right)$$ + (20)

Although we assume a single diffusion coefficient in each cell in this paper, we note that our scheme can accommodate a different diffusion coefficient for each vertex. In particular, Eq. (20) becomes

$$\int_{V}^{}$$ - D^-1$$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{H}$}$$ ^. D$$\mbox{$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{\nabla}$}$}$$$$\phi$$  dV $$\approx$$

$$\left(\vphantom{ \hat{H}^{LBD} \cdot \mathbf{S}^{LBD} \hat{F}^{LBD} }\right. \hat{{H}}^{{LBD}}_{} \hat{{F}}^{{LBD}}_{} \left.\vphantom{ \hat{H}^{LBD} \cdot \mathbf{S}^{LBD} \hat{F}^{LBD} }\right) \left(\vphantom{ \hat{H}^{RBD} \cdot \mathbf{S}^{RBD} \hat{F}^{RBD} }\right. \hat{{H}}^{{RBD}}_{} \hat{{F}}^{{RBD}}_{} \left.\vphantom{ \hat{H}^{RBD} \cdot \mathbf{S}^{RBD} \hat{F}^{RBD} }\right)$$ +

$$\left(\vphantom{ \hat{H}^{LTD} \cdot \mathbf{S}^{LTD} \hat{F}^{LTD} }\right. \hat{{H}}^{{LTD}}_{} \hat{{F}}^{{LTD}}_{} \left.\vphantom{ \hat{H}^{LTD} \cdot \mathbf{S}^{LTD} \hat{F}^{LTD} }\right) \left(\vphantom{ \hat{H}^{RTD} \cdot \mathbf{S}^{RTD} \hat{F}^{RTD} }\right. \hat{{H}}^{{RTD}}_{} \hat{{F}}^{{RTD}}_{} \left.\vphantom{ \hat{H}^{RTD} \cdot \mathbf{S}^{RTD} \hat{F}^{RTD} }\right)$$ +

$$\left(\vphantom{ \hat{H}^{LBU} \cdot \mathbf{S}^{LBU} \hat{F}^{LBU} }\right. \hat{{H}}^{{LBU}}_{} \hat{{F}}^{{LBU}}_{} \left.\vphantom{ \hat{H}^{LBU} \cdot \mathbf{S}^{LBU} \hat{F}^{LBU} }\right) \left(\vphantom{ \hat{H}^{RBU} \cdot \mathbf{S}^{RBU} \hat{F}^{RBU} }\right. \hat{{H}}^{{RBU}}_{} \hat{{F}}^{{RBU}}_{} \left.\vphantom{ \hat{H}^{RBU} \cdot \mathbf{S}^{RBU} \hat{F}^{RBU} }\right)$$ +

$$\left(\vphantom{ \hat{H}^{LTU} \cdot \mathbf{S}^{LTU} \hat{F}^{LTU} }\right. \hat{{H}}^{{LTU}}_{} \hat{{F}}^{{LTU}}_{} \left.\vphantom{ \hat{H}^{LTU} \cdot \mathbf{S}^{LTU} \hat{F}^{LTU} }\right) \left(\vphantom{ \hat{H}^{RTU} \cdot \mathbf{S}^{RTU} \hat{F}^{RTU} }\right. \hat{{H}}^{{RTU}}_{} \hat{{F}}^{{RTU}}_{} \left.\vphantom{ \hat{H}^{RTU} \cdot \mathbf{S}^{RTU} \hat{F}^{RTU} }\right)$$ + (21)

Although we assume a scalar diffusion coefficient in this paper, we note that our scheme can accommodate a tensor diffusion coefficient. Specifically, with a tensor diffusion coefficient at each vertex, Eq. (21) becomes

$$\int_{V}^{}$$ - D^-1$$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{H}$}$$ ^. D$$\mbox{$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{\nabla}$}$}$$$$\phi$$  dV $$\approx$$

$$\left(\vphantom{ \hat{H}^{LBD} \cdot \mathbf{G}^{LBD} \hat{F}^{LBD} }\right. \hat{{H}}^{{LBD}}_{} \hat{{F}}^{{LBD}}_{} \left.\vphantom{ \hat{H}^{LBD} \cdot \mathbf{G}^{LBD} \hat{F}^{LBD} }\right) \left(\vphantom{ \hat{H}^{RBD} \cdot \mathbf{G}^{RBD} \hat{F}^{RBD} }\right. \hat{{H}}^{{RBD}}_{} \hat{{F}}^{{RBD}}_{} \left.\vphantom{ \hat{H}^{RBD} \cdot \mathbf{G}^{RBD} \hat{F}^{RBD} }\right)$$ +

$$\left(\vphantom{ \hat{H}^{LTD} \cdot \mathbf{G}^{LTD} \hat{F}^{LTD} }\right. \hat{{H}}^{{LTD}}_{} \hat{{F}}^{{LTD}}_{} \left.\vphantom{ \hat{H}^{LTD} \cdot \mathbf{G}^{LTD} \hat{F}^{LTD} }\right) \left(\vphantom{ \hat{H}^{RTD} \cdot \mathbf{G}^{RTD} \hat{F}^{RTD} }\right. \hat{{H}}^{{RTD}}_{} \hat{{F}}^{{RTD}}_{} \left.\vphantom{ \hat{H}^{RTD} \cdot \mathbf{G}^{RTD} \hat{F}^{RTD} }\right)$$ +

$$\left(\vphantom{ \hat{H}^{LBU} \cdot \mathbf{G}^{LBU} \hat{F}^{LBU} }\right. \hat{{H}}^{{LBU}}_{} \hat{{F}}^{{LBU}}_{} \left.\vphantom{ \hat{H}^{LBU} \cdot \mathbf{G}^{LBU} \hat{F}^{LBU} }\right) \left(\vphantom{ \hat{H}^{RBU} \cdot \mathbf{G}^{RBU} \hat{F}^{RBU} }\right. \hat{{H}}^{{RBU}}_{} \hat{{F}}^{{RBU}}_{} \left.\vphantom{ \hat{H}^{RBU} \cdot \mathbf{G}^{RBU} \hat{F}^{RBU} }\right)$$ +

$$\left(\vphantom{ \hat{H}^{LTU} \cdot \mathbf{G}^{LTU} \hat{F}^{LTU} }\right. \hat{{H}}^{{LTU}}_{} \hat{{F}}^{{LTU}}_{} \left.\vphantom{ \hat{H}^{LTU} \cdot \mathbf{G}^{LTU} \hat{F}^{LTU} }\right) \left(\vphantom{ \hat{H}^{RTU} \cdot \mathbf{G}^{RTU} \hat{F}^{RTU} }\right. \hat{{H}}^{{RTU}}_{} \hat{{F}}^{{RTU}}_{} \left.\vphantom{ \hat{H}^{RTU} \cdot \mathbf{G}^{RTU} \hat{F}^{RTU} }\right)$$ + (22)

where

$$\left[\vphantom{ {\bf A}^{LBD} }\right. \bf A^{{LBD}}_{} \left.\vphantom{ {\bf A}^{LBD} }\right]^{t}_{} \left[\vphantom{ \mathbf{D}^{LBD} }\right. \left.\vphantom{ \mathbf{D}^{LBD} }\right]^{{-1}}_{} \bf A^{{LBD}}_{}$$ (23)

and D^LBD is the Left-Bottom-Down diffusion tensor in the Cartesian basis. The remaining G -matrices are defined analogously. The diffusion tensor must be symmetric positive-definite to ensure that its inverse exists and that the coefficient matrix for our diffusion scheme is symmetric positive-definite.

Finally, we approximate the divergence volumetric integral:

$$\int_{V}^{} \phi \mbox{$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{\nabla}$} \cdot$} \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{H}$} \approx \phi^{{C}}_{} \left[\vphantom{ h^{L} + h^{R} + h^{B} + h^{T} + h^{D} + h^{U} }\right. \left.\vphantom{ h^{L} + h^{R} + h^{B} + h^{T} + h^{D} + h^{U} }\right]$$ (24)

Equations (11), (20), and (24) are certainly not unique, but they are fairly straightforward. For instance, Eq. (11) represents a face-centered second-order approximation to a surface integral. Equation (20) represents a vertex-based volumetric integral consisting of a dot-product contribution from each pair of vertex vectors. Equation (24) is a particularly simple second-order approximation which gives all of the weight to the cell-center value of $\phi$ while using a surface-integral formulation for $\mbox{$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{\nabla}$} \cdot$}$$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{H}$}$ that is analogous to the surface-integral used in Eq. (11).

Substituting from Eqs. (11), (20), and (24) into Eq. (10), we obtain the discrete version of Eq. (10):

$$\phi^{{L}}_{}$$hL + $$\phi^{{R}}_{}$$hR + $$\phi^{{B}}_{}$$hB + $$\phi^{{T}}_{}$$hT + $$\phi^{{D}}_{}$$hD + $$\phi^{{U}}_{}$$hU + D-1$$\left(\vphantom{ \hat{H}^{LBD} \cdot \mathbf{S}^{LBD} \hat{F}^{LBD} }\right.$$$$\hat{{H}}^{{LBD}}_{}$$ . SLBD$$\hat{{F}}^{{LBD}}_{}$$$$\left.\vphantom{ \hat{H}^{LBD} \cdot \mathbf{S}^{LBD} \hat{F}^{LBD} }\right)$$VLBD + D-1$$\left(\vphantom{ \hat{H}^{RBD} \cdot \mathbf{S}^{RBD} \hat{F}^{RBD} }\right.$$$$\hat{{H}}^{{RBD}}_{}$$ . SRBD$$\hat{{F}}^{{RBD}}_{}$$$$\left.\vphantom{ \hat{H}^{RBD} \cdot \mathbf{S}^{RBD} \hat{F}^{RBD} }\right)$$VRBD + D-1$$\left(\vphantom{ \hat{H}^{LTD} \cdot \mathbf{S}^{LTD} \hat{F}^{LTD} }\right.$$$$\hat{{H}}^{{LTD}}_{}$$ . SLTD$$\hat{{F}}^{{LTD}}_{}$$$$\left.\vphantom{ \hat{H}^{LTD} \cdot \mathbf{S}^{LTD} \hat{F}^{LTD} }\right)$$VLTD + D-1$$\left(\vphantom{ \hat{H}^{RTD} \cdot \mathbf{S}^{RTD} \hat{F}^{RTD} }\right.$$$$\hat{{H}}^{{RTD}}_{}$$ . SRTD$$\hat{{F}}^{{RTD}}_{}$$$$\left.\vphantom{ \hat{H}^{RTD} \cdot \mathbf{S}^{RTD} \hat{F}^{RTD} }\right)$$VRTD + D-1$$\left(\vphantom{ \hat{H}^{LBU} \cdot \mathbf{S}^{LBU} \hat{F}^{LBU} }\right.$$$$\hat{{H}}^{{LBU}}_{}$$ . SLBU$$\hat{{F}}^{{LBU}}_{}$$$$\left.\vphantom{ \hat{H}^{LBU} \cdot \mathbf{S}^{LBU} \hat{F}^{LBU} }\right)$$VLBU + D-1$$\left(\vphantom{ \hat{H}^{RBU} \cdot \mathbf{S}^{RBU} \hat{F}^{RBU} }\right.$$$$\hat{{H}}^{{RBU}}_{}$$ . SRBU$$\hat{{F}}^{{RBU}}_{}$$$$\left.\vphantom{ \hat{H}^{RBU} \cdot \mathbf{S}^{RBU} \hat{F}^{RBU} }\right)$$VRBU + D-1$$\left(\vphantom{ \hat{H}^{LTU} \cdot \mathbf{S}^{LTU} \hat{F}^{LTU} }\right.$$$$\hat{{H}}^{{LTU}}_{}$$ . SLTU$$\hat{{F}}^{{LTU}}_{}$$$$\left.\vphantom{ \hat{H}^{LTU} \cdot \mathbf{S}^{LTU} \hat{F}^{LTU} }\right)$$VLTU + D-1$$\left(\vphantom{ \hat{H}^{RTU} \cdot \mathbf{S}^{LTU} \hat{F}^{RTU} }\right.$$$$\hat{{H}}^{{RTU}}_{}$$ . SLTU$$\hat{{F}}^{{RTU}}_{}$$$$\left.\vphantom{ \hat{H}^{RTU} \cdot \mathbf{S}^{LTU} \hat{F}^{RTU} }\right)$$VRTU = $$\phi^{{C}}_{}$$$$\left[\vphantom{ h^{L} + h^{R} + h^{B} + h^{T} + h^{D} + h^{U} }\right.$$hL + hR + hB + hT + hD + hU$$\left.\vphantom{ h^{L} + h^{R} + h^{B} + h^{T} + h^{D} + h^{U} }\right]$$   .

(25)

Note that Eq. (25) defines the discrete inner products, discussed in Reference 8, that are associated with the adjoint relationship between the divergence and gradient operators. We can now use this relationship to solve for the flux operator components by requiring that the resulting discretized identity hold for all discrete $\hat{{H}}$ and $\phi$ values. In particular, the equation for the face-area component of $\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$}$ on any given cell face is obtained from Eq. (25) simply by setting the same face-area component of $\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{H}$}$ on that face to unity and setting the remaining face-area components of $\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{H}$}$ on all other faces to zero. For instance, we obtain the equation for f^L from Eq. (25) by setting h^L to unity and all the other face-area components of $\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{H}$}$ , i.e., h^R , h^B , h^T , h^D , h^U , to zero:

$$\phi^{L}_{}$$ +	D-1$$\left(\vphantom{ s_{L,L}^{LBD} f^L + s_{L,B}^{LBD} f^B + s_{L,D}^{LBD} f^D }\right.$$sL, LLBDfL + sL, BLBDfB + sL, DLBDfD$$\left.\vphantom{ s_{L,L}^{LBD} f^L + s_{L,B}^{LBD} f^B + s_{L,D}^{LBD} f^D }\right)$$VLBD
+	D-1$$\left(\vphantom{ s_{L,L}^{LTD} f^L + s_{L,T}^{LTD} f^T + s_{L,D}^{LTD} f^D }\right.$$sL, LLTDfL + sL, TLTDfT + sL, DLTDfD$$\left.\vphantom{ s_{L,L}^{LTD} f^L + s_{L,T}^{LTD} f^T + s_{L,D}^{LTD} f^D }\right)$$VLTD
+	D-1$$\left(\vphantom{ s_{L,L}^{LBU} f^L + s_{L,B}^{LBU} f^B + s_{L,U}^{LBU} f^U }\right.$$sL, LLBUfL + sL, BLBUfB + sL, ULBUfU$$\left.\vphantom{ s_{L,L}^{LBU} f^L + s_{L,B}^{LBU} f^B + s_{L,U}^{LBU} f^U }\right)$$VLBU
+	D-1$$\left(\vphantom{ s_{L,L}^{LTU} f^L + s_{L,T}^{LTU} f^T + s_{L,U}^{LTU} f^U }\right.$$sL, LLTUfL + sL, TLTUfT + sL, ULTUfU$$\left.\vphantom{ s_{L,L}^{LTU} f^L + s_{L,T}^{LTU} f^T + s_{L,U}^{LTU} f^U }\right)$$VLTU	= $$\phi^{C}_{}$$   ,

(26)

where s_{L, L}^LBD denotes the (L, L) element of the matrix $\mathcal {S}$^LBD defined by Eq. (19). The remaining S -matrix elements are defined analogously. We obtain the equation for f^R from Eq. (25) by setting h^R to unity and all the other face components of $\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{H}$}$ to zero:

$$\phi^{R}_{}$$ +	D-1$$\left(\vphantom{ s_{R,R}^{RBD} f^R + s_{R,B}^{RBD} f^B + s_{R,D}^{RBD} f^D }\right.$$sR, RRBDfR + sR, BRBDfB + sR, DRBDfD$$\left.\vphantom{ s_{R,R}^{RBD} f^R + s_{R,B}^{RBD} f^B + s_{R,D}^{RBD} f^D }\right)$$VRBD
+	D-1$$\left(\vphantom{ s_{R,R}^{RTD} f^R + s_{R,T}^{RTD} f^T + s_{R,D}^{RTD} f^D }\right.$$sR, RRTDfR + sR, TRTDfT + sR, DRTDfD$$\left.\vphantom{ s_{R,R}^{RTD} f^R + s_{R,T}^{RTD} f^T + s_{R,D}^{RTD} f^D }\right)$$VRTD
+	D-1$$\left(\vphantom{ s_{R,R}^{RBU} f^R + s_{R,B}^{RBU} f^B + s_{R,U}^{RBU} f^U }\right.$$sR, RRBUfR + sR, BRBUfB + sR, URBUfU$$\left.\vphantom{ s_{R,R}^{RBU} f^R + s_{R,B}^{RBU} f^B + s_{R,U}^{RBU} f^U }\right)$$VRBU
+	D-1$$\left(\vphantom{ s_{R,R}^{RTU} f^R + s_{R,T}^{RTU} f^T + s_{R,U}^{RTU} f^U }\right.$$sR, RRTUfR + sR, TRTUfT + sR, URTUfU$$\left.\vphantom{ s_{R,R}^{RTU} f^R + s_{R,T}^{RTU} f^T + s_{R,U}^{RTU} f^U }\right)$$VRTU	= $$\phi^{C}_{}$$   .

(27)

We obtain the equation for f^B from Eq. (25) by setting h^B to unity and all the other face components of $\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{H}$}$ to zero:

$$\phi^{B}_{}$$ +	D-1$$\left(\vphantom{ s_{B,L}^{LBD} f^L + s_{B,B}^{LBD} f^B + s_{B,D}^{LBD} f^D }\right.$$sB, LLBDfL + sB, BLBDfB + sB, DLBDfD$$\left.\vphantom{ s_{B,L}^{LBD} f^L + s_{B,B}^{LBD} f^B + s_{B,D}^{LBD} f^D }\right)$$VLBD
+	D-1$$\left(\vphantom{ s_{B,R}^{RBD} f^R + s_{B,B}^{RBD} f^B + s_{B,D}^{RBD} f^D }\right.$$sB, RRBDfR + sB, BRBDfB + sB, DRBDfD$$\left.\vphantom{ s_{B,R}^{RBD} f^R + s_{B,B}^{RBD} f^B + s_{B,D}^{RBD} f^D }\right)$$VRBD
+	D-1$$\left(\vphantom{ s_{B,L}^{LBU} f^L + s_{B,B}^{LBU} f^B + s_{B,U}^{LBU} f^U }\right.$$sB, LLBUfL + sB, BLBUfB + sB, ULBUfU$$\left.\vphantom{ s_{B,L}^{LBU} f^L + s_{B,B}^{LBU} f^B + s_{B,U}^{LBU} f^U }\right)$$VLBU
+	D-1$$\left(\vphantom{ s_{B,R}^{RBU} f^R + s_{B,B}^{RBU} f^B + s_{B,U}^{RBU} f^U }\right.$$sB, RRBUfR + sB, BRBUfB + sB, URBUfU$$\left.\vphantom{ s_{B,R}^{RBU} f^R + s_{B,B}^{RBU} f^B + s_{B,U}^{RBU} f^U }\right)$$VRBU	= $$\phi^{C}_{}$$   .

(28)

We obtain the equation for f^T from Eq. (25) by setting h^T to unity and all the other face components of $\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{H}$}$ to zero:

$$\phi^{T}_{}$$ +	D-1$$\left(\vphantom{ s_{T,L}^{LTD} f^L + s_{T,T}^{LTD} f^T + s_{T,D}^{LTD} f^D }\right.$$sT, LLTDfL + sT, TLTDfT + sT, DLTDfD$$\left.\vphantom{ s_{T,L}^{LTD} f^L + s_{T,T}^{LTD} f^T + s_{T,D}^{LTD} f^D }\right)$$VLTD
+	D-1$$\left(\vphantom{ s_{T,R}^{RTD} f^R + s_{T,T}^{RTD} f^T + s_{T,D}^{RTD} f^D }\right.$$sT, RRTDfR + sT, TRTDfT + sT, DRTDfD$$\left.\vphantom{ s_{T,R}^{RTD} f^R + s_{T,T}^{RTD} f^T + s_{T,D}^{RTD} f^D }\right)$$VRTD
+	D-1$$\left(\vphantom{ s_{T,L}^{LTU} f^L + s_{T,T}^{LTU} f^T + s_{T,U}^{LTU} f^U }\right.$$sT, LLTUfL + sT, TLTUfT + sT, ULTUfU$$\left.\vphantom{ s_{T,L}^{LTU} f^L + s_{T,T}^{LTU} f^T + s_{T,U}^{LTU} f^U }\right)$$VLTU
+	D-1$$\left(\vphantom{ s_{T,R}^{RTU} f^R + s_{T,T}^{RTU} f^T + s_{T,U}^{RTU} f^U }\right.$$sT, RRTUfR + sT, TRTUfT + sT, URTUfU$$\left.\vphantom{ s_{T,R}^{RTU} f^R + s_{T,T}^{RTU} f^T + s_{T,U}^{RTU} f^U }\right)$$VRTU	= $$\phi^{C}_{}$$   .

(29)

We obtain the equation for f^D from Eq. (25) by setting h^D to unity and all the other face components of $\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{H}$}$ to zero:

$$\phi^{D}_{}$$ +	D-1$$\left(\vphantom{ s_{D,L}^{LBD} f^L + s_{D,B}^{LBD} f^B + s_{D,D}^{LBD} f^D }\right.$$sD, LLBDfL + sD, BLBDfB + sD, DLBDfD$$\left.\vphantom{ s_{D,L}^{LBD} f^L + s_{D,B}^{LBD} f^B + s_{D,D}^{LBD} f^D }\right)$$VLBD
+	D-1$$\left(\vphantom{ s_{D,R}^{RBD} f^R + s_{D,B}^{RBD} f^B + s_{D,D}^{RBD} f^D }\right.$$sD, RRBDfR + sD, BRBDfB + sD, DRBDfD$$\left.\vphantom{ s_{D,R}^{RBD} f^R + s_{D,B}^{RBD} f^B + s_{D,D}^{RBD} f^D }\right)$$VRBD
+	D-1$$\left(\vphantom{ s_{D,L}^{LTD} f^L + s_{D,T}^{LTD} f^T + s_{D,D}^{LTD} f^D }\right.$$sD, LLTDfL + sD, TLTDfT + sD, DLTDfD$$\left.\vphantom{ s_{D,L}^{LTD} f^L + s_{D,T}^{LTD} f^T + s_{D,D}^{LTD} f^D }\right)$$VLTD
+	D-1$$\left(\vphantom{ s_{D,R}^{RTD} f^R + s_{D,T}^{RTD} f^T + s_{D,D}^{RTD} f^D }\right.$$sD, RRTDfR + sD, TRTDfT + sD, DRTDfD$$\left.\vphantom{ s_{D,R}^{RTD} f^R + s_{D,T}^{RTD} f^T + s_{D,D}^{RTD} f^D }\right)$$VRTD	= $$\phi^{C}_{}$$   .

(30)

Finally, we obtain the equation for f^U from Eq. (25) by setting h^U to unity and all the other face components of $\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{H}$}$ to zero:

$$\phi^{U}_{}$$ +	D-1$$\left(\vphantom{ s_{U,L}^{LBU} f^L + s_{U,B}^{LBU} f^B + s_{U,U}^{LBU} f^U }\right.$$sU, LLBUfL + sU, BLBUfB + sU, ULBUfU$$\left.\vphantom{ s_{U,L}^{LBU} f^L + s_{U,B}^{LBU} f^B + s_{U,U}^{LBU} f^U }\right)$$VLBU
+	D-1$$\left(\vphantom{ s_{U,R}^{RBU} f^R + s_{U,B}^{RBU} f^B + s_{U,U}^{RBU} f^U }\right.$$sU, RRBUfR + sU, BRBUfB + sU, URBUfU$$\left.\vphantom{ s_{U,R}^{RBU} f^R + s_{U,B}^{RBU} f^B + s_{U,U}^{RBU} f^U }\right)$$VRBU
+	D-1$$\left(\vphantom{ s_{U,L}^{LTU} f^L + s_{U,T}^{LTU} f^T + s_{U,U}^{LTU} f^U }\right.$$sU, LLTUfL + sU, TLTUfT + sU, ULTUfU$$\left.\vphantom{ s_{U,L}^{LTU} f^L + s_{U,T}^{LTU} f^T + s_{U,U}^{LTU} f^U }\right)$$VLTU
+	D-1$$\left(\vphantom{ s_{U,R}^{RTU} f^R + s_{U,T}^{RTU} f^T + s_{U,U}^{RTU} f^U }\right.$$sU, RRTUfR + sU, TRTUfT + sU, URTUfU$$\left.\vphantom{ s_{U,R}^{RTU} f^R + s_{U,T}^{RTU} f^T + s_{U,U}^{RTU} f^U }\right)$$VRTU	= $$\phi^{C}_{}$$   .

(31)

Equations (26) through (31) can be expressed in matrix form as follows:

$$\hat{{\mathcal{F}}} \Delta \hat{{\Phi}}$$ (32)

where

$$\hat{{\mathcal{F}}} \left(\vphantom{ f^L, f^R, f^B, f^T, f^D, f^U }\right. \left.\vphantom{ f^L, f^R, f^B, f^T, f^D, f^U }\right)^{{t}}_{}$$ (33)

and

$$\Delta \hat{{\Phi}} \left(\vphantom{ \phi^C -\phi^L, \phi^C - \phi^R, \phi^C - \phi^B, \phi^C - \phi^T, \phi^C - \phi^D, \phi^C - \phi^U }\right. \phi^{C}_{} \phi^{L}_{} \phi^{C}_{} \phi^{R}_{} \phi^{C}_{} \phi^{B}_{} \phi^{C}_{} \phi^{T}_{} \phi^{C}_{} \phi^{D}_{} \phi^{C}_{} \phi^{U}_{} \left.\vphantom{ \phi^C -\phi^L, \phi^C - \phi^R, \phi^C - \phi^B, \phi^C - \phi^T, \phi^C - \phi^D, \phi^C - \phi^U }\right)^{{t}}_{}$$ (34)

To obtain a matrix that gives the face-center components of the flux operator in terms of the face-center and cell-center intensities, one need simply invert the 6 x 6 matrix in Eq. (32):