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

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

Since it is not practical to perform this inversion algebraically, we perform it numerically. Thus we cannot give an explicit expression for the matrix W . Nonetheless, it can be shown that it is an SPD matrix (see the Appendix.) In addition, if we assume a rectangular mesh, W becomes diagonal and can be trivially inverted. For instance, under this assumption, Eq. (26) becomes:

$$\phi^{L}_{} \left(\vphantom{ \Delta y \Delta z }\right. \Delta \Delta \left.\vphantom{ \Delta y \Delta z }\right)^{{-2}}_{} {\frac{{\Delta x \Delta y \Delta z}}{{2}}} \phi_{C}^{}$$ (36)

where we have also assumed that the indices i , j , k , correspond to the spatial coordinates x , y , z , respectively. Solving Eq. (36) for f^L , we obtain

$${\frac{{2D}}{{\Delta x}}} \left(\vphantom{ \phi^L - \phi^C }\right. \phi^{L}_{} \phi^{C}_{} \left.\vphantom{ \phi^L - \phi^C }\right) \Delta \Delta$$ (37)

which is exact for $\phi$ linearly-dependent upon x .

Having derived Eq. (35), we can construct the discrete equation for the cell-center intensity in every cell. Each such equation represents a discretization of Eq. (3), i.e., a balance equation for the cell. Furthermore, each balance equation uses a discretization for the divergence of the flux that is identical to that used in Eq. (25). In some sense, this is the point at which we obtain a diffusion operator by combining our discrete divergence and flux operators. Specifically, the equation for $\phi^{C}_{}$ is:

$${\frac{{\partial \phi^C}}{{\partial t}}}$$ (38)

where V denotes the total volume of the cell, the face-area flux components are expressed in terms of the intensities via Eq. (35), and Q^C denotes the source or driving function evaluated at cell-center. We have chosen not to discretize the time derivative in Eq. (38) simply because essentially any standard discretization, e.g., the backward-Euler and Crank-Nicholson schemes [9], can be applied in conjunction with our spatial discretization. Equation (38) contains all of the intensities in cell (i, j, k) . Thus it has a 7-point stencil.

Now that we have defined the equations for the cell-center intensities, we must next define equations for the face-center intensities. Our local indexing scheme admits two intensities and two face-area flux components at each face on the mesh interior. In particular, there is one intensity and one flux component from each of the cells that share a face. For instance, the cell face with global index (i + ${\frac{{1}}{{2}}}$, j, k) is associated with the two intensities, $\phi^{R}_{{i,j,k}}$ and $\phi^{L}_{{i+1,j,k}}$ , and the two face-area flux components, f^R_{i, j, k} and f^L_{i+1, j, k} . We previously obtained the flux components in terms of the intensities by forcing Eq. (25), a discrete version of Eq. (4), to be satisfied on each individual cell for all discrete scalars and vectors. We now obtain equations for the interior-mesh face-center intensities by requiring that this identity be satisfied over the entire mesh for all discrete scalars and vectors.

When Eq. (25) is summed over the entire mesh, the two volumetric integrals are naturally approximated in terms of a sum of contributions from each individual cell. However, a valid approximation for the the surface integral in Eq. (25) will occur if and only if contributions to the surface integral from each individual cell cancel at all interior faces, thereby resulting in an approximate integral over the outer surface of the mesh. By inspection of Eq. (25) it can be seen that this will be achieved by requiring both continuity of the intensity and continuity of the face-area flux component at each interior cell face. In particular, we require that

$$\phi^{R}_{{i,j,k}} \phi^{L}_{{i+1,j,k}} \equiv \phi_{{i+\frac{1}{2},j,k}}^{}$$ (39)

$$\phi^{T}_{{i,j,k}} \phi^{B}_{{i,j+1,k}} \equiv \phi_{{i,j+\frac{1}{2},k}}^{}$$ (40)

$$\phi^{U}_{{i,j,k}} \phi^{D}_{{i,j,k+1}} \equiv \phi_{{i,j,k+\frac{1}{2}}}^{}$$ (41)

f^R_{i, j, k} + f^L_{i+1, j, k} = 0   , (42)

