Gridding and Discretization For Divergence Form ( Semiconductor-Like ) PDEs

We develop the box method for solving partial differential equations in the divergence (or conservation law) form. We first use graphic theoretical approaches to generate grids, i.e. partition the governing domain into boxes. Then we apply Green’s theorem box-wisely and the constant-j assumption edge-pair-wisely to obtain the discretization system of equations, which lead to the numerical solution to the problem. We show that the box method is inherently more efficient than the traditional finite element method for linear (or convection-diffusion) problems in and 2 dimensions. We also present some potential thoughts on implementing the box method for problems in higher dimensions and with nonlinearity.


INTRODUCTION
In this paper we will discuss gridding and discretiza- tion for a PDE (or PDE system) of the form.
We are mainly interested in the cases where has a dominant convection term and the region on which the PDE is defined has complex geometry.Rather than the traditional finite element methods, we use two graphs for partitioning (gridding) the region, namely, the Voronoi diagram and the Delaunay trian- gulation, which are straight line duals from the geometric point of view.
Our approach to discretization (in 2D) will be to use the "box method" (finite volume and finite box methods [1, 5]) using Voronoi boxes for appropriate grid points in the region specifying the PDE.By using the divergence theorem, this transforms (1.1) to the operational equation fo j.nds+ f f -O, (1.5)   Bi where B is a box containing grid point p/with bound- ary OB i.After approximating the integrals in (1.5) on each box B associated with the unknown u i, we obtain the coupled system of nonlinear (or linear) algebraic equations representing the discretization.We have found that discretization and gridding are tightly coupled to each other in a way that discretiza- tion imposes a set of specifications onto gridding while limitations of gridding schemes may require adjustments on discretization strategies, it is important to understand how the trade-offs are made among the generality of the problem, the numerical quality of the discretization, and the complexity of grid con- struction.
In order to discretize (1.5), we need to approximate the integrals.We will derive a discretization based on the constant j assumption, namely, that j is constant pairwise on the edges of the Delaunay triangles dual to the sides of box B i. Suppose that B i, B 1, and B are neighbors pairwisely, and therefore PTP] and PTP form an edge pair.Depending on whether the circumcenter of Applp k lies inside or outside the triangle, we need to approximate the integrals in (1..5) in different ways (Figure 1).
In particular, there may be more errors in case 2 of Figure because we approximate integrals outside the edge pair by using the interior information.Let J [kil].andF[kil denote approximations of the part of the line integral of j and the area integral off in (1.5)   that is contained in the edge pair (p/pl,pg)).In the case where o and 13 are constant, we have the follow- ing result J[kil] --Iltzllnj + IItllnirj (" +" for case and "-" for case 2)  hlhi Ihilll)   ("+" for case and "-" for case 2) H-(hil hik)T n-h Ilhll hrct-1 , and B() ,/(e -1).If z and [3 vary smoothly over space, we can evaluate them at the circumcenter of Ap/pl (or the edge midpoints), and the above form can still be applied to achieve first order accuracy.This allows us to assemble for each box the parts of the line integral and the area integral which, in turn, relate the variable U to the other u corresponding to Delaunay neighbors, say p, of the box mesh point Pi.We see that the constant j discretization is a natural generalization of the Scharfetter-Gummel [1, 5, 6]  discretization used in semiconductor modeling.
When j is linear, as in (1.4), the constant j discreti- zation produces a linear system to solve with interest- ing matrix properties which are closely related to those of the gridding output.For instance, a lower bound on the angles of the Delaunay triangulation gives a lower bound on the conditional number of the matrix, and an upper bound on the angles is crucial to the accuracy of approximation in (2.2).These consid- erations motivate the following gridding scheme.

GRIDDING
Suppose the region on which the PDE in (1.1) is speci- fied is a connected open set E2 of 2. We partition E via the Voronoi scheme based on a set of grid points

P {Pl
Pn}" For each grid point Pi, we find its Voronoi box B where each point is closer to Pi than to any other grid points in P. The collection of the Voronoi boxes, v {B Bn}, is called the Voronoi diagram of P. In addition, the straight-line dual planar graph of the Voronoi diagram yields the Delaunay tri- angulation D of P, with the grid points in P being the vertices of the triangles.The edges of the Delaunay tri- angles are perpendicularly bisected by the sides of the Voronoi boxes.The union of the Voronoi diagram and the Delaunay triangulation yields a grid on D.
Most of the previous results on grid generation have been oriented to the finite element method.Two schemes that have been widely studied are quadtree based techniques [2] and Delaunay based techniques [3, 4].Our grid however has to meet a more specific set of conditions that are required by our discretiza- tion scheme.
1.The angles of all the triangles of D must be bounded both below and above by constants 0min and 0max respectively.The constant 0min must be large enough to ensure that a well conditioned algebraic system can be obtained from the discre- tization; 0max should be small (ideally < 90) so that the effect of case 2 in Figure is negligible.
2. The edges of all the triangles of D must be bounded above by h 0, so that the discretization can reach certain required accuracy.

No Voronoi boxes in v cross
We call a grid that satisfies these conditions a quasi- uniform grid.Based on some of the previous work [2,   3, 4], we have implemented our own scheme for gen- erating quasiuniform grids on a region with general boundary geometries.Here is a brief description of our algorithm: 1. Decompose the given region f into subregions by using the quad-tree based techniques [2].2. Generate points in each subregion according to its relative position to the boundary of D. 3. Apply the filtering process of Miller et al [3] so that the Delaunay triangulation of the remaining points meets condition above.4. Generate the Voronoi diagram and the Delaunay triangulation for the remaining points.
Let P be the set of output grid points of this algorithm.We can show that the size of P, IPI, is optimal, i.e., P 1< c h (3.1) for some constant c > 0, where S(E) denotes the area of f2.Furthermore, the points of P are uniformly dis- tributed on 2 with density approximately h 2 In the cases where and [, or the source density function, f, are ill behaved at some points, we can locally refine the grid so that these points will be included in P and conditions 1-3 above for the grid will still hold.
4. EXAMPLES 5. CONCLUSION Here we show some of our results on solving two simple cases of (1.1) with j oVu + [u.The first case (see Figure 2) is an extended 1D problem where f is specified as a square region with vertices located at (_+1, _+1).Here o is a 2 x 2 identity matrix, [ =( 10] andf--0" The bundary cnditins are0 all Dirichlet and fit the function p(x,y) which can be solved exactly from In the second case (see Figure 3), we solve the Poisson equation VZuon the disk (x, y)lx 2 + y2 < }.The boundary conditions are Dirichlet and identically zero.We see that our results fit well to the analytical solu- tions for both cases.For the dominant convection case, our box method gives a much better approximation to the analytical solution than the finite difference or the finite element methods.We have shown that the constant j discretization is efficient and convenient for solving the 2D boundary value problems with j linear in terms of Vu and u.We are extending our study to the cases where j is nonlin- ear and generalizing our 2D schemes to 3D.