Discretization of Anisotropic Convection-Diffusion Equations, Convective M-Matrices and Their Iterative Solution

We derive the constant-j box method discretization for the convection-diiusion equation, rj = f, with j = ?ru + u. In two dimensions, is a 2 2 symmetric, positive deenite tensor eld and is a two-dimensional vector eld. This derivation generalizes the well-known Scharfetter-Gummel discretization of the continuity equations in semiconductor device simulation. We deene the anisotropic Delaunay condition and show that under this condition and appropriate evaluations of and , the stiiness matrix, M, of the discretization is a convective M-matrix. We then examine classical iterative splittings of M and show that convection (even convection dominance) does not degrade the rate of convergence of such iterations relative to the purely diiusive (= 0) problem under certain conditions.


INTRODUCTION
We consider the two-dimensional convection- diffusion equation V.j=f, j=-aVu+flu (1.1)   with standard boundary conditions.Here a,/3, and f are spatially dependent on the defining domain f; the diffusion coefficient a is a two by two symmetric and positive-definite tensor field, and the convection coefficient /3 is a two dimension- al vector field.Equation (1.1) arises in semiconduc- tor device simulation as the (electron and hole) continuity equations and allows the mobility and diffusivity coefficient functions a greater generality.
In 2 we derive a box method discretization for Eq.(1.1) on a meshed domain fM which is a union of triangles.We use the constant-j assumption on triangle edge-pairs in order to resolve the vector *Corresponding author.Tel.: (919) 660-6510, e-mail: djr@cs.duke.educurrent density j= [jl, j2]T" The use of edge-pairs is essential when diffusivity of the media is aniso- tropic.The edge-pair approach allows us to ap- proximately integrate Eq. (1.1) on a "box" via the divergence theorem in the standard manner.The various parts of the approximate integration give a linear relation between a nodal value ui of the unknown function u and its nodal values at the various triangle neighbors and ultimately the com- pletely assembled linear system, in particular, the stiffness matrix M.
The use of the constant-j assumption dates from the seminal paper [29] by Scharfetter and Gummel (SG), and is well-appreciated in the semiconductor simulation community.Bank, Rose and Fichtner [5] show in III.C, Eq. ( 55), that a essential ingredient of the SG constant-j derivation can be interpreted as an edge evaluation scheme for e.This property also allows a finite element version of the SG discretization (see [38,16,2,3]).
In a box method (or finite volume) scheme, it becomes natural to apply the constant-j idea directly generalizing the discussion in [5], III.D.
As mentioned above, in order to resolve both components ofj, we are led to apply the constant-j assumption on edge-pairs.In the case where a aI2 x 2 is a scalar, one could directly apply the constant-j idea by assuming the projection of j on the edge is constant.If the vector field j is needed, however, this projection assumption is insufficient.
Interestingly both the edge-pair constant-j (SG) discretization and the finite element (SG) discre- tization give exactly the same stiffness matrix M for comparable evaluation schemes for a and/3.
We feel that our approach also has value when meshes are no longer triangular or tetrahedral in 3D, and we will suggest a generalization in 6.
Our box method seems peculiar, at first, since each triangle edge has two currents which may not be conservative (/j,t+ It,jO, in general, for edge (j, l)).This is due to the discretization pa- radigm which only requires current conservation on the whole box, and hence allows (small) de- viations locally.However, there are conditions on the triangular mesh and evaluation schemes for a and / which allow the then conservative cur- rents to be viewed edgewise.These conditions are also needed for matrix to be an M-matrix.
In 4 we examine, in some detail, the M-matrix nature of the system matrix M. We show, for example, that any irreducible column diagonally dominant M-matrix can be interpreted as a gen- eralized resistive network where the current on edge (j, l) is defined as/j= auj-but, with a and b both being positive.Clearly when a b, the general re- sistors become ordinary resistors.The quantity c (a-b)/2 is identified as the discrete version of convection.Furthermore, the discrete analogies of the curl and divergence of the (scaled) convection vector field can represented in terms of the discrete edge conductance and edge convection as well as the incidence and cycle matrices associated with the convection-directed graph.
Our analysis in 5 shows that under certain conditions on the scaled convection vector field b=aand the vector field fl itself, standard classical iterative methods converge faster when there is convection (fl-0) than when there is not (/3 0) in Eq. (1.1).Of course, as convection begins to dominate, it becomes necessary to refine the mesh to obtain comparable accuracies of the discrete solution to the corresponding solution of (1.1).What is perhaps surprising is that, for fix- ed accuracy, the necessary number of iterations is bounded (or even slightly decreases) as convection increases.
It is historically interesting that our iterative analysis in Section 5 could be regarded as an ex- tension of the Kahan-Varga theory of succes- sive over-relaxation (SOR); see Varga [36], 4.4, and Kahan [22, 23].In Kahan-Varga theory one gives up the consistent ordering assumption (Property A) and replaces it with an assumption which, in our context, follows from the discrete curl-free condition (3.13).While consistent order- ings are still best in this theory, when they exist, all orderings give boundedly related rates of convergence (for Jacobi, Gauss-Seidel and SOR).
We show how these convergence rates change with convection and relate the convergence rate of a convective problem directly to the same problem with convection identically zero.The discrete curl- free condition also implies that there are orderings which put the predominance of the discrete convection in the upper or lower triangle of the stiffness matrix, but a detailed look at our analysis shows that our bounds do not depend on these convection-related orderings; instead, any order- ing preference for a convective problem is simp- ly inherited from the zero-convection (diffusive) problem.Generalizations to more sophisticated iterative methods are clearly possible.
Finally, in 6 we present briefly another applied problem from computational physiology where the anisotropy of diffusivity (2 2 nature of c) is important, and make some additional concluding remarks.

Mesh Description
A discretization method for the model problem (1.1) in two spatial dimensions usually requires a mesh, in which the defining domain f is de- composed into a number of smaller simple cells.
A mesh can be either structured, such as grid graphs, or unstructured, such as triangular meshes.Our box method can be applied to both types of meshes.We will, however, present the box meth- od with the use of unstructured triangular meshes to demonstrate its greater generality.In the rest of the paper, we assume that a mesh is a triangula- tion over the defining domain, 9t, a compact and connected two-dimensional region on which the partial differential Eq. (1.1) is defined.Let P= {Pl,P2,...,PlPI) denote the set of mesh points, E--{el, e2,...,elEI) the set of mesh edges, and T {tl, t2,..., tiT I) the set of triangles.More spe- cifically, a point, an edge and a triangle are represented as 2. THE CONSTANT-j BOX METHOD In this section, we present the constant-j box method for discretization of Eq. (1.1) in two spa- tial dimensions.In two (and higher) dimensions, the main difficulties for discretization are how to handle the domain geometry and how to resolve the anisotropy of diffusion.We use unstruc- tured meshes to in order to handle arbitrary do- main geometry in two dimensions.Details of mesh generation and adaptation in two dimensions can be found in [32], Chapter 5.An important fea- ture of our box method in two dimensions is the use of edge-pair constant-j assumption which na- turally resolves the anisotropy of diffusion, hence that name "constant-j box method".
In 2.6, we show the equivalence between the box method and a specific finite element method on the anisotropic convection-diffusion equations in two dimensions.Based on this equivalence, we are able to obtain an error estimate for the box method using finite element analysis as in [38].
Pi (xi,yi) E 2, ej (uj, v), < u, v < np, for edge PujPvj, tk (ak, b, c), < a, b, ck < np, for triangle Apakpbkpck, respectively.We require the mesh to be a conforming triangulation in a sense that it preserves the boundary and interfaces of the domain f; that is, the domain boundary and the interfaces can be recovered by a subset of edges in E. The inter- faces in the defining domain are usually used to separate areas with different materials or block flows.
The mesh serves as the base of approximation.The unknown state function u is approximated at mesh points by a vector, uh, by means of u(pi) be- ing approximated by ui, the ith component of Uh.Since the values of u at the Dirichlet boundary points are already known, the approximate vector uh needs to be defined only on the reduced set of the mesh points, Pr, a point set that excludes the Dirichlet points from P. In special problems where flux blocking interfaces are included, a mesh point on a interface (end points excluded) is duplicated in Pr, with one on either side of the interface, and is treated as a zero-flux Neumann boundary point; (see [32], Appendix I).unknown variable is needed, only the boxes that correspond to points in the reduced set Pr will be needed for discretization.These boxes will be referred to as the effective boxes.