f^T_{i, j, k} + f^B_{i, j+1, k} = 0   , (43)

f^U_{i, j, k} + f^D_{i, j, k+1} = 0   , (44)

where the indices in Eqs. (39) through (44) take on all values associated with interior cell faces, and the flux components in Eqs. (42) through (44) are expressed in terms of intensities via Eq. (35). One would expect that the continuity of the face-area flux components expressed by Eqs. (42) through (44) would require that the difference of the components be zero rather than the sum of the components. However, one must remember that each of the components is defined with respect to an area vector that is equal in magnitude but opposite in direction to that of the other component.

Equations (39) through (41) establish that there is only one intensity unknown associated with each interior-mesh cell face. Thus, as shown in Eqs. (39) through (41), each such intensity can be uniquely referred to using a global mesh index. The equations for these intensities are given by Eqs. (42) through (44). For instance, Eq. (42) is the equation for $\phi_{{i+\frac{1}{2},j,k}}^{}$ . In general, Eq. (42) contains only and all of the intensities in cells (i, j, k) and (i + 1, j, k) . Thus it has a 13-point stencil. The only intensity shared by these two cells is $\phi_{{i+\frac{1}{2},j,k}}^{}$ . Thus in a certain sense it can be said that $\phi_{{i+\frac{1}{2},j,k}}^{}$ is ``chosen'' to obtain continuity of the face-area flux components on cell-face (i + ${\frac{{1}}{{2}}}$, j, k) . The properties of Eqs. (43) and (44) are completely analogous to those of Eq. (42).

If the mesh is orthogonal, Eqs. (42) through (44) simplify to such an extent that they relate each interior-mesh face-center intensity to the two cell-center intensities adjacent to it. This enables the face-center intensities to be explicitly eliminated, resulting in the standard 7-point cell-centered diffusion discretization. This is completely analogous to the 2-D case discussed in detail in [1]. However, if the mesh is non-orthogonal, the face-center intensities cannot be eliminated, and Eqs. (42) through (44) must be included in the diffusion matrix. In this case, these equations must be reversed in sign to obtain a symmetric diffusion matrix:

- f^R_{i, j, k} - f^L_{i+1, j, k} = 0   , (45)

- f^T_{i, j, k} - f^B_{i, j+1, k} = 0   , (46)

- f^U_{i, j, k} - f^D_{i, j, k+1} = 0   . (47)

Having defined the equations for the cell-center and interior-mesh face-center intensities, we need only define the equations for the face-center intensities on the outer mesh boundary to complete the specification of our diffusion discretization scheme. Cell faces on the outer boundary are associated with only one cell. Thus there is only one face-center intensity and one face-area flux component associated with each such face. The equation for each boundary intensity is very similar to that for each interior-mesh face-center intensity in that it expresses a continuity of the face-normal flux component. The only difference in the boundary equations is that the analytic boundary condition for the diffusion equation is used to define a ``ghost-cell'' face-normal flux component that must be equated to the standard face-normal flux component defined by Eq. (35). A ghost cell is a non-existent mesh cell that represents a continuation of the mesh across the outer mesh boundary. For instance, assuming that the left face of cell 1, j, k is on the outer boundary of the mesh and its remaining faces are on the interior of the mesh, the ghost cell ``adjacent'' to cell 1, j, k carries the indices 0, j, k .

The analytic diffusion boundary condition of interest to us is the so-called ``extrapolated'' boundary condition. This condition is of the mixed or Robin type and can be expressed as follows:

$$\phi \mbox{$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{\nabla}$}$} \phi \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{n}$} \phi^{e}_{}$$ (48)

where d^e is called the extrapolation distance, $\phi^{e}_{}$ is called the extrapolated intensity (a specified function), and $\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{n}$}$ denotes an outward-directed unit normal vector. Equation (48) is satisfied at each point on the outer surface of the problem domain. Of course, the values of the parameters, d^e and $\phi^{e}_{}$ , may change as a function of position. One obtains a vacuum boundary condition when $\phi^{e}_{}$ = 0 , a source condition when $\phi^{e}_{}$ is non-zero, and a reflective (Neumann) condition when $\phi^{e}_{}$ = $\phi$ . The extrapolated boundary condition is said to be a Marshak condition whenever d^e = 2D .

