The purpose of this paper is to present a local support-operators diffusion
discretization for arbitrary 3-D hexahedral meshes. We use the standard
finite-element definition for hexahedra [2]. The method that we present is a
generalization of a similar scheme for 2-D *r* - *z*
quadrilateral meshes that was
developed by Morel, Roberts, and Shashkov [1]. The diffusion equation that we
seek to solve can be expressed in the following general form:

We define a cell-centered diffusion discretization scheme as one that numerically preserves the integral of Eq. (1) over each spatial cell. In particular, substituting from Eq. (2) into Eq. (1) and integrating that equation over a cell volume, we obtain:

whereIf one considers only non-orthogonal meshes with material discontinuities, existing vertex-centered diffusion discretizations are generally more advanced than cell-centered discretizations. This is primarily so because of the enormous success of Galerkin finite-element methods [2] and variants of those methods. Nonetheless, there are applications for which cell-centered schemes appear to yield superior accuracy relative to vertex-centered schemes. For instance, when coupling radiation diffusion calculations with cell-centered hydrodynamics calculations, a cell-centered diffusion scheme is highly desirable because it avoids certain difficulties associated with vertex-centered diffusion schemes [4]. Our new scheme has been developed with coupled radiation-diffusion/hydrodynamics applications in mind.

The discretization scheme that we have developed is cell-centered, but it has intensity unknowns at both cell centers and face centers. It can be applied on both structured and unstructured meshes consisting of combinations of arbitrary hexahedra and arbitrary degenerate hexahedra (i.e., wedges, pyramids, and tetrahedra). It yields second-order accurate solutions for the intensities on both smooth and non-smooth meshes even when material discontinuities are present, and it generates a sparse symmetric positive-definite coefficient matrix.

The literature relating to cell-centered diffusion discretization schemes for
arbitrary hexahedra is not particularly extensive. One of the earliest relevant
papers appeared about ten years ago. In particular, Rose developed a cell-centered
hexahedral-mesh discretization scheme for the Laplacian operator.[5] The
diffusion operator that we consider degenerates to the Laplacian operator when the
diffusion coefficient is everywhere unity. Unlike our scheme, which has only the
normal component of the current on each cell face, Rose's scheme has three
components of the flux on each cell. Furthermore, the flux is continuous across each
cell face in Rose's scheme, whereas only the normal component of the flux is
continuous in our scheme. A central aspect of Rose's method is the preservation of
an integral expression that is referred to as an energy principle. Our method is
actually based upon the preservation of an integral identity. The energy principle
used by Rose is not the same as the integral identity that we use, but they are
related. In particular, the principle used by Rose can be derived from the diffusion
equation together with the integral identity that we use. Rose presented a proof
that his hexahedral-mesh method converges with second-order accuracy, but he
provided computational results only for a 1 - *D*
version of his method. Arbogast, et
al., [6] have recently developed a cell-centered expanded mixed
finite-element method for solving the tensor diffusion equation on general meshes
(including hexahedral meshes.) Their method has only cell-center intensity unknowns
if both the mesh and the diffusion tensor are smooth, but additional face-center
intensities are required wherever the mesh or the diffusion tensor is non-smooth.
The coefficient matrix generated by their method is always symmetric positive
definite (SPD). The method of Arbogast, et al., actually shares some of the best
properties of the standard mixed finite-element method and the hybrid mixed
finite-element method. Standard mixed finite-element diffusion methods have only
cell-center intensities, but this is achieved at the cost of solving a
computationally expensive saddle-point linear system. The saddle-point system can be
avoided by using the hybrid mixed finite-element approach, which generates a
symmetric positive-definite coefficient matrix at the expense of additional
face-center unknowns. The method of Arbogast, et al., yields an SPD coefficient
matrix like the mixed hybrid method but can sometimes require far fewer unknowns.
Although they proved several convergence theorems for their hexahedral-mesh method,
Arbogast, et al., provided computational results only for a 2-D version of their
method.

Our local support-operators method is similar to hybrid mixed finite-element methods in that it is cell-centered, it has both cell-center and cell-face intensities, and it produces a coefficient matrix that is symmetric positive-definite. However, our scheme is fundamentally a finite-difference technique since basis functions never appear in our formalism. The similarities between our method and mixed hybrid finite element methods suggest that there may be a connection between them, but we cannot presently demonstrate such a connection.

To summarize, the following combination of characteristics appear to be unique to our support-operators diffusion discretization scheme:

- It is a cell-centered discretization for arbitrary hexahedral meshes.
- It gives second-order convergence of the intensity on both smooth and non-smooth meshes both with and without material discontinuities.
- It generates a sparse SPD coefficient matrix.
- It is equivalent to the standard 7-point cell-center diffusion discretization scheme when the mesh is orthogonal.

The remainder of this paper is organized as follows. We first explain the central theme of our local support-operators method, and apply it to an arbitrary hexahedral mesh in Cartesian geometry. We next describe an approximate version of our scheme that we use as a preconditioner in conjunction with a conjugate-gradient solution technique [7]. Finally, computational results are given, followed by a summary and recommendations for future work.