In this section we describe the supportoperators method. It is convenient at this point to define a flux operator given by  D . The diffusion operator of interest is given by the product of the divergence operator and the flux operator:  D . The supportoperators method is based upon the following three facts:
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 logicallyrectangular mesh. However, our discretization scheme can be used with unstructured meshes as well. The assumption of a logicallyrectangular 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 halfintegral global indices, e.g., (i + , j + , k + ) ; and face centers carry mixed global indices composed of both integral and halfintegral indices, e.g., (i + , j, k) . The global indices for four of the vertices associated with cell (i, j, k) are illustrated in Fig. 1.
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.
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.The vector and matrix notation used from this point forward in this paper is as follows. Each vector is denoted by an uppercase symbol and the components of that vector are denoted by the corresponding lowercase symbol. An arrow is placed over the uppercase symbol if the vector is physical, while a chevron is placed above the uppercase symbol if the vector is algebraic. Each matrix is denoted by a boldface uppercase symbol and the elements of that matrix are denoted by the corresponding lowercase symbol.
The intensities (scalars) are defined to exist at both cell center: , and on the face centers: , , , , , . 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 .
Vectors are defined in terms of facearea 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 with the outwarddirected area vector located at the center of the left face of cell i, j, k . The other facearea components are defined analogously. The area vector is defined as the integral of the outwarddirected unit normal vector over the face, i.e.,
where is a unit vector that is normal to the face at each point on the face. The average outwarddirected unit normal vector for the face is defined as follows: where  denotes the magnitude (standard Euclidean norm) of . Equation (6) can be used to convert facearea flux components to facenormal components if desired, e.g.^{ . }  =  ^{ . } ,  
=  .  (7) 
As explained in Reference [8], the adjoint relationship between the gradient and divergence operators is embodied in the following integral identity:
where is an arbitrary scalar function, is an arbitrary vector function, V denotes a volume, V denotes its surface, and denotes the outwarddirected unit normal associated with that surface. The vector has the same mesh locations as the flux vector , but is not necessarily equal to  D . We stress that the function at this point represents an arbitrary scalar function, and not necessarily the solution of the diffusion equation. The next step in our supportoperators 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:Next we approximate the flux volumetric integral:
D^{1}^{LBD . }^{LBD}V^{LBD}  +  D^{1}^{RBD . }^{RBD}V^{RBD}  
D^{1}^{LTD . }^{LTD}V^{LTD}  +  D^{1}^{RTD . }^{RTD}V^{RTD}  
D^{1}^{LBU . }^{LBU}V^{LBU}  +  D^{1}^{RBU . }^{RBU}V^{RBU}  
D^{1}^{LTU . }^{LTU}V^{LTU}  +  D^{1}^{RTU . }^{RTU}V^{RTU} ,  (12) 
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 facearea flux vectors defined by Eq. (9). This is achieved by first transforming the facearea vectors to Cartesian vectors and then taking the dot product. Rather than explicitly define the matrix that transforms facearea 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 LeftBottomDown vertex vectors. We denote the matrix that transforms facearea vectors to Cartesian vectors as A^{LBD} . Its inverse is the matrix that transforms Cartesian vectors to facearea vectors:
where denotes a LeftBottomDown facearea flux vector, and denotes a LeftBottomDown Cartesian flux vector, and where a^{L}_{x} denotes the xcomponent of the area vector associated with the left face. The remaining components of the matrix are defined analogously. Transforming the facearea vector for the LeftBottomDown vertex, we obtain:
D^{1}^{ . }S^{LBD}V^{LBD}  +  D^{1}^{ . }S^{RBD}V^{RBD}  
D^{1}^{ . }S^{LTD}V^{LTD}  +  D^{1}^{ . }S^{RTD}V^{RTD}  
D^{1}^{ . }S^{LBU}V^{LBU}  +  D^{1}^{ . }S^{RBU}V^{RBU}  
D^{1}^{ . }S^{LTU}V^{LTU}  +  D^{1}^{ . }S^{RTU}V^{RTU} .  (20) 
D^{LBD1}^{ . }S^{LBD}V^{LBD}  +  D^{RBD1}^{ . }S^{RBD}V^{RBD}  
D^{LTD1}^{ . }S^{LTD}V^{LTD}  +  D^{RTD1}^{ . }S^{RTD}V^{RTD}  
D^{LBU1}^{ . }S^{LBU}V^{LBU}  +  D^{RBU1}^{ . }S^{RBU}V^{RBU}  
D^{LTU1}^{ . }S^{LTU}V^{LTU}  +  D^{RTU1}^{ . }S^{RTU}V^{RTU} ,  (21) 
^{ . }G^{LBD}V^{LBD}  +  ^{ . }G^{RBD}V^{RBD}  
^{ . }G^{LTD}V^{LTD}  +  ^{ . }G^{RTD}V^{RTD}  
^{ . }G^{LBU}V^{LBU}  +  ^{ . }G^{RBU}V^{RBU}  
^{ . }G^{LTU}V^{LTU}  +  ^{ . }G^{RTU}V^{RTU} ,  (22) 
Finally, we approximate the divergence volumetric integral:
Equations (11), (20), and (24) are certainly not unique, but they are fairly straightforward. For instance, Eq. (11) represents a facecentered secondorder approximation to a surface integral. Equation (20) represents a vertexbased volumetric integral consisting of a dotproduct contribution from each pair of vertex vectors. Equation (24) is a particularly simple secondorder approximation which gives all of the weight to the cellcenter value of while using a surfaceintegral formulation for that is analogous to the surfaceintegral used in Eq. (11).
Substituting from Eqs. (11), (20), and (24) into Eq. (10), we obtain the discrete version of Eq. (10):

(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 and values. In particular, the equation for the facearea component of on any given cell face is obtained from Eq. (25) simply by setting the same facearea component of on that face to unity and setting the remaining facearea components of 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 facearea components of , i.e., h^{R} , h^{B} , h^{T} , h^{D} , h^{U} , to zero:

(26) 

(27) 

(28) 

(29) 

(30) 

(31) 
Having derived Eq. (35), we can construct the discrete equation for the cellcenter 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 is:
where V denotes the total volume of the cell, the facearea flux components are expressed in terms of the intensities via Eq. (35), and Q^{C} denotes the source or driving function evaluated at cellcenter. We have chosen not to discretize the time derivative in Eq. (38) simply because essentially any standard discretization, e.g., the backwardEuler and CrankNicholson 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 7point stencil.Now that we have defined the equations for the cellcenter intensities, we must next define equations for the facecenter intensities. Our local indexing scheme admits two intensities and two facearea 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 + , j, k) is associated with the two intensities, and , and the two facearea 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 interiormesh facecenter 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 facearea flux component at each interior cell face. In particular, we require that
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 facearea 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 interiormesh 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 . 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 13point stencil. The only intensity shared by these two cells is . Thus in a certain sense it can be said that is ``chosen'' to obtain continuity of the facearea flux components on cellface (i + , 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 interiormesh facecenter intensity to the two cellcenter intensities adjacent to it. This enables the facecenter intensities to be explicitly eliminated, resulting in the standard 7point cellcentered diffusion discretization. This is completely analogous to the 2D case discussed in detail in [1]. However, if the mesh is nonorthogonal, the facecenter 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:
Having defined the equations for the cellcenter and interiormesh facecenter intensities, we need only define the equations for the facecenter 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 facecenter intensity and one facearea flux component associated with each such face. The equation for each boundary intensity is very similar to that for each interiormesh facecenter intensity in that it expresses a continuity of the facenormal flux component. The only difference in the boundary equations is that the analytic boundary condition for the diffusion equation is used to define a ``ghostcell'' facenormal flux component that must be equated to the standard facenormal flux component defined by Eq. (35). A ghost cell is a nonexistent 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 socalled ``extrapolated'' boundary condition. This condition is of the mixed or Robin type and can be expressed as follows:
where d^{e} is called the extrapolation distance, is called the extrapolated intensity (a specified function), and denotes an outwarddirected 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 , may change as a function of position. One obtains a vacuum boundary condition when = 0 , a source condition when is nonzero, and a reflective (Neumann) condition when = . The extrapolated boundary condition is said to be a Marshak condition whenever d^{e} = 2D .We begin the derivation of the ghostcell facearea flux component by substituting from Eq. (2) into Eq. (48):
where ^{g} is the flux vector associated with a ghost cell. Next we recognize that the outwarddirected unit normal vector for a ghostcell must be identical to an inwarddirected unit normal vector on the outer surface of the problem domain. Thus where ^{g} denotes a ghostcell outwarddirected unit normal vector. Substituting from Eq. (50) into Eq. (49), we obtain: Next we solve Eq. (51) for the outwarddirected flux component associated with a ghost cell: 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 and the facearea flux component on that face is f^{L}_{1, j, k} . Evaluating Eq. (52) at the center of face (, j, k) and multiplying the resulting expression by the magnitude of the outwarddirected areavector on that face associated with cell 1, j, k , we obtain the desired expression for the ghostcell facearea flux component: where the extrapolated intensity and the extrapolation distance are assumed to carry the ghostcell index.We next obtain the equation for by requiring that the Right and Left facearea flux components for cells (0, j, k) and (1, j, k) , respectively, sum to zero:
Note that Eq. (54) is identical to Eq. (45) with the latter equation evaluated at i = 0 . Thus Eqs. (45) through (47) provide all facecenter intensity equations with the caveat that when an intensity is on the outer mesh boundary, the associated ghostcell 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 7point stencil. This completes the specification of our diffusion discretization scheme.To summarize,
It is interesting to note that the equation for a cellcenter intensity contains a timederivative of that intensity, but the equations for the facecenter equations do not contain any form of time derivative. Thus in timedependent calculations, one must have initial values for the cellcenter intensities, but initial values are not required for the facecenter intensities. Thus only cellcenter 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 positivedefinite. We demonstrate this latter property in the Appendix.