Box Formation
To apply the box method, we need to construct a secondary mesh, a set of boxes, on top of the triangular mesh.We partition the defining domain into a set of non-overlapping boxes, {B1,..., each around a mesh point.We tentatively form the boxes by selecting a point, called a box vertex in each triangle.We then connect the box vertex with the three mid-points of the triangle edges (thinner lines in Fig. 2.1) to form three box edges (thicker lines in Fig. 2.1).Any two of the three box edges form a boundary patch (shaded lines in Fig. 2.1) which will become a part of the box boundary that surrounds a vertex (p; in Fig. 2.1) of the triangle.Typically the circumcenter (intersection of the perpendicular bisectors) or the centroid of the tri- angle is chosen.And the boxes formed in this way are called Voronoi boxes and barycentric boxes, respectively.An interior box, one that covers an interior mesh point, is bounded by box edges only, while a boundary box, one that covers a boundary mesh point, is bounded by both box edges and (two) boundary edges of the defining domain.Hence, the boundary boxes are also called "half boxes".Since only one equation for each After the boxes formed, we integrate the model Eq.(1.1) on each effective box Be, V .jdaf f da (2.1) for 1,..., n.The divergence theorem in vector calculus (see [24]) implies that v .jaa where OBe is the boundary of box Bi, and u is the unit outward normal vector of the box boundary..3) Bi The model problem (1.1) is now equivalently converted into a set of integral equations (2.3), or the local integral forms, on all effective boxes.
The area integral of the source term on the right-hand-side of Eq. (2.3) is approximated by applying the following numerical quadrature rule [10] that yields the desired degree of accuracy: f da wikfk, (2.4) where Ai is the index set of the neighbors (mesh points that are directly connected to Pi) of Pi, Wik is the quadrature weight which determined by the quadrature rule, and fk is the nodal value of the source term f at pk.Notice that for higher degree quadrature schemes, values of the integrand at points other than the mesh points (midpoints of triangle edges for example) may be needed.The consistency condition requires that the quadrature approximation (2.4) be exact when the integrand f is a constant, or Z wik area(B/).
(2.5) kEAiOi Often the value of wi is taken as for k, and 0 otherwise.This procedure is often referred to as the lumping, and is equivalent to assuming the integrand of the area integral is box- wise constant.Calculating the box areas in terms of the triangular mesh and the boxes is straight- forward..4.Edge-pair Constant-j Assumption Approximation of the line integral is more com- plex.It is carried out at two levels.First at the edge- pair level, we assume the scaled current density n= a-lj constant at a pair of edges, or an edge- pair, ((i, k), (i, l)) for each triangle (i, k, l) incident to the mesh point Pi (see Fig. 2.2).Projecting onto the two edges of the edge-pair, we have + bu) hiktli k, hit (--U -}-bu) hiTlik, where b a-lfl is called the scaled convection.If we associate with each edge a constant vector approximating the scaled convection b, e.g., bik for edge (i,k) and bil for edge (i,/), the above set of two ordinary differential equations along the two edges can be solved exactly as hi1li k n(--/ik)U n(/ik)Uk, hi;lik n(-/il)Ui n(/il )Ul, ( T -1 Aik hikbik hikOik /ik, hilOii il, are the scaled edge convection scalars, or edge convection in short; aik, ik, ail, and /il the evaluations of a and /3 at edges (i,k) and (i,l) respectively.Equation (2.8) is a direct result from solving the two one-dimensional first order ordi- nary differential equations in (2.7)  (2.10) being the edge-pair matrix.
Remark 2.1 We note that the constant-j assump- tion is just one of many techniques that can be used to approximate the projection of (scaled) current density along an edge.For example, one may replace the Bernoulli function used in Eq.
(2.8) by an appropriate function to obtain an upwinding type of discretization along an edge.
As long as the function used to compute the coefficients of u in Eq. (2.8) is always nonnegative, the monotonicity of the final discretization will not depend on the convection of the problem as we shall see later.
Remark 2.2 In the case where the scaled convec- tion is curl-flee (see 3.4), we are able to solve the above set of ordinary differential equations (2.7) exactly without assuming the scaled convection b is constant on either of the edges.The solution is still in the form of Eq. (2.8) with "ik .3(?i) (2.11) where b -7b for some potential function b.This result can also be extended to the box method in three dimensions.
Remark 2.3 When the model (1.1) has no con- vection, i.e., /30, the scaled current density is simply the negative gradient of the unknown state function, -Vu.Furthermore, the three constant vectors, namely lik, ikl, and tkti, approximating the negative gradient on the three edge-pairs of a triangle are identical, which means that the box method approximates the gradient by a single constant vector on each triangle.This obser- vation is also valid in three dimensions.
Equation (2.9) enables us to approximate the edge-pair flux as uTj ds ,, fmo u T ani ds. (2.13) The quasi-discrete form of Eq. (2.13) can be written in two equivalent ways: (2. is the pair of right and left quasi-discrete edge conductance scalars that are associated with the right edge (i, k) and the left edge (i, l) of the edge- pair ((i,k), (i, l)) respectively, and f--6 lITads (2.17) is the quasi-discrete edge-pair conductance matrix.
Approximating the quasi-discrete edge conduc- tance pair and the edge-pair conductance matrix are two equivalent approaches that allow the edge-pair flux in Eq. (2.13) to be completely dis- cretized.However, each approach provides a different interpretation.Here, we only present the derivation of the discrete edge conductance pair approximating (2.16).For derivation of the discrete edge-pair conductance matrix approxi- mating (2.17), see [32], 3.2.4 and also [33].
The quasi-discrete edge-pair conductance scalars (2.16) can be discretized by assuming the diffusion tensor a is constant in some manner.
We may choose to use different constants approxi- mating a for the two edge-pair conductance scalars in Eq. (2.14), or we may choose them to be equal.These two choices lead to the following two different evaluation schemes for the diffusion coefficient.
1. Per edge conductivity evaluation we evaluate the conductivity (diffusivity) tensor a by two edge-associated constant 2 2 matrices, aik at edge (i, k) and ail at edge (i, l), to approximate the right and the left quasi-discrete edge con- ductance scalars in Eq. (2.16) as area integral on the right-hand-side of Eq. (2.3) still depends on the coordinates of box vertex o.Let us define the generalized dot product, cross product and angle between two vectors u and v with respect to a symmetric and positive-definite where s and u are respectively the length and the unit normal vector for a box edge; see Figure 2.2.edge-pair conductivity evaluation we evalu- ate the conductivity tensor a by a single con- stant 2 2 matrix, auk, associated with the edge-pair ((i,k), (i, l)) to approximate both the right and the left quasi-discrete edge conduc- tance scalars as A convenient special case here is to let the three edge-pairs of a single triangle all give the same evaluation result, which then is called the per triangle evaluation scheme.
(2. 25)   which can be thought as the normal dot product, cross product and angle measured in the inner product space with the inner product (2.23).Then the edge conductance values can be interpreted as g,R/ v/det(aev)cot(Ozv) R is the generalized angle between edge where O;v I is the generalized angle vector hil and hl, and between edge vector hi and ht.
With the scaled edge convection and the edge conductance well-defined, we are able to approxi- mate the edge-pair flux in Eq. (2.13) by the edge- pair current -Jr-glLik(B(--, il)Ui B() il)Ul), (2.28)   where Oev can be substituted with any chosen evaluation of the diffusion tensor.The above Eqs.(2.21), (2.22) show that the discrete edge con- ductance scalars do not depend on the choice of box vertex o.They depend on only the local mesh geometry, the triangle edge vectors, and the evaluation of the diffusion tensor.However, the which can be split into the right edge current glik(B(--/ik)Ui B(/ik)Uk) (2.29)   and the left edge current glil(B(-Ail)Ui-B(Aitc)uI).
(2.30) such that the net current at Pi is expressed as the sum of edge currents Iik for edges (i, k) incident to pi.Based on their locations, the mesh edges can be divided into boundary and internal edges.An internal edge (i, k) is shared by two edge-pairs that are incident to point Pi (Fig. 2.3 left), and its edge current is the sum of the right edge current for the edge-pair to its left and the left edge current for the edge-pair to its right: A boundary edge (i,k) is attached to only one edge-pair incident to point pi (Fig. 2.3 right), and its edge current is either the right or the left edge current of the edge-pair it is attached to, depending on edge orientation with respect to the defin- ing domain f: Iik IkLo. (2.33) Furthermore, since the scaled edge convection is computed by evaluating the diffusion and con- vection coefficients at edges, it is invariant on an edge regardless which edge-pair that contains the edge is being considered.As a result, we are able to write the edge current from pi to pk as By considering all the nodal net currents al- together as the box method approximation of the current density divergence over the entire defining domain f, we have a linear system Muh fh + b (2.37)   with the stiffness matrix M approximating the convection-diffusion operator V. (-aV +/) in the model Eq.(1.1), the vector fh approximating the per box area integral of the source term f, and the vector b approximating Dirichlet and Neumann boundary conditions assigned to the boundary of the defining domain.
We can construct the stiffness matrix M through edge assembly.Based on their incidence relation hik "ik Pi FIGURE 2.3 An internal edge (i, k) shared by two triangles (left) and a boundary edge (i, k) attached to only one triangle.
with Dirichlet boundary points, we divide the edge set into non-grounded edges E and ground- ed edges Eo.
For a floating edge (j,/), both / and I O. are used for integration around point py and point p.Therefore, its contribution to the stiffness matrix M is the 2 2 edge stamp gjB(-j) -gjB(At) ] (2.38) sit= _goB(Aj goB(_AIj) jr with use of the notation with ei, and ej,n being the ith and jth columns of the n n identity matrix respectively.If Pt is the Dirichlet point for grounded edge (j,l), the grounded edge (j,l) contributes to the stiffness matrix M by a edge stamp sit gjt [B(-Ajt ]j  (2.42) Then it is clear that the stiffness matrix can be written in the edge assembly form as M--y sj,  The edge-assembly view of the stiffness matrix M provides a circuit interpretation of the box method discretization, and we pursue this inter- pretation in more detail in 3.3.There we will orient each grid edge (j, l) such that (j, l) will be- come an edge in the directed graph G(M) with Aj > 0 using Aj =-A O.

Equivalence to a Finite Element Method
For non-convective problems, the equivalence be- tween the box method and the piecewise linear finite element method in one and two dimensions has been well described in [4, 18, 20].The use of piecewise linear base functions in finite element methods can be considered exactly equivalent to the constant-j assumption of the box method (see Remark 2.3).More recently, Xu and Zikatanov [38] presented the edge-averaged finite element (EAFE) scheme which gives the same left-hand- side discretization as does the box method for two-dimensional isotropic convection-diffusion problems when the evaluation schemes for the coefficients a and/3 used in the two methods are consistent.See also Gatti, Micheletti and Sacco [16] for a similar discussion.
In this section, we generalize the EAFE scheme to the anisotropic EAFE scheme such that it will be equivalent to the box method for general anisotropic convection-diffusion problems again assuming consistent evaluation schemes.
Suppose we are solving the following boundary value problem V (-aVu + flu) f (2.45)   on the defining domain f with zero Dirichlet boundary condition.The weak formulation of the above problem is to find a function u in the Soblev space H(f), such that a(u, v) =_ f(au flu)q-(7v)da =f(v) (2.46) for every trial function v E H (f). Let Vr denote the space of piecewise linear functions on trian- gles in the triangular mesh T, with a set of basis functions defined as vi(Pj)---6ij i,j=l,...,n. (2.47) Then the finite element formulation is to find a piecewise linear function Uh E Vr that approxi- mates the true solution to Eq. (2.45) such that ah(u, vh) Z t" (oVu -/3u) q-(Vv)da =fn tET (2.48)   with fh--f(Vh), for all trial functions Vh6 Vr.
The lehand-side integral in (2.48) on a triangle (i, k, l) can be written as where det(cQ (hice_lhik) are the edge conductance scalars for edges hkt, hli and hig respectively (cf.Fig. 2.2), Eqs. (2.21), (2.22)).Here we have used the identity (2.19) in [2].Using the edge average approximation Eq. (3.1!) in [38], this integral (2.49) can be approxi- mately written as Z ge't"'e6e(e-eUh)te(Vh)' (2.53) ect where e takes edges hkl, hti and hik, 6e is the edge differencing operator, and the edge-associated quantities g, and ,,/are defined as gkl,t -1 ,cl,tda :r(f ) (2.54) 4It] 2 hti det(a)cda hi, ' h;-lds, (2.55) for edge ht, and cyclically for edges hi and hi.In reality, these integrals may have to be computed approximately.One simple approach is to assume = t and fl= fll are constant along edge (k, l) in these integrals, which leads to the following approximation of (2.54), (2.55), (2.56): This assumption corresponds to the per edge evaluation scheme of the constant-j box method.
.61) tDe=(i,k) It is easy to verify that the expression of mjt above is exactly the same as the (2, 1) entry of the edge stamp (3.3).For Pi not connected to a Dirichlet boundary node, we have mii Z Imj, I, (2.62) (i,k)eE which also coincides with the result from the box method discretization.Equivalence between the box method and the above anisotropic EAFE scheme at boundary nodes also can be easily verified.If we choose to use the finite element approach to obtain the approximation of the right-hand-side (area integral) in Eq. (2.3), then the box method and the anisotropic EAFE meth- od are completely equivalent.
Using the same finite element analysis provided in [38], we obtain the following a priori error bound.
THEOREM 2.4 Let u be the solution of the problem (2.45).Assume that for all E T, a is symmetric and positive-definite whose eigenvalues are uniformly bounded away from zero, all components of a and are in W ' (t), and the scaled current density =_ -Vu + bu [wl'2(t)]2.Then the following esti- mate holds

MATRIX PROPERTIES
In this section, we will identify a set of conditions, which depends only on the edge conductance and convection evaluation schemes as well as mesh properties, that give rise to the M-matrix proper- ties of the discretization system of (1.1).Further- more, we define discrete analogies of the curl and divergence of the (scaled) convection vector field in terms of the discrete edge conductance and convection.We show that certain conditions on the discrete curl and the discrete divergence imply that the stiffness matrix M is symmetrizable or row diagonally dominant.Although these convection- related properties do not affect the monotonicity of the box method discretization, they do arise in other discretization methods to ensure monotoni- city, and they appear as important conditions in our convective iteration analyses in 5. ) Ilu uhlll, < Ch (2.63) tT for sufficiently small h.
Remark 2.5 We notice that the above result is based on the exact computation of the quasi- discrete quantities-edge conductance ge,t, edge convection potential e and its exponential edge average 7e-which are in the form of integrals of the diffusion and convection coefficients.Therefore, the error analysis needs to be interpreted carefully when the coefficients become irregular.
The above error estimator is called a priori because it is generally true for a class of problems in the form of Eq. (2.45) with only smoothness restrictions on the diffusion and convection co- efficients.[31] presented a tight one-dimensional a priori error analysis on the box method for constant-coefficient convection-diffusion equations.It shows that the maximum norm of the solution error is bounded by O(hZlbl) where b /3/a is the scaled convection.

M-matrix
A real matrix A is said to be nonnegative, denoted by A _> 0, if all its entries ajl >_ O. Let Z nn denote the class of Z-matrices, square matrices with non- positive off-diagonal elements, i.e., Z nxn {A [aft] gnxnla fl <_ O, j l}.
A general definition of M-matrices is stated as follows.
DEFINITION 3.1 ([7] Chapter 6) A real square ma- trix A is called an M-matrix if and only if A can be written in the form, A=sI-B, s>O, B>0, (3.1)   where s _> p(B); p(B) is the spectral radius of B.
The spectral radius of a square matrix is defined as the largest modulus of all its eigenvalues.A nonsingular M-matrix requires s > p(B), hence A -1 exists.Equivalently, a nonsingular M-matrix is a Z-matrix with nonnegative inverse.A sym- metric nonsingular Z-matrix is an M-matrix if and only if it is positive-definite; such a matrix is called a Stieltjes matrix.For a Z-matrix A, there are other equivalent conditions for A to be a nonsingular M-matrix [7], such as There exists a positive diagonal matrix D1 such that DIA is strictly column diagonally domi- nant; There exists a positive diagonal matrix D2 such that ADz is strictly row diagonally dominant; There exists positive diagonal matrices D1 and D2 such that D1AD2 is both strictly row and column diagonally dominant.
When A is irreducible (equivalently, the graph G(A) defined in 4.1 is connected), then the above conditions can be changed to irreducibly row (column) diagonally dominant; see Varga [36], page 23, and Axelsson [1], Definition 4.4.An irreducible M-matrix has A -1 > 0, i.e., each entry of A-1 is positive.

Nonnegative Column Sum and Edge Current Conservation
From the electrical circuit viewpoint, it is natural to assume that the current is conserved along the edge (j, l).Using Eq. (2.34), such an assumption corresponds to Ifl + Io.O e: gjt g O, ( which we call the edge conservation condition. LEMMA 3.2 If we use the per edge or the per triangle evaluation scheme to evaluate the diffusion tensor field c when computing the edge conductance scalars, then for each edge (i,k)EE, we have gik=gki, i.e., the edge conservation condition is satisfied.
Notice that we always evaluate c and/3 edge- wise when computing the edge convection Ajt.Then the 2 2 edge stamp that is used to assem- ble the stiffness matrix M in Eq. (2.43) can be simplified as B(-/jl) gfl -B(--jl) (3.3)   for edge (j,/), noticing that Ayt=-'0.Since the column sum for each edge stamp is zero, it is obvious that the column sum of the stiffness ma- trix is either zero when the column corresponds to a non-grounded node or positive otherwise.Finally, we show the following result.
LEMMA 3.3 An irreducible square matrix with nonnegative column sums and at least one positive column sum is an M-matrix if and only if it is a Z-matrix.
Proof If a square matrix A is a Z-matrix and has nonnegative column sums, then A is column diagonally dominant.Furthermore, if A has at least one positive column sum, A is irreducibly column diagonally dominant.Applying the Mmatrix conditions given in 3.1, we have A is an M-matrix.The other part of the proof is obvious since an M-matrix is always a Z-matrix.

M-matrix and Nonnegative Edge Resistors
From classical partial differential equations theory, we know that the elliptic operator V.(-aV+/3) satisfies the maximum principle, which states that the maximum and the minimum of the solution u to the boundary value problem, X7. (-aVu +/u) O, occur only at the boundary of the defining do- main (f. [35,21]).The discrete analogy of such property will be that the stiffness matrix M is inverse-monotone or simply monotone (cf.[26]), i.e., Mv>O= v>O, for any vector v.A sufficient condition for M to be monotone is that M is an M-matrix, (cf.Defini- tion 3.1).The M-matrix property of the stiffness matrix M is equivalent to the nonnegativity of edge conductance scalars, i.e., gic >_ 0 (3.4)   which is consistent with the diffusive nature of the elliptic operator 7. (-a7 + fl) being approxi- mated.Since, in addition, a monotone scheme is numerically stable (no spurious oscillations), monotonicity of the stiffness matrix M is desired for discretization schemes on elliptic partial dif- ferential equations.We have seen from simula- tion examples that discretized systems without the M-matrix property may cause artificial spurious wave fronts in reaction-diffusion systems such as the heart models. I addition to capturing the monotonicity of the physical problem and ensur- ing numerical stability, many iterative methods for solving M-matrix systems are guaranteed to converge.Furthermore, detailed analyses on such convergence is often possible as we shall see in 5.

Main Results
We now examine conditions on the mesh and ways to evaluate conductivity tensor a that will en- sure that the stiffness matrix M is an M-matrix.Both the algebra of edge-assembly (2.43) and the edge conservation condition (3.2) imply that the M-matrix property of the stiffness matrix M is completely determined by the nonnegativity of the discrete conductance scalars.
LEMMA 3.4 Under the edge conservation condition (3.2), the stiffness matrix M derived by the box method is an irreducible M-matrix (cf.Definition 3.1) /f and only if the edge conductance scalar gik is nonnegative for each (directed) edge (i,k) E. Proof Using Lemmas 3.2 and 3.3, we only need to show that the stiffness matrix M is a Z-matrix.
By looking at the formulae in Eqs.(2.38)-(2.43), it is clear that the edge conductance gi is nonnegative if and only if the corresponding off- diagonal entry in the stiffness matrix M is non- positive, due to the nonnegativity of the Bernoulli function B(x) (see Appendix A).Then M is a Z-matrix by definition, m Remark 3.5 The stiffness matrix is irreducibly dia- gonally dominant if there is at least one Dirichlet node.In this case the M is a Stieltjes matrix (cf. 3.1) if the convection vanishes identically in Eq. (1.1).
Using the formulae for edge conductance scalars we derived in 2, we have the following results.
THEOREM 3.6 If we use per edge evaluation scheme in the box method discretization, the stiffness matrix M is an M-matrix if and only if we have > 0 x ho.) 2 for each boundary or interface edge (i,k) (cf.Fig.
With respect to the edge conductivity evaluation cik, the inequality (3.5) is called the anisotropic boundary Delaunay condition for triangle Apipjpk that is attached to a boundary edge (i, k), and the inequality (3.6) is called the anisotropic Delaunay condition for the pair of triangles ApiPjPk and Apkppi that share an edge (i,k).When the dif- fusion tensor c becomes a scalar, we have the following known corollary [38, 3].
COROLLARY 3.7 When the diffusion tensor c in Eq. (1.1) degenerates to a scalar, then the stiffness matrix of the per edge evaluation scheme is an M-matrix if and only if any two triangles (i,j,k) and (k,l,i) that share an edge (i,k) satisfy the generic Delaunay condition or simply the Delaunay condition < PiPjPk q-< PkPlPi <_ 7r. (3.7) and any triangle (i,j,k) that is attached to a boundary or &terface edge (i, k) satisfies the generic boundary Delaunay condition or simply the bound- ary Delaunay condition < pipjpk < -. (3.8) As of this point, we have shown that a sufficient condition for the stiffness matrix M to be an M- matrix is that (i), we evaluate the conductivity c at edges and that each pair of triangles sharing a common non-interface edge in the underlying mesh satisfies the anisotropic Delaunay condition (3.6); and that (ii), each triangle attached to a boundary or an interface or a flux barrier satisfies the anisotropic boundary Delaunay condition (3.5).Such a triangular mesh is called an edgewise anisotropic Delaunay mesh with respect to the edge conductivity evaluation.
A special case of the edge conductivity evalua- tion scheme is to assign the same conductivity value for edge (i,k) shared by triangles Apipjpk and Apkptpi (ef.Fig. 2.3 left), or edge (j, l) shared by triangles Apipjpt and Apkplp, for edge (i,k) or (/',l) not being an interface edge.Such eval- uation scheme is also called the per quadrilateral conductivity evaluation scheme.A triangular mesh that satisfies the anisotropic Delaunay condition and the anisotropic boundary Delaunay condition locally and respects the quadrilateral conducti- vity evaluation scheme is called a pairwise aniso- tropic Delaunay mesh.
We see that the anisotropic Delaunay con- dition is only a condition on the "generalized angles" of the triangles, and not on their sizes.The significance lies in the fact that one is al- lowed to use an arbitrarily coarse mesh, if the mesh is pairwise anisotropic Delaunay, to get a qualitatively correct solution, which may be im- portant in iterative methods where a hierarchy of meshes is used.

Anisotropic Delaunay Mesh
We have shown that our discretization method is guaranteed a priori to be "monotone" for (1.1) if the mesh satisfies the anisotropic Delaunay con- dition.While we are aware of that implemen- tations of anisotropic meshes exist [11, 34, 8], we have also noticed that these implementations are mainly used for a posteriori error control of the numerical solutions; and these implementations do not always guarantee the exact anisotropic Delaunay condition (3.5) and (3.6) that is required by the constant-j box method for monotone discretization.Therefore, we chose to implement our own anisotropic Delaunay mesh generator to meet the condition for monotone discretization.Our mesh generator uses a generalized version of the Delaunay refinement algorithm which is described in [32], Chapter 5.For a given continu- ous diffusion tensor field on a closed and bounded domain (subdomain) with arbitrary geometry, the output mesh is guaranteed to satisfy the anisotrop- ic Delaunay condition.The density of the mesh mostly depends on how severely the anisotropy of the diffusivity changes in the domain, and it does not depend on the strength of the anisotropy.We believe that such restriction on mesh density is likely to occur in any implementation of the anisotropic Delaunay mesh generation.

A Simple Monotonicity-pvesevving Fix
In cases where a given mesh does not have all its edges satisfying the anisotropic Delaunay condition, the stiffness matrix may lose mono- tonicity, a fundamental property of the elliptic operator _7.c7 which the discretization should always inherit.Such matrices, when applied to the heart model equations, may even cause arti- ficial wave fronts.Here we propose the following simple post-processing procedure that will guar- antee that the output stiffness matrix is an Mmatrix, and therefore the discretization system is monotone: i: max{gi:, 0}. (3.9) By Lemma 3.4, it is obvious that the matrix assembled by discrete edge conductance gk is an M-matrix.(3.10) if the Einstein relation on mobilities and diffusi- vities holds; see Selberherr [30], and Fichtner, Rose and Bank [14].
We generalize the property of the electric field and the Einstein relation to the curl-free condition on the scaled convection vector field b a-l/3: V b=0 on the defining region f.Based on theo- rems from advanced calculus [24], we have the following equivalent conditions for b being curl- free on f.
For any connected open subset 9t0 of f, the line integral of -b, the tangential component of b, along the closed boundary loop of f0 is zero: There exists a differentiable real valued function, b, defined on f, such that the vector field b is the negative gradient of the function : b -Vb.(3.11)We consider the triangular mesh over the defin- ing domain f as a planar graph G (P, E).For each directed edge (j, l) in E, we assign a scalar p-pII Ajl --bds b(pj) b(pi), (3.12)   which is the integral of tangential component of b along the edge (j,l).We say that (G,A), a combination of the mesh and the scalars assigned to all the mesh edges, is a discretization of the vector field b.Then the condition is stated as follows.
discretized curl-free CONDITION 3.8 (Discrete Curl-flee) The discre- tization (G, A) of a vector field b satisfies the dis- cretized curl-free condition if and only iffor each cycle (il, i2,..., it, il) of the graph G, we have .13) Remark 3.9 If the graph G is a triangulation, the discrete curl-free condition (Condition 3.8) can be equivalently stated as Aik n t-Akl n t-Ali 0 (3.14)   for each triangle (i,k, l) in the triangulation.This is obvious because A triangle (i,k, 1) is a cycle in the graph G. Any cycle of the graph G encloses a set of tri- angles, and the sum of edge convection along the cycle is the same as the sum of triangle edge convection sums, which are the sums of edge convection around all triangles enclosed in the cycle.
We chose to state the discrete curl-free condition in the form of Condition 3.8 so that we would be able to apply the related results even when a non-triangular mesh is used.
There are two obvious ways to compute the edge convection A's on a mesh so that the dis- cretized curl-free can be satisfied when the scaled convection b in the physical problem is curl-free.
If the convection potential function b is giv- en, either analytically or discretely (with use of a background mesh), then the edge convec- tion At for edge (j, 1) can be computed as Ajt b(pj) b(pt), ( .11) and (2.12)).In the former case, the edge convection is computed as the exact line integral of the scaled convection vector field along the edge (j,/): fo IlhJ'l' --bds fo 'q-(-p)ds 3 pj 3 p

/j
If the scaled convection field b is given instead of the convection potential b, and if the line integral of the scaled convection vector field can be computed only by approximation, then computing the edge convection scalars b us- ing numerical integration, say the midpoint rule, such as in (3.12) will not guarantee the discrete curl-flee condition.Instead, we use the following approach: 1. Construct a spanning tree that spans the entire mesh.For each edge (j, l) in the span- ning tree with node j the parent of node l, assign to the edge the numerical integration of the scaled convection, e.g., Ajl h)Xlbl((pj+pl)/2). ( Also let bj bl + Ajl.We can simply set the b value at the root of the spanning tree to zero.
2. Assign an edge convection value A to each of the remaining edges in E such that the cycle formed by adding the edge to the spanning tree satisfies the discretized curl-free condition.
Based on both the curl-free condition on b and the discrete curl-free condition on A, it is obvious that the edge convection computed on edges that are not in the spanning tree should have same order of accuracy as those that are in the spanning tree.If the edge convection scalars in the box method discretization satisfy the discretized curl-free condition (3.13), we have the following result for the stiffness matrix M. THEOREM 3.10 If both the edge conservation condi- tion (3.2) and the discrete curl-free condition (3.13) are satisfied in a box method discretization of the model Eq.(1.1), the stiffness matrix M can be symmetrized by a positive column scaling.
Proof Using Eqs. (2.43) and (3.3), the stiffness matrix M can be written as g(l 3j) -g(3j gjl --B(l-Cj) B(j M (j,l)EUEo when both the edge conservation condition (3.Using Eq. (3.18), we also see that the stiffness matrix M is diagonally similar to a symmetric mat- 19) is a symmetric and positive-definite matrix, where the diagonal matrix D is defined in Eq. (3.17).
We conclude our discussion of the discrete curl- free condition with the following remarks.Remark 3.11 If the convection potential b is explicitly given in the model Eq.(1.1) as in the semiconductor device models, the symmetrized stiffness matrix S can be interpreted as the box method discretization of the symmetrized problem V. (-e-aVv) =f, v eu, (3.20)   using the third evaluation scheme for the exponential function ediscussed in [5].
Remark 3.12 Notice that Theorem 3.10 is still true even if the sign of the edge conductance gjt values or the ground conductance values in G g becomes negative.In fact, the discrete curl-free condition (3.13) is exactly equivalent to the com- bination of the "conservation law" and the "reversibility law" defined by.Parter and Young [25].The construction of the diagonal scaling matrix D /2) in [25] is also equivalent to our construction of the values of in the spanning tree manner we have described.

M-MATRICES AND DISCRETE CONVECTION-DIFFUSION
In this section, we present a circuit interpreta- tion of the stiffness matrix of the constant-j box method.We examine, in some detail, the M-mat- rix nature of the stiffness matrix M, using graph theory to discuss the combinatorial interpreta- tion of the curl-free condition.We also show how to interpret any convective M-matrix as a generalized resistive network.
an undirected graph G (V, E), where V {l;1,..., Vn} ( is a set of n vertices, corresponding to the rows (columns) of M, and E is a set of undirected edges: .2) Notice that E does not contain any grounded edges, i.e., edges that are connected to nodes with Dirichlet boundary conditions.We always assume that the graph G is connected, i.e., the matrix M is irreducible.Now we assign directions to the edges.The pair (j, 1) for edge ek is oriented in such a way that Im ll ImoI,(4.3) and we say that ek is an .out-edgefor vertex vj and an in-edge for vertex vt.See Example 4.6 below.We call this directed graph the convection- directed graph and often still use the notation G. Furthermore, we define the vertex-edge inci- dence matrix, or simply the incidence matrix (see Deo  [12]), Anm as A(j,k) -1 0 if ek is an out-edge for vj, if ek is an in-edge for vj, otherwise.
We can order the edges, and partition the edge set E into ET, (n-1) edges in a spanning tree of the undirected graph G, and Eo (m-n + 1) edges that are left, and write the incidence matrix A as be an n n convective M-matrix.We can define accordingly.As we know from graph theory, adding an edge in Ec to the spanning tree creates a fundamental cycle in the undirected graph G.
Therefore, we can define the cycle matrix (see Rose [27]), B(n-m+ 1)xm, as and Bo(i,k) -1 if edge ek edge in ET is in ith cycle, and has the same direction as ith edge in Ec; if edge ek edge in ET is in ith cycle, and has the opposite direction as ith edge in Ec 0 otherwise.
Each row of the cycle matrix corresponds to a fundamental cycle in G with nonzero entries indi- cating the edges forming the cycle.We define a generalized resistor connecting node j and node as a ordered pair of positive real numbers (a,b) such that the current from j to ! is described as while the current from to j is with ak >_ be.We disassemble the matrix M by edges, i.e., a 2 2 stamp (defined in Eq. (2.39)) -a b d -a b en (4.9) for edge ek=(j, 1).The edge stamp is simply the nodal admittance matrix for two nodes j and connected by a generalized resistor [a,b].Now the matrix M can be written as I 0 -I d. (4.7) Then we define a generalized resistive network as a connected graph where each edge represents a generalized resistor.In addition, there is at least one generalized resistor for some node v E V which is connected to a boundary set B. The graph Go (VU B, EU E0) contains G with E0 the set of grounded edges.Let A0 be the incidence matrix for Go, which adds additional columns to A. A column vector ej,n (-G-) is added if node j is connected by an out (in) grounded edge.
THEOREM 4.1 Any convective M-matrix is the nodal admittance matrix [12] for a generalized resistive network where the current Ik from node j to node is Ifl akuj bku! dk(Uj Ul) + Ck(Uj + Ul) being the diagonal dominance for jth column.
Since the convective M-matrix M is always column diagonally dominant, the ground conductance values are nonnegative, m Remark 4.2 Even though we always have the entry with bigger absolute value (ak) in the lower triangle for each edge stamp, the left-hand-side of (4.9), it may appear in the upper triangle of the matrix M when j > l.
Remark 4.3 In general, the diagonal matrix G g can be assembled as Proof Consider a convective M-matrix M and its graph G(V,E).For edge ek=(j,l) in E, let ak--Imyll and bk Imzyl be the absolute values of the corresponding off-diagonal elements of M, Gg Z [fk]j,, fk ag or f bg (4.11) ek Eo where ee is a grounded edge connecting node j of the graph G with node in the boundary set B.
Hence, we can write M more compactly than in (4.10) as M Z [ --a k a k -bkbk (4.12) k jl when this causes no confusion.If M alone is spe- cified, then only the regular resistor g is known for ek E E0 and can be interpreted in the general- ized resistor form (g,g).However, if both M and the generalized resistors on the grounded edges are specified, then G g is further decompos- ed as in Eq. (4.11).This is always the cases when M arises from discretization.respectively, the edge stamp can be split into a diffusive piece and a convective piece: Remark 4.5 If the graph G (V, E) of a convective M-matrix M can be embedded into a triangula- tion, called the induced triangular mesh, then we can construct the diffusion and convection coeffi- cients in Eq. (1.1) so that M is exactly the stiff- ness matrix of the box method discretization of the constructed partial differential equation using the induced triangular mesh.This is because that for each edge in G(V, E), we can define the fol- lowing two quantities (gk,/k) as ak bk )k ln ( ak ) (4.15) gk In ak In bk -k when ak > bk, and gk ak bk, Ak 0, (4.16)   in the box method discretization, then the edge stamps (3.3) used to assemble the stiffness matrix are exactly the same as the edge stamps (4.9) used to assemble M, and M can be written as B(-jt) -B(fl) 1 (4.17) gk -B(-/fl) B()jl) jr" M: ek EEt3Eo Notice that k representing a generalized resistor is always nonnegative because ek=(j,l) is oriented in this way when the convection-directed graph was constructed.
We illustrate the above results with the follow- ing example.
By examining the diagonal dominance for each column of the matrix M, we have G g diag([0, 0, 0, 0, 1]).(4.20) when ak b.It is clear that if the edge convection and the edge conductance scalars take these values Finally the generalized resistive network is illu- strated in the following Figure Using the graph representation of a convective M-matrix that corresponds to a stiffness matrix arising from the box method discretization, we restate the discrete curl-free condition (3.13) as Therefore, we have BA BA-rb 0 O.
On the other hand, suppose condition (4.21)   holds.We can solve the vector b in the follow- ing n n system, which can be made upper tri- angular by reordering the edges in the spanning tree edge set Ear: is the vector of the edge convection scalars as defined in Eq. (4.15).A convective M-matrix that satisfies the discrete curl-free condition is called a conservative convective M-matrix.Graph theory implies the following result.
THEOREM 4.7 The discrete curl-free condition (4.21) holds if and only if there exists an n-dimen- sional vector [1,   ,,]q-, such that Aq-b A, ( where A is the &cidence matrix for the convection- directed graph of M. Proof Suppose there exists an n-dimensional vec- tor b=[bl,...,bn] -r and Eq.(4.23) holds.From the definitions of matrices A and B and basic graph theory, we have BA-r BoA-+ A 0. COROLLARY 4.9 The convection-directed graph de- fined above for a conservative convective M-matrix M is acyclic if there is at least one convective edge ek, i.e., ;kk > O, in each fundamental cycle.Proof We represent a cycle in the graph G(M) by a vector c of length m in the following manner: if edge ek is in the cycle and is directed counterclockwise along the cycle; -1 if edge ek is in the cycle and is directed clockwise along the cycle; 0 otherwise.

(4.28)
Suppose there is a directed cycle in G(M).Then its representation c must be either nonnegative or nonpositive.Without loss of generality, we assume c is nonnegative.Since there is at least one con- vective edge in each fundamental cycle, there is at least one convective edge in the cycle that c represents, which implies that c-> 0. (4.29) On the other hand, using the graph theory, we have that there exists a vector f of length p, the number of fundamental cycles, consisting of O's and + l's, such that c -c fq-B.(4.30) Therefore, using Eq.(3.13), we have c q-A fq-BA 0, which contradicts (4.29), and hence the graph G(M) is acyclic, m COROLLARY 4.10 For a conservative convective M-matrix M, there exists a permutation P, such that the permuted matrix [mj,t, PMP -r has Imj, < Im;,, for j' < , i.e., the large absolute entries are all in the lower triangle.A similar statement holds for a e' (the backward permutation) such that [mj,t, P MP T has Imj,t,[ >_ Im [7,l   (4.32)   for j' < l'.Proof Since M is a conservative convective M- matrix, the convection-directed graph G(V, E) is acyclic by Corollary 4.9.Using basic graph the- ory, there exists a topological ordering (labelingi of the nodes in the convection-directed graph G, such that, after all the nodes in G labeled ac- cording to such ordering, an edge (j',I)EE implies that j'< '.Let permutation matrix P represent the relabeling of the nodes in G, such that Pej ej, and Pet et,.Since the direction of an edge is defined by Eq. (4.3), which is in- dependent of labeling, an edge (j',I')EE', or (j, l) E, is equivalent to [ml[ <_ Im0[, (4.33)   or Eq.(4.31).Therefore, Eq. (4.31) implies j < .
Since the graph G is acyclic, there also exists a topological ordering of the nodes, such that an edge (j', l') E' implies that j' > l'.Using this reordering, the permuted matrix has (4.31)   implying j> , or equivalently (4.32) implying j < .m Based on Corollary 4.10, we say that the per- mutation of the original stiffness matrix PMP -v is topologically ordered if and only if its upper- digraph, the directed graph corresponding to the upper-triangle of PMP -, is isomorphic to the convection-directed graph G or its transpose G-v.The transpose of a directed graph is a directed graph with the same nodes but reversed edge.Existence of topological orderings allows us to think of j < in (4.12) and (4.17), which is some- times convenient as in 5.1 below.
Recall that consistent orderings play a role in SOR theory [36].We give following relation be- tween topological orderings and consistent order- ings of a convective M-matrix M. THEOREM 4.11 If the convection-directed graph G of a convective M-matrix M is acyclic, then either all topologically ordered matrices PMP -v are con- sistently ordered or none of them are.Proof From the definition of topological order- ing, it is clear that the upper-digraphs for all topological ordered matrices are either isomorphic to G or to G-r.Since the cycle matrix of a directed graph is invariant under isomorphism, we can let B be the cycle matrix for the upper-digraphs isomorphic to G and B -r be the cycle matrix for those isomorphic to G q-. Since the edges of Gare exactly the reverse of those in G, we have B -v=-B.
If One of the topologically ordered matrix is consistently ordered, then we have B-re Be 0 (4.34)  according to the theorem by Rose [27].Applying the same theorem again, we have that all topologically ordered matrices are consistently ordered.
If one of the topologically ordered matrix is consistently ordered, then we have -Bq-e Be O. (4.35)   Applying the theorem in [27] in the opposite direction, we have that no topologically ordered matrices are consistently ordered.operator acting on a constant is always non- negative if and only if the divergence of the con- vection V./3 is nonnegative: V. (-aVl +/31) V./3 > O. (4.36) This leads us to define the following discrete ver- sion of the nonnegative divergence condition: CONDITION 4.12 (Discrete Nonnegative Diver- gence) dj 2e,nAoc > 0, (4.37) for <j < n.Here c=[cl,... ,Cm] -v, Ck=gkA:/2, is the vector of edge convection elements defined in (4.13).
.38)The notation (j,l)gj represents all (directed) edges (j,l) that are incident to node j..39) along the boundary of box By, with u the (out- ward) unit normal vector at the box boundary.In other words, we have clj (V lpj + O(h))Area(Bj), Assuming V./3 is continuous in each box, we can apply the intermediate value theorem and have V. da V./lp,Area(By) (4.41)   for some p' Bj.Furthermore, if V./3 has finite derivatives in Bj, we have X7./3lp,-V /3lp O(diameter(Bj)) O(h), (4.42) or / f V flda V /3lpArea(Bj) O(h)Area(Bj).

(4.44)
Using the same technique that was used to prove inequality (6.6) in Xu and Zikatanov [38], .45)In general, when the convection vector field in (1.1) fails to have nonnegative divergence, we may rewrite the equation with a scaled variable v such that u= r/v, r/ 0, and have fir/v) =f, (4.46)   and the discretization of (4.52) corresponds to the symmetric matrix S. The function e sym- metrizes the differential operator.
But we can also expand (4.50) as or V. (-&v +/v) f (4.47) with Theoretically, we can find v. (-v + ) _> o (4.48) so that the new convection field/3 has nonnega- tive divergence.In a discrete analogy, if we want to find a column scaling D for an M-matrix M such that MD is row diagonally dominant, we just need to find a vector d > 0 such that xx xx + xx u +/3xx =f" (4.53)When a > 0 and (d/dx)3 >_ 0, the sum of the first two terms, say 1 of (4.53) gives rise to a sym- metric and positive-definite operator which is perturbed by the first order operator /3(d/dx)u.
When the scaled convection b=a-l is curl- free, then we only need to choose r/= e .Similarly, when the convective M-matrix M satisfies the dis- crete curl-free condition, then we can choose D Dj -1 with Ds defined in (3.17).

CONVECTIVE ITERATION
In this section we discuss iterative methods for solving a discretized version of Eq. (1.1), Mu =f (5.1) with M=[m] defined in (2.43).We assume that the diseretized problem has at least one Diriehlet node, and therefore, M is irreducibly column dia- gonally dominant, since, otherwise, the matrix M on the left-hand-side of (5.1) is singular.More generally we consider iterative methods for solv- ing the system (5.1) when M is a convective M- matrix as defined in (4.12).
We assume that the reader is familiar with the basic theory of iterative methods of the form where M is split as M P-Q.For reference see Axelsson [1], Berman and Plemmons [7] and Varga [36].Recall that the rate of convergence of u k in (5.2) to u of (5.1) is determined by p(R) max I1, (5.3)   where R=p-1Q and Rv=Av; that is p(R) is the maximum modulus of the set of eigenvalues of R. The sequence u k converges to u if and only if p(R) < 1; the smaller that p(R) is, the better the convergence.Furthermore, R ln Pl (5.4)   is often called the asymptotic rate of convergence.
For any M E Z (Eq.(3.1)), M is an M-matrix if and only if any regular splitting(or weak-regular splitting) is convergent.For example, (P, Q) of (5.2) is a regular splitting, by definition, if p-1 > 0 and Q >_ 0. If M= P-Q is a regular split- ting of an M-matrix M, then we sometimes study the spectral radius of the matrix H= M-1Q which is related to p(R) as p(/4) p(R) + p(H); see Varga [36], Theorem 3.13.We will see that, under certain assumptions, it is possible to exploit convection to enhance con- vergence.By this statement we mean that usually p(Mc) < p(MD) where Mc corresponds to (1.1) with a nonzero convection (/3-0) and MD cor- responds to the same discretization but with fl 0 identically.To be more precise, we parameter- ize the scaled convection b=a-1/3 in (1.1) as tb for E[0, ).Therefore, the stiffness matrix M is parameterized by a single nonnegative number which scales the edge convection scalar, A, for all edge stamps, including the grounded ones as B(--tAk) --B(tAk)l (5.6) M(t)= gtc-B(-tAI) B(t) 1" ekEUEo Our analysis considerably extends, enhances, and corroborates the special case, asymptotic (/3--.oe) results in [37] by Wang and Xu.This is important since as /34 oe, the discretization parameter h must have h0 to insure that the discrete solution bears some resemblance to the true solution of the continuous problem (1.1).We show how the computational work W, amount of arithmetic operations, depends on the error in approximating the continuous solution by the discrete solution and the scaled convec- tion b ty-1/, i.e., we examine W(err, b; q, d) as defined in (5.103).See Theorem 5.13 and Corollaries 5.14 and 5.15.Our analysis requires some assumptions on the tensor field a and vector field /3 in (1.1).We re- quire the scaled convection to satisfy V b=0 (curl(b)=0) or, equivalently that b=-Vb for some scalar potential function b.Furthermore, we require that the discretized problem preserves a discrete analog of this condition; see 3.4 Condition 3.8.Using this assumption, we have shown that the convection-directed graph G(M) is acyclic (Corollary 4.9), and hence that there are (re-)orderings or permutations, P, such that the entries in the lower triangular part of PMP have special properties (Corollary 4.10).Interest- ingly and significantly the curl-flee condition is equivalent to the existence of a diagonal matrix D (see (3.17)) such that S=MD is symmetric (a Stieltjes matrix) as is (equivalently) Ms D-(1/2)SD-(1/2) D-(1/2)MDI/2); (5.7) see Theorem 3.10.Also the nonnegative diver- gence condition on/3 (4.37) reappears in our itera- tive analysis, and, when both V x b=0 and V.
/3>_ 0 hold, we obtain our most definitive result (Theorem 5.5).Using Eq. (3.19) we derive a bound on p(t) rela- tive to p(0) for the parameterization mentioned in (5.6), p being the spectral radius of the specific splitting under consideration.This bound is dependent only on the ordering of the purely diffusive problem corresponding to p(O).Hence, given the curl-free condition, topological orderings seem to play no role; we have discussed them in 4.2 since they may play a more important role when M can no longer be symmetrized as above.
Our discussion on iterative solving deals al- most exclusively with the Jacobi, Gauss-Seidel, and SOR point iterative splittings.We extend a block splitting example given by Wang and Xu [37].Some generalizations to block methods, preconditioning analysis, incomplete factorization [6], multigrid interaction and other advanced iter- ative techniques are immediate or straightforward while others will be quite subtle.Similarly, while our analysis does not always assume V./3_>0, weakening the X7 b--0 condition and carrying out the corresponding analyses and algorithmics are likely to be challenging.
as Ms DsMD-/1 (5.11)The symmetrized stiffness matrix is a tridiagonal matrix in the form of Ms TriDiag(-D, 2C,-D), (5.12) with D D(hb) defined in Appendix A. Now order- ing does not matter.Suppose the length of the domain is 1, then we have h 1In.Consider the Gauss-Seidel iteration matrix RGs P-IQs for Ms (and also for M), with Ps =1 TriDiag(-D, 2C, 0), Qs=-l TriDiag(0, 0, D) Notice that the iteration matrices for M and for Ms are similar to each other.The eigenvalues and the eigenvectors of matrix R6s can be computed exactly as 5.1.1.Toeplitz Tridiagonal Case Consider the box method discretization of Eq. (5.8) in one dimension.Suppose an equidistant mesh is used with (n-1) points between 0 and 1.
The hyperbolic function sech(x) can be viewed as the ratio between D(2x) and C(2x), which is also called the Bernoulli E-function (see Appendix A).
Remark 5.1 The eigenvalue problem for the above Gauss-Seidel iteration matrix can be converted to the eigenvalue problem for the Gauss-Seidel iteration matrix of the standard tridiagonal matrix Mo tridiag(-1,2, 1) using simple diagonal similarity transforms.
Let Mo Po-Qo and RGS,O PIQo denote the Gauss-Seidel iteration for Mo.Define and the spectral radius of the SOR iteration matrix using the optimal over-relaxation factor is PSOR + v/l-p j2 v/1 E2(hb) cos2(hTr) + V/1 E2(hb) cos2(hTr) (5.18) 5.1.2.Block Toeplitz Tridiagonal Case Next, we consider the box method with the use of an (n-1)x (n-1) square grid to the simplified model problem (5.8) in two dimensions with b [b, 0]a constant vector field.Notice that convec- tion only appears in horizontal edges, not in vertical edges.If we use a square grid as a mesh, then we obtain Ms Tc (R) I + I (R) To blocking.The n n tridiagonal matrices Tc and Tz are defined as Tc TriDiag(-D, 2C, -D), To TriDiag(-1,2, 1), respectively.Since the convergence rates for the Gauss-Seidel and the optimal SOR methods can be derived from pj using the consistent order- ing theorem, it is sufficient that we only provide the spectral radius of the Jacobi iteration ma- trix.For crosswind blocking, we have for the block iterative methods We can compare the asymptotic rate convergence R defined in (5.4) as h 0 as a function of b.For Jacobi, Gauss-Seidel and SOR R increases, i.e., b 2 (5.19) We again compare this with the convergence rate for corresponding purely diffusive problems, and have lim R'sR -- (5.20) (5.24) lim R'sr -- (5.25) 87r" h0 R,SOR However, the downwind blocking scheme is much worse that the crosswind blocking scheme as shown in Figure 5.1, because it takes almost no advantage of the convection.We find E(hb) cos(hTr) p(Rj) E(hb) + D(hb)(1 cos(hTr)) (5.26) E(hb) cos( E(hb) + O(hb)(1 cos(hTr)) (5.27) Rc,J Ro,GS R,SOR Our analysis here corroborates that of Wang and Xu [37].
See 3.4 and 4.2 for construction of the con- vection potential b from the edge convection A.
Notice that M(0) M(0) represents the stiffness matrix for a purely diffusive problem-the model problem (1.1) with /3=0 on the defining do- main 9t.Also notice that M(1) M is the original stiffness matrix.
We start with some results for systems that have nonnegative convection-divergence (see Condition 4.12), and then consider two further results which relax this condition.

General Analyses
In this section, we bound the iteration spectral radius of basic iterative methods for general con- servative convective M-matrix systems.We con- sider the parameterized edge assembly form of the parameterized stiffness matrix in Eq. (5.6).The symmetrized stiffness matrix (cf.(3.17) and 5.2.1.Nonnegative Convection-divergence Case We shall see that under the nonnegative diver- gence condition (Condition 4.12), the spectral radii of both the Jacobi and the Gauss-Seidel itera- tion matrices decrease monotonically if/3 in the model Eq.(1.1) does not vanish identically, and 2.5 2 ' "'%'*o,,.cos(nh) 0.9 "" Crosswind Blocking 25 ----.furthermore, we derive a bound on the spectral radii that generalizes our analysis for the special Toeplitz tridiagonal case.
Let Dd denote the diagonal matrix of the dis- crete convection-divergence of the original stiffness matrix M: Dd = Diag(dl,... ,dn) diag(Me); (5.31) cf Eq. (4.37).The discrete convection-divergence for the parameterized stiffness matrix M(t) is Da(t) diag(M(t)e) tDa, (5.32)   and the edge assembly of M(t) in Eq. ( 5.29) can be equivalently written as Ms(t) ---Dd + Z gk ek EEtOEo C(tAk) -D(tAk) l -D(tA) C(tAk) k (5.33)where the identity (A.5) is used.We further consider the edge assembly forms of the Jacobi and Gauss-Seidel splitting matrices and their deri- vatives with respect to the parameter t.Based on the form of Ms in (5.33), we have b(tA) D(tA)(1 C(tA)), (5.39)   the derivative of the above splitting matrices can be written as M satisfies the nonnegative discrete convection- divergence condition (4.37), or Da >_ O, and has at least one convective edge, say Ak:fi0 for some eke E U Eo, then the spectral radii of the Jacobi and the Gauss-Seidel iterative methods for the parameterized system M(t)uh=fh decrease mono- tonically as > O.
Proof First notice that M(t) and Ms(t) are simi- lar, and therefore, their Jacobi (Gauss-Seidel) iteration matrices have the same eigenvalues.
Next, we consider the derivative of the iteration matrix as d -1 /s,J,GS (Ps,J,GS Qs,J,GS) ----Ps,J,GsPs,J,GS s,J,GS--Ps,J,GSOS,J,GS. (5.44)The derivative of its spectral radius can be written as --[-Ys'J'GsRs'J'GSXs'J'GS (5.45)   /J,GS T Ys,J,GSXS,J,GS where Xs,J,GS and Ys,J,GS are the right and left eigenvectors for PJ,GS respectively for Rj,GS.See Appendix B for detailed derivation of iS.From (5.44) and (5.45), we have Yf,J,GsPj1,GS ps,J,Gsxs,J,Gs /J,GS _<-T PJ,GS.(5.46)   Ys,J,GSXs,J,GS According to the Perron-Frobenius Theorem [36], Theorem 2.1 and [36], Exercise 6 on page 75, both xs,j,GS and Ys,J,GS are vectors with all positive entries.Since Ps,J,GS is an irreducibly diagonally dominant M-matrix, p-1 is nonnegative and s,J,GS has all positive diagonal entries.Since there is at least one convective edge in G(M) and Da>_O, .-1 the matrix Ps,J,GS is nonnegative and has at least one positive diagonal entry when > 0; see Eq. (5.40) or (5.42).For > 0, we see that /J,GS(t) < 0, (5.47)   and therefore PJ,GS(t) decreases strictly monotoni- cally, m Next, we strengthen the condition on convective edges such that each node is connected by at least one convective edges, or A* > 0 with A* defined as A* =_ min max jt.
(5.48) <<_n (j,l)gj (j,I) EEUEo Proof The proof is complete after we separate vari- ables and integrate both sides in the above inequality, m Remark 5.4 Using Theorem 5.7 and Remark 5.9 below, the bound on the spectral radius PJ,GS can be improved as pJ,GS(t) _< min(pj,Gs(O)C(*t) -r, Cj,Gs/t) (5.52)   for some constant Cj,GS when A* > 0.
Eqs. (5.42) and (5.43), we have Similarly, from 1(D2(IA[min)) /bs'GS(t) --> 7 1- C([/[min) es,GS(t), (5.56) )s,s(t) < 1-C(tlAlmin Qs,6s(t). (5.57) Applying the form of s,J,GS in Eq. ( 5.44), we have for both Jacobi and Gauss-Seidel that s,J,GS ---7 1-C(iA[mint) -+-C([Almint)t 1] RsJGS'' Furthermore, we apply the expression for the derivative of the spectral radius of the iteration matrix in Eq. (5.45) and have (5.39),we are able to separate variables (p and t) on both sides on inequality (5.59), then integrate both sides, and finally complete the proof, m Notice that the bound given by the above theo- rem decreases like the Bernoulli E-function if and only if IAImin > 0, or all edges in the correspond- ing convection-directed graph are convective.If the convection vector field /3 does not vanish anywhere in the defining domain f, then a mesh can be constructed (adjusted) such that no mesh edges are orthogonal to the convection vector field, and hence the discretized system will meet this condition.
Further notice that the spectral radius bound for the Gauss-Seidel iteration method is not as tight as that for the Jacobi iteration method because the positive off-diagonal entries in P6s were simply discarded in the above proof.In the best situation, where the consistent orderings exist and are used, as we have seen in the Toeplitz and block Toeplitz cases, the spectral radius of the Gauss-Seidel itera- tion matrix decreases as E2(At).In more general situations, one may use the Jacobi bound provided by Theorem 5.5 in combination of Theorem 4.8 in Varga [36] to derive possibly a tighter bound for the Gauss-Seidel case.

No Assumption on Divergence
We start by proving the bounds on the 2-norm of the structural matrix of an undirected graph G(V, E).Let Q0 be the strictly upper-triangular matrix such that if(j,l) EE j<l, (5.60) Qo(j, I) 0 otherwise.
Then we assign directions to the edges in G to make it the upper-digraph of Q-+ Q0, i.e., assign edge (.j,l) the direction from j to if j < l; otherwise, assign the direction from to j.Let Outj and Inj denote the set of nodes that are connected with j by its out-edges and in-edges respectively.
(5.65) /6Outj Notice that each node in Out] has at least one in- edge (j, l).Therefore, all entries of (QQo) are non- negative, and the jth row sum can be written as and f C(At)-D(t) ifek isanon-groundededge, C(It)/2 ifek isagroundededge. (5.70) is yet another member of the family of Bernoulli functions (see Appendix A).The constants Cs and Cas depend only on the mesh.
Proof Let Hj and Hs denote MTI(Q + Qs) and Mj-1Qs respectively.We have p(HJ) [IHjII2 _< IIMj-IIzlIQ /QII2.(5.71) We will give bounds on IIMj-IIz and IIQ / QII2 respectively as follows.(5.66) V J The second bound of the lemma is also proved by applying the Gershgorin Theorem.m We show in the following that as the param- eter is sufficiently large, the convergence rate de- creases as the magnitude of convection decreases uniformly in the domain.Here, we still assume curl-free condition on both the continuous and the discrete convection.

PJ + p(Hj)
Applying the second part of Lemma 5.6, we obtain for the Gauss-Seidel iteration that v/lOutmaxl IInmaxl cj CGS iOutmaxl + ilnmax (5.79)   which certainly can be controlled by the mesh.The proof is complete by using the fact that both Jacobi and Gauss-Seidel splittings and regular splittings, and hence pj and PGS _< 1. m Remark 5.8 When t[/k]min>> 1,M tends to a diagonal matrix due to the Bernoulli D-function, and the condition number is roughly bounded by the ratio between the largest and the smallest diagonal elements of M.
(M) IOutmaxllAlmax (5.80) ,, with A* defined in (5.48).Remark 5.9 If we assume that all edge convection scalars uniformly bounded away from zero" min > 0, (5.81) then for E [0, oe), the spectral radius of the Jacobi or Gauss-Seidel iteration matrix for M(t) decays exponentially as oc.This follows easily from the properties of the Bernoulli functions.Even if there are some zero-convection edges but A* > 0, the bound given by the above theorem still decays as O(1/t) according to the property of function F. We next study the behavior of the spectral ra- dius of Jacobi or Gauss-Seidel related iteration matrices Hj and Has as t---,0.First, we consi- der the derivative of the Gauss-Seidel spectral ra- dius HGs at t=0.Using (5.33) and the fact that  (5.8) It is clear that when the discrete convection-diver- gence dj. is nonnegative for < j < n and nonzero at least at one node, then iS(Has) < 0. However, when dj. are negative for all < j < n, then p(HGs) increases as increases from 0. Since, we know the spectral radius will eventually reach the ex- ponential decay region when is sufficiently large, there must exists a positive number to such that p(t) < p(0) for all > to, since p(t) is a continuous function.We will give a bound on the value of to which will depend on the divergence of convection, the condition number of the stiffness matrix, and the mesh properties.We first note the following lemma; the proof is straightforward.
Then it suffices to take > Cm (-V'/3[)ma x + O(h), where Cm 6C2/(gmin COS2(0a*-(h*,/3*))) is a quantity that can be controlled by the mesh.Of course, if the (-V./3[)max turns out to be suf- ficiently negative, we only need to take > 0. m Notice that the value of to largely depends on the characters of the diffusion and the convec- tion coefficients.For V./3 > 0, to 0. The value 7 can be controlled by meshing techniques so it should be close to 2. Obviously, when the non- negativity of the convection-divergence is given, then the spectral radius of the Gauss-Seidel iteration matrix decreases monotonically for > 0.
Furthermore, we have the following analysis for /-/c;s and tb(Hc;s).
Remark 5.12 These bounds we obtained in Theorem 5.5 and Theorem 5.7 can be applied to the SOR methods in a similar manner using Ob-< PSOR < V/ab- See Varga [36], Theorem 4.9.For example, the bounds on PSOR using Eq.(5.53) in Theorem 5.5 can be derived as (5.97)

Operation Count
We count the total number of operations, such as additions and multiplications, that are needed for (iteratively) solving linear systems arising from discretization of convection-diffusion equations.
It is clear that the total number of operations is the product of the number of iterations and the work per iteration, or W-mk. (5.98) The number of iterations k is roughly the recipro- cal of the asymptotic convergence rate Ro In p for given tolerance, i.e., In err (5.99) The number of the operations per iteration is proportional to the number of nonzeros in the stiffness matrix, which, in turn, is proportional to the number of edges in the mesh, i.e., rn oc h -d (5.100) with d the spatial dimensionality.
Since the discretization error depends on both the mesh size h and the scaled convection b (cf. [31], Appendix B), it is not appropriate to fix the mesh size while solving a spectrum of problems with different scaled convection.Our analysis that follows is with respect to a given tolerance that bounds both the discretization error and the itera- tive solving error.We first assume that the discreti- zation error has the form of err hq(C1 -+cllb[I), (5.101) for some positive constant numbers C1, C2 and q.We further assume that q > 1.Then in order to ensure that the discretization error is less than the tolerance we need to set the mesh size as err ) (l/q) h= Cl + C2llbl] (5.102)Therefore, we can express the total amount of work W in terms of err, b, q and d as [lnerr ( Cl -+-C2b ) (d/q) W(err, b; q, d) [In p(err, b; q)l err (5.103)  while treating h as an implicit variable.If we consider the work for solving a convective problem relative to solving a corresponding diffusive problem, we have W(err, lib[l; q, d) [ln p(err, 0; q)[ W(err, O; q, d) [In p(err, ( C1) where we assume that [A [min Ch h b for some constant Ch.Notice that both terms on the right-hand-side of the above equation are always between zero and one.We pick an arbitrary num- ber b0 > O.When 0 _< b <bo, we have p(h(err, llbll; q), b) _< p(h(err, llbll; q), 0) where It is obvious that the convective-diffusive ratio of the work is the product of the ratio of the num- ber of iterations k= Ck/[ln(p) and the ratio of the number of operations per iteration m.Given (5.102), the increase of m(b)/m(O) is inevitable as b increases.The main advantage of using iterative methods for solving convective problems is stated in the following theorem and corollaries.Here, we only consider the Jacobi and the Gauss-Seidel iterative methods since similar results can be obtained for the SOR method using Remark 5.12.
THEOREM 5.13 Assume the conditions in Theorem 5.5.The ratio of number of iterations k (b) / k (0) is bounded above by a number which does not depend on the spatial dimensionality when q >_ in (5.101).
Proof Since k(b) ln p(err, 0; q)l k(0) lnp(err, llbll;q) (5.105) and. the numerator depends only on err and q, it is sufficient for us to show that p(err, b [[; q) is uniformly bounded above by one for given err and q.
We see that the iteration spectral radius p(h, b) depends on both the mesh size and the magni- tude of convection.Using the bound given in Theorem 5.5, we have err ) (l/q) hi C1 + C2bo (err) (l/q) h2 l It is clear that the bound maXhl<h<_h2P(h,O) is bounded above by one, and it only depends on the choice of b0 and the iteration spectral radius of the diffusive problem as well as err and q.
When b _> b0, we have p(h(err, Ilbll; q),b) < E(Chh(err, Ilbll; q)llbll) < E(Chh(err, b0; q)bo) (5.106) since h(err, b II;q)II b monotonically increases as bl[ increases for q > 1.From the property of the Bernoulli-E function, it is clear that E(Chh(err, bo;q)bo) is strictly less than one, and it only depends on b0, err, and q.Since the spec- tral radius (p) does not depend on the spatial dimensionality, it is clear that the bound does not either, m Examining (5.106) in view of the properties of the Bernoulli-E function gives a more precise re- sult as follows.

CONCLUDING REMARKS
In this section, we comment briefly on some aspects of our work which flow beyond the mainstream of our results.Remark 6.1 The constant-j box method has been successfully extended to the investigation of wavefront propagation in anisotropic cardiac tis- sue; see [33].By using a triangulated mesh that is formed with consideration of the underlying j -aVu (6.1)   as any variation of ion concentration inside or out- side the cell is considered to be negligibly small.
For diseased myocardium, it is possible to have local accumulation of ions near an injury leading to spatial gradients in ion concentrations.These concentration variations will also contribute to the overall ion flux and will need to be explicitly taken into account in the model.The advantage of using the box scheme described in this work is that any reformulation of the model equations to handle both diffusion and drift of ions can be incorporated in a straightforward manner with no loss in computational performance.
The drift-diffusion equations arise more natu- rally when modeling ion fluxes at the indivi- dual channel level; see Eisenberg [13], and Gardner, Jerome and Eisenberg [15].In conven- tional heart modeling, the channel currents are lumped to create a macroscopic decryption of the current kinetics.There are, however, a number of pharmacological or even genetic therapies for wavefront anomalies (arrhythmias) that will need to target at the level of the ion channel.The effects of these, drugs on the overall macroscopic mod- els may be hard to predict unless the models themselves are initially constructed from a chan- nel perspective.This construction process may in- volve solving both the Poisson and drift-diffusion equations in three dimensions over the relevant D.J. ROSE et al. portion of a cell.The governing equations are analogous to those used in modeling semiconduc- tor devices.The results from this work suggest a robust numerical approach that will enable the design of a drug or gene therapy and its action on a given cell much in the same way that one might design a specific portion of a semiconduc- tor element.This "re-engineering" of the cardiac myocyte may lead to the equivalent of creating a semiconductor heart.Remark 6.2 We note that the constant-j box method discretization paradigm presented here is not limited to triangular meshes.Indeed, the edge-pair constant-j assumption is easily seen to apply to quadrilateral meshes, for example.Even the "terminating line" quadrilateral meshes as in Selberherr [30], page 177, present no fundamental difficulties.For example, once the box is chosen as in Figure 6.1 (adapted from [30], Fig. 6.2-3), edge-pairs are then chosen to determine how to compute j on various parts of the box boundary.
(Some of the edges of these edge-pairs, e.g., (0, 4) and (0,5), may not be part of the original quadrilateral mesh.)This provides an alternative to using finite difference and allows flexibility in both the choice of the box for an unknown and the corresponding edge-pairs.Extension to 3D meshes is also straightforward in principle.Of course, the special properties of the stiffness matrix are no longer guaranteed.Remark 6.3 Formally, extending the constant-j box method to three dimensions is straightforward where the two-dimensional edge-pairs are replaced by edge-triples (Fig. 6.2) which will be used to derive the local constant-j's.In three di- mensions, box vertices need to be specified for all edges, triangles, and tetrahedra in the mesh.The edge conductance values will depend on local tetrahedron vertices and box vertices, as well as evaluation of the 3 x 3 diffusion tensor.Unlike the two-dimensional case, the box-vertex dependency does not naturally vanish in the expres- sion for edge conductance.Finding a clean mesh condition for the M-matrix discretization property is unresolved when diffusion is anisotropic and varies over space.As in 2D, these conditions may be needed to serve as guidelines for mesh generation in three dimensions.However, in an isotropic case, i.e., c a(x, y, z)I3 x 3, the M-matrix property of the stiffness matrix can be obtained when the mesh is a Delaunay tetrahedralization.Then we have Ax px. (B.1) Taking the derivative with respect to on both sides, we have ]x + ASc pSc + x, (B.2) or y-)x + y-ASc y-pSc + y-C x.
Since y is the left Perron eigenvector of A, we have y-VA=y-p.Therefore, we have y-4x [9-y-x" (B.4)

FIGURE 2 . 1
FIGURE 2.1 Box B; for mesh point Pi as an interior point (left) and as boundary point (right).
2) and the discrete curl-free condition (3.13) are satis- fied.The values of b are obtained through either of the methods described above.Let us define a positive diagonal matrix D as e

4. 1 .
Graph and Circuit Interpretation of Convective M-matricesWe define convective M-matrices as m-matrices which are irreducible, column diagonally domi- nant, and structurally symmetric.Let M=[mt]

Remark 4 . 4
By defining the edge diffusion and the cycle matrix of the graph G repre- senting M

(4. 25 )
Then by applying(4.25) in the first step, (4.26) in the second step, and (4.21) and (4.5) in the third step, have shown that Eq. (4.23) holds with the b obtained from(4.26).m Remark 4.8 The proof of Theorem 4.7 provides an algorithm to construct the vector b if the edge convection vector A is given.The algorithm is simp- ly backsolving Eq. (4.26), setting %= 0.

1 4. 3 .
Nonnegative Divergence and Row Diagonal Dominance At the continuous level, the convection-diffusion
(5.21)   if we use the crosswind blocking(vertical ordering)    or M To (R) + (R) Tc(5.22)   if we use the downwind(horizontal ordering)

2
(t)  =-Dd + Z gk --D(tAk) If the convection-directed graph G(M) for a given conservative convective M-matrix

Figure A. 1
for a plot of the Bernoulli E-function.

FIGURE 5 . 2
FIGURE 5.2 Upper bounds on convective-diffusive ratio of the numbers of iterations k, R ln(p)[.

FIGURE 6 . 1 FIGURE 6 . 2
FIGURE 6.1 Terminating line mesh point with two standard edge-pairs (and box segments) and three non-standard edge-pairs (and box segments).

FIGURE A. 3
FIGURE A.3 The Ski-Slope function (IB(z)t,z E C).
It is easy to verify that