We begin the derivation of the ghost-cell face-area flux component by substituting from Eq. (2) into Eq. (48):

$$\phi {\frac{{d^e}}{{D}}} \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$} \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{n}$} \phi^{e}_{}$$ (49)

where $\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$}$^g is the flux vector associated with a ghost cell. Next we recognize that the outward-directed unit normal vector for a ghost-cell must be identical to an inward-directed unit normal vector on the outer surface of the problem domain. Thus

$$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{n}$} \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{n}$}$$ (50)

where $\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{n}$}$^g denotes a ghost-cell outward-directed unit normal vector. Substituting from Eq. (50) into Eq. (49), we obtain:

$$\phi {\frac{{d^e}}{{D}}} \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$} \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{n}$} \phi^{e}_{}$$ (51)

Next we solve Eq. (51) for the outward-directed flux component associated with a ghost cell:

$$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$} \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{n}$} {\frac{{D}}{{d^e}}} \left(\vphantom{ \phi^e - \phi }\right. \phi^{e}_{} \phi \left.\vphantom{ \phi^e - \phi }\right)$$ (52)

Now let us assume that the left face of cell 1, j, k is on the outer boundary of the mesh with its remaining faces on the mesh interior. The ghost cell whose right face is identical to the left face of cell 1, j, k carries the indices 0, j, k . The intensity on the left face of cell (1, j, k) is $\phi_{{\frac{1}{2},j,k}}^{}$ and the face-area flux component on that face is f^L_{1, j, k} . Evaluating Eq. (52) at the center of face (${\frac{{1}}{{2}}}$, j, k) and multiplying the resulting expression by the magnitude of the outward-directed area-vector on that face associated with cell 1, j, k , we obtain the desired expression for the ghost-cell face-area flux component:

$${\frac{{D_{1,j,k}}}{{d^e_{0,j,k}}}} \left(\vphantom{ \phi_{\frac{1}{2},j,k} - \phi^e_{0,j,k} }\right. \phi_{{\frac{1}{2},j,k}}^{} \phi^{e}_{{0,j,k}} \left.\vphantom{ \phi_{\frac{1}{2},j,k} - \phi^e_{0,j,k} }\right) \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{A}$}$$ (53)

where the extrapolated intensity and the extrapolation distance are assumed to carry the ghost-cell index.

We next obtain the equation for $\phi_{{\frac{1}{2},j,k}}^{}$ by requiring that the Right and Left face-area flux components for cells (0, j, k) and (1, j, k) , respectively, sum to zero:

- f^R_{0, j, k} - f^L_{1, j, k} = 0   . (54)

Note that Eq. (54) is identical to Eq. (45) with the latter equation evaluated at i = 0 . Thus Eqs. (45) through (47) provide all face-center intensity equations with the caveat that when an intensity is on the outer mesh boundary, the associated ghost-cell flux component must be defined via the boundary condition rather than Eq. (35). Note that Eq. (54) couples all of the intensities within a cell and therefore has a 7-point stencil. This completes the specification of our diffusion discretization scheme.

To summarize,

• The face-area flux components for each cell are expressed in terms of the intensities within that cell via Eq. (35).
• The discrete equation for each cell-centered intensity is given in Eq. (38).
• The equations for the interior-mesh face-centered intensities are given in Eqs. (45) through (47).
• The equation for a face-center intensity on the outer mesh boundary is given by Eqs. (53) and (54) when the boundary face is the Left face of a cell. Analogous equations for the other five cases are easily derived using Eqs. (45) through (47) and Eq. (53).

It is interesting to note that the equation for a cell-center intensity contains a time-derivative of that intensity, but the equations for the face-center equations do not contain any form of time derivative. Thus in time-dependent calculations, one must have initial values for the cell-center intensities, but initial values are not required for the face-center intensities. Thus only cell-center intensities must be saved from one time step to the next.

We have already shown that our diffusion matrix is sparse. It is also symmetric positive-definite. We demonstrate this latter property in the Appendix.