- Email: [email protected]
Magnetostatic solution by hybrid technique and fast multipole method
ARTICLE IN PRESS
Physica B 403 (2008) 368–371 www.elsevier.com/locate/physb
Magnetostatic solution by hybrid technique and fast multipole method G. Gruossoa, M. Repettob, a
Politecnico di Milano, Dipartimento di Elettronica e Informazione, I-20133 Milano, Italy Politecnico di Torino, Dipartimento di Ingegneria Elettrica, C.so Duca Abruzzi 24, I-10129 Torino, Italy
b
Abstract The use of fast multipole method (FMM) in the solution of a magnetostatic problem is presented. The magnetostatic solution strategy is based on finite formulation of electromagnetic field coupled with an integral formulation for the definition of boundary conditions on the external surface of the unstructured mesh. Due to the hypothesis of micromagnetic problem, the resulting matrix structure is sparse and integral terms are only on the RHS. Magnetic surface charge is used as source of these integral terms and is localized on the faces between tetrahedra. The computation of the integral terms can be performed by analytical formulas for the near field contributes and by FMM for far field ones. r 2007 Elsevier B.V. All rights reserved. Keywords: Micromagnetics; Fast multipole method; Cell method; Hybrid techniques
1. Hybrid magnetostatic solution The use of an hybrid form of magnetostatic computation based on finite formulation of electromagnetic fields (FFEF), more often called cell method, and on a Green function applied to magnetization sources, has been presented [1,2]. In this formulation an integral boundary condition on the exterior surface of the mesh is related to the magnetization sources. Integral boundary conditions are used to terminate Ampere’s law for all boundary edges expressing external magnetic field circulation as a potential difference evaluated by the Green function integrals at the center of each boundary face. Under the micromagnetic hypothesis, the magnetized region has permeability m0 and, for each time step of the dynamic analysis, magnetization distribution is supposed to be known. The absence of magnetic polarization, that is of a magnetization-dependent field like M ¼ wH, where w is the magnetic susceptibility, eliminates magnetization contribution depending on configuration variables like magnetic vector potential. Thus all integral terms are known and can be carried on to the RHS. Another consequence of this fact is that there are no integral relations between matrix rows Corresponding author. Tel.: +39 011 5647140.
E-mail address: [email protected] (M. Repetto). 0921-4526/$ - see front matter r 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.physb.2007.08.051
and thus matrix structure remains sparse as in the usual formulation of the cell method. This numerical approach has been already applied to demagnetizing field computation on a structured hexahedral grid [3] and is now applied to a tetrahedral mesh. Magnetization is considered to be uniform in each tetrahedron so that magnetic sources are located only on triangular faces. The computation of magnetizationdependent integrals is performed by fully analytical formulas, nevertheless this phase is one of the most timeconsuming parts of each magnetostatic solution. From these considerations stem the idea of using fast multipole method (FMM) to speed up the integral boundary conditions computation.
2. Fast multipole scheme for RHS computation 2.1. Multipole expansion Magnetic scalar potential is computed by an integral operator applied to magnetization-dependent sources. The scalar potential c is expressed by Z cðrÞ ¼ rM Gðr; r0 Þ dOM , (1) O
ARTICLE IN PRESS G. Gruosso, M. Repetto / Physica B 403 (2008) 368–371
where G ¼ 1=ð4pjr r0 jÞ is the usual Green’s function of three-dimensional static problems and rM is the magnetic charge density obtained as rM ¼ rMrSM, and OM is the volume domain containing the magnetization. This term can be evaluated by considering two contributions to the potential: one due to near field and one due to the far field c ¼ cnear+cfar. The first term [2] can be calculated by means of Eq. (1), while the last one can be obtained by resorting to an expansion in multipoles [4]: cfar ðr; W; jÞ ¼
l X l X Mm l Y m ðW; jÞ, rnþ1 l l¼0 m¼l
369
cube well separated from them. The cost of this algorithm is approximately Nf log Ns, where Nf is the number of evaluation points and Ns the number of the sources. Generally this approach is not enough to improve the computational efficiency and multilevel algorithms are introduced [7]; in the present formulation, the number of field points, located only on the surface of the active volumes, is much lesser than the number of sources located in the whole active volumes, therefore only the O(N log N) algorithm has been considered.
(2) 2.3. Data structures
with
where ri, Wi, and ji are spherical coordinates of ith charge qi ¼ siAi, and r, W, and j are the coordinates of field point. Ym l are the normalized spherical harmonics of degree l and order m given by the following: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2l þ 1 ðl mÞ! m m P ðcos WÞejmj , Y l ðy; jÞ ¼ (4) 4p ðl þ mÞ! l
In this section a data structure based on octree [8] subdivision of the domain combined with a bit interleaving techniques [9] will be described. This technique allows an efficient indexing of the cube in the hierarchical tree and then a robust and fast method to find relationship between cubes at each level of the octree. Briefly some characteristic aspects of the method will be outlined: considering a box at level l the whole number of children of this box is 8, numbered from 0 to 7, therefore the index of one of these boxes is unique and it is given by the following relationship:
with the condition
IDXðn; lÞ ¼ ðN 1 ; N 2 ; . . . ; N j Þ,
Mm l ¼
N X
qi rli Y m l ðWi ; ji Þ,
(3)
i¼1
m mn Y m l ðy; jÞ ¼ ð1Þ Y l ðy; jÞ,
(5)
where Pm l ðcos WÞ are the Legendre polynomials of degree l and order m [4]. 2.2. O(N log(N)) fast multipole scheme Several fast multipole schemes have been developed in order to increase the computational efficiency of the method [5–8]. All these methods are based on the same basic principle: in order to make systematic use of multipole expansions, a hierarchy of boxes which refines the computational domain into smaller and smaller regions will be defined. At refinement level 0, the whole computational domain will be considered. Refinement level k+1 is obtained recursively from level k by subdivision of each box into eight equal parts. This yields a natural tree structure, where the eight boxes at the k+1 level, obtained by subdivision of a box at level k, are considered its children. By starting from the generated tree it is possible, for all sources which are ‘‘sufficiently far’’ from a given field point, using the multipole expansion, to summarize their total interaction. According to error estimation, a convenient definition of ‘‘sufficiently far’’ (or ‘‘well separated’’) is a relationship between the diagonal of each cube and the distance ofp field ffiffiffiffiffiffi point from the center of the considered cube: r0 o2 ð3Þw, where r0 is the distance of field point from the local origin and w the length of the edge of the cube. The first algorithm consists in evaluating the far field at all field points by considering the contribution of each
(6)
where l is the considered level, Nj the index of the box at level j given by Nj ¼ 0, y, 7 with j ¼ 1, y, l and n is the global number of the cubes at level l with n ¼ N1(8)(l1)+ N2(8)(l2)+?+Nl1(8)+Nl. For example, by referring to the construction shown in Fig. 1,the cube N1 ¼ 0 and N2 ¼ 5 at level 2 is the cube n ¼ 5 of the level 2 and so on; in this way the uniqueness of the indexing is guaranteed. In this way it is possible to obtain the father box (PRT) for each pair (n, l)by means of PRTðn; l 1Þ ¼ ðN 1 ; N 2 ; . . . ; N j0 Þ,
(7)
with j0 ¼ 1, y, l1. From the previous statement it is possible to compute all the children (CHD) at level l+1: CHDðn; l þ 1Þ ¼ ðN 1 ; N 2 ; . . . ; N jþ1 Þ,
(8)
where Nj+1 ¼ 0, y, 7. In other words the use of octree makes obtaining parent and children indices very convenient. Indeed the above operations are nothing but shift operations in the bit representation of n. Performing a right
Fig. 1. Octree subdivision of the domain and cube hierarchy.
ARTICLE IN PRESS G. Gruosso, M. Repetto / Physica B 403 (2008) 368–371
370
bit-shift operation on n by 3 bits, one can obtain the index of the parent. One can list all indices of the children boxes of n by a left bit-shift operation on n by 3 bits and adding all possible combinations of 3 bits. In order to decode this structure in bit representation, it is necessary to introduce a spatial ordering. Firstly it is possible to define a proper scaling of the point Pj ¼ (xj, yj, zj) as ¯ j ¼ Pj Pmin , P D
3. Results The results, obtained on a cylindrical domain with an impressed uniform magnetization, are outlined only in terms of cpu-time comparison of one magnetostatic solution obtained with three different schemes:
(9)
where Pmin ¼ (xmin, ymin, zmin) is the lower corner of the box enclosing all the domain and D is the length of biggest edge of the box. In this way the system will be scaled into a box with the lower corner in the origin and with edges equal to 1. ¯ j ¼ ðx¯ j ; y¯ j ; z¯ j Þ can The scaled coordinate of the jth point P be represented in binary form as
LU factorization of the matrix and computation of RHS by exact formulas; GMRES with incomplete LU preconditioning and computation of RHS by exact formulas; same GMRES with FMM computation of RHS.
As can be seen from Fig. 2, the combined use of GMRES and FMM sensibly reduce the computational time of one solution. Despite the O(n3) computational burden of LU factorization the O(n) cost of back-substitution can make
x¯ j ¼ ð0:a1 a2 a3 . . . Þ2 , y¯ j ¼ ð0:b1 b2 b3 . . . Þ2 , z¯ j ¼ ð0:c1 c2 c3 . . . Þ2 , where ai ¼ 0, 1, bi ¼ 0, 1, and ci ¼ 0, 1. Instead of having 3 indices for each point it is possible to form a single binary index that represents the same point by means of ordered interleaving of the digits (bit interleaving), so the index becomes P ¼ ð0:a1 b1 c1 a2 b2 c2 a3 b3 c3 . . . Þ2 ¼ ð0:N 1 N 2 . . . N j Þ2 ,
(10)
where Nj ¼ (aj bj cj)2. By starting from these considerations, it is easy to perform the following operations: 1. Finding the IDX containing a given point. For a generic scaled point obtained from Eq. (10) at level l, by referring to a hierarchical structure ordered as in Fig. 1, the following relationship is true: P¯ ¼2 BoxððN 1 N 2 . . . N l Þ8 Þ where the value of the Box function is given by Eq. (6). 2. Finding the center of a given point. The relation between the coordinate of a point and its box index enables easy finding of the coordinates of the center for a given box index n at level l, simply converting this information in binary form, using Eq. (6) and then decomposing the result into the three coordinates:
Fig. 2. Cpu time comparison by using different algorithms for the calculation of the RHS.
n ¼ ða1 b1 c1 a2 b2 c2 . . . al bl cl Þ2 , x ¼ ð0:a1 a2 . . . al Þ2 , y ¼ ð0:b1 b2 . . . bl Þ2 , z ¼ ð0:c1 c2 . . . cl Þ2 . Combining this algorithm with the multipole calculation it is possible to reduce the computational effort due to the proposed formulation.
Fig. 3. Cpu time cost of all solution phases: int, computation of integral coeffcients; mat, matrix computation and assembling; sol, linear system solution; tnot, RHS computation.
ARTICLE IN PRESS G. Gruosso, M. Repetto / Physica B 403 (2008) 368–371
the first case appealing in the iterative use of magnetostatic solution. The computational cost of all solution phases is reported in Fig. 3. In the following activity the effect of using multilevel FMM on the computational cost of the procedure will be investigated. References [1] [2] [3] [4]
E. Tonti, IEEE Trans. Magn. 38 (2002) 333. G. Giuffrida, C. Gruosso, M. Repetto, IEEE Trans. Magn. 42 (2006) 1503. C. Giuffrida, C. Ragusa, M. Repetto, Physica B 372 (1–2) (2006) 299. J.D. Jackson, Classical Electrodynamics, Wiley, New York, 1998.
371
[5] L. Greengard, V. Rokhlin, J. Comput. Phys. 73 (1987) 325. [6] K. Nabors, J. White, IEEE Trans. Comput.-Aided Des. 10 (1991) 1447. [7] R. Beatson, L. Greengard, A short course on fast multipole methods, in: Wavelets, Multilevel Methods and Elliptic PDEs, Oxford Science Publications, Oxford, 1997. [8] A. Buchau, C. Huber, W. Rieger, W. Rucker, IEEE Trans. Magn. 36 (4) (2000) 680. [9] N.A. Gumerov, R. Duraiswami, E.A. Borovikov, Data structures, optimal choice of parameters, and complexity results for generalized multilevel fast multipole methods in d dimensions, Technical Report UMIACS-TR-2003-28, Institute for Advanced Computer Studies, University of Maryland, 2003.
Physica B 403 (2008) 368–371 www.elsevier.com/locate/physb
Magnetostatic solution by hybrid technique and fast multipole method G. Gruossoa, M. Repettob, a
Politecnico di Milano, Dipartimento di Elettronica e Informazione, I-20133 Milano, Italy Politecnico di Torino, Dipartimento di Ingegneria Elettrica, C.so Duca Abruzzi 24, I-10129 Torino, Italy
b
Abstract The use of fast multipole method (FMM) in the solution of a magnetostatic problem is presented. The magnetostatic solution strategy is based on finite formulation of electromagnetic field coupled with an integral formulation for the definition of boundary conditions on the external surface of the unstructured mesh. Due to the hypothesis of micromagnetic problem, the resulting matrix structure is sparse and integral terms are only on the RHS. Magnetic surface charge is used as source of these integral terms and is localized on the faces between tetrahedra. The computation of the integral terms can be performed by analytical formulas for the near field contributes and by FMM for far field ones. r 2007 Elsevier B.V. All rights reserved. Keywords: Micromagnetics; Fast multipole method; Cell method; Hybrid techniques
1. Hybrid magnetostatic solution The use of an hybrid form of magnetostatic computation based on finite formulation of electromagnetic fields (FFEF), more often called cell method, and on a Green function applied to magnetization sources, has been presented [1,2]. In this formulation an integral boundary condition on the exterior surface of the mesh is related to the magnetization sources. Integral boundary conditions are used to terminate Ampere’s law for all boundary edges expressing external magnetic field circulation as a potential difference evaluated by the Green function integrals at the center of each boundary face. Under the micromagnetic hypothesis, the magnetized region has permeability m0 and, for each time step of the dynamic analysis, magnetization distribution is supposed to be known. The absence of magnetic polarization, that is of a magnetization-dependent field like M ¼ wH, where w is the magnetic susceptibility, eliminates magnetization contribution depending on configuration variables like magnetic vector potential. Thus all integral terms are known and can be carried on to the RHS. Another consequence of this fact is that there are no integral relations between matrix rows Corresponding author. Tel.: +39 011 5647140.
E-mail address: [email protected] (M. Repetto). 0921-4526/$ - see front matter r 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.physb.2007.08.051
and thus matrix structure remains sparse as in the usual formulation of the cell method. This numerical approach has been already applied to demagnetizing field computation on a structured hexahedral grid [3] and is now applied to a tetrahedral mesh. Magnetization is considered to be uniform in each tetrahedron so that magnetic sources are located only on triangular faces. The computation of magnetizationdependent integrals is performed by fully analytical formulas, nevertheless this phase is one of the most timeconsuming parts of each magnetostatic solution. From these considerations stem the idea of using fast multipole method (FMM) to speed up the integral boundary conditions computation.
2. Fast multipole scheme for RHS computation 2.1. Multipole expansion Magnetic scalar potential is computed by an integral operator applied to magnetization-dependent sources. The scalar potential c is expressed by Z cðrÞ ¼ rM Gðr; r0 Þ dOM , (1) O
ARTICLE IN PRESS G. Gruosso, M. Repetto / Physica B 403 (2008) 368–371
where G ¼ 1=ð4pjr r0 jÞ is the usual Green’s function of three-dimensional static problems and rM is the magnetic charge density obtained as rM ¼ rMrSM, and OM is the volume domain containing the magnetization. This term can be evaluated by considering two contributions to the potential: one due to near field and one due to the far field c ¼ cnear+cfar. The first term [2] can be calculated by means of Eq. (1), while the last one can be obtained by resorting to an expansion in multipoles [4]: cfar ðr; W; jÞ ¼
l X l X Mm l Y m ðW; jÞ, rnþ1 l l¼0 m¼l
369
cube well separated from them. The cost of this algorithm is approximately Nf log Ns, where Nf is the number of evaluation points and Ns the number of the sources. Generally this approach is not enough to improve the computational efficiency and multilevel algorithms are introduced [7]; in the present formulation, the number of field points, located only on the surface of the active volumes, is much lesser than the number of sources located in the whole active volumes, therefore only the O(N log N) algorithm has been considered.
(2) 2.3. Data structures
with
where ri, Wi, and ji are spherical coordinates of ith charge qi ¼ siAi, and r, W, and j are the coordinates of field point. Ym l are the normalized spherical harmonics of degree l and order m given by the following: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2l þ 1 ðl mÞ! m m P ðcos WÞejmj , Y l ðy; jÞ ¼ (4) 4p ðl þ mÞ! l
In this section a data structure based on octree [8] subdivision of the domain combined with a bit interleaving techniques [9] will be described. This technique allows an efficient indexing of the cube in the hierarchical tree and then a robust and fast method to find relationship between cubes at each level of the octree. Briefly some characteristic aspects of the method will be outlined: considering a box at level l the whole number of children of this box is 8, numbered from 0 to 7, therefore the index of one of these boxes is unique and it is given by the following relationship:
with the condition
IDXðn; lÞ ¼ ðN 1 ; N 2 ; . . . ; N j Þ,
Mm l ¼
N X
qi rli Y m l ðWi ; ji Þ,
(3)
i¼1
m mn Y m l ðy; jÞ ¼ ð1Þ Y l ðy; jÞ,
(5)
where Pm l ðcos WÞ are the Legendre polynomials of degree l and order m [4]. 2.2. O(N log(N)) fast multipole scheme Several fast multipole schemes have been developed in order to increase the computational efficiency of the method [5–8]. All these methods are based on the same basic principle: in order to make systematic use of multipole expansions, a hierarchy of boxes which refines the computational domain into smaller and smaller regions will be defined. At refinement level 0, the whole computational domain will be considered. Refinement level k+1 is obtained recursively from level k by subdivision of each box into eight equal parts. This yields a natural tree structure, where the eight boxes at the k+1 level, obtained by subdivision of a box at level k, are considered its children. By starting from the generated tree it is possible, for all sources which are ‘‘sufficiently far’’ from a given field point, using the multipole expansion, to summarize their total interaction. According to error estimation, a convenient definition of ‘‘sufficiently far’’ (or ‘‘well separated’’) is a relationship between the diagonal of each cube and the distance ofp field ffiffiffiffiffiffi point from the center of the considered cube: r0 o2 ð3Þw, where r0 is the distance of field point from the local origin and w the length of the edge of the cube. The first algorithm consists in evaluating the far field at all field points by considering the contribution of each
(6)
where l is the considered level, Nj the index of the box at level j given by Nj ¼ 0, y, 7 with j ¼ 1, y, l and n is the global number of the cubes at level l with n ¼ N1(8)(l1)+ N2(8)(l2)+?+Nl1(8)+Nl. For example, by referring to the construction shown in Fig. 1,the cube N1 ¼ 0 and N2 ¼ 5 at level 2 is the cube n ¼ 5 of the level 2 and so on; in this way the uniqueness of the indexing is guaranteed. In this way it is possible to obtain the father box (PRT) for each pair (n, l)by means of PRTðn; l 1Þ ¼ ðN 1 ; N 2 ; . . . ; N j0 Þ,
(7)
with j0 ¼ 1, y, l1. From the previous statement it is possible to compute all the children (CHD) at level l+1: CHDðn; l þ 1Þ ¼ ðN 1 ; N 2 ; . . . ; N jþ1 Þ,
(8)
where Nj+1 ¼ 0, y, 7. In other words the use of octree makes obtaining parent and children indices very convenient. Indeed the above operations are nothing but shift operations in the bit representation of n. Performing a right
Fig. 1. Octree subdivision of the domain and cube hierarchy.
ARTICLE IN PRESS G. Gruosso, M. Repetto / Physica B 403 (2008) 368–371
370
bit-shift operation on n by 3 bits, one can obtain the index of the parent. One can list all indices of the children boxes of n by a left bit-shift operation on n by 3 bits and adding all possible combinations of 3 bits. In order to decode this structure in bit representation, it is necessary to introduce a spatial ordering. Firstly it is possible to define a proper scaling of the point Pj ¼ (xj, yj, zj) as ¯ j ¼ Pj Pmin , P D
3. Results The results, obtained on a cylindrical domain with an impressed uniform magnetization, are outlined only in terms of cpu-time comparison of one magnetostatic solution obtained with three different schemes:
(9)
where Pmin ¼ (xmin, ymin, zmin) is the lower corner of the box enclosing all the domain and D is the length of biggest edge of the box. In this way the system will be scaled into a box with the lower corner in the origin and with edges equal to 1. ¯ j ¼ ðx¯ j ; y¯ j ; z¯ j Þ can The scaled coordinate of the jth point P be represented in binary form as
LU factorization of the matrix and computation of RHS by exact formulas; GMRES with incomplete LU preconditioning and computation of RHS by exact formulas; same GMRES with FMM computation of RHS.
As can be seen from Fig. 2, the combined use of GMRES and FMM sensibly reduce the computational time of one solution. Despite the O(n3) computational burden of LU factorization the O(n) cost of back-substitution can make
x¯ j ¼ ð0:a1 a2 a3 . . . Þ2 , y¯ j ¼ ð0:b1 b2 b3 . . . Þ2 , z¯ j ¼ ð0:c1 c2 c3 . . . Þ2 , where ai ¼ 0, 1, bi ¼ 0, 1, and ci ¼ 0, 1. Instead of having 3 indices for each point it is possible to form a single binary index that represents the same point by means of ordered interleaving of the digits (bit interleaving), so the index becomes P ¼ ð0:a1 b1 c1 a2 b2 c2 a3 b3 c3 . . . Þ2 ¼ ð0:N 1 N 2 . . . N j Þ2 ,
(10)
where Nj ¼ (aj bj cj)2. By starting from these considerations, it is easy to perform the following operations: 1. Finding the IDX containing a given point. For a generic scaled point obtained from Eq. (10) at level l, by referring to a hierarchical structure ordered as in Fig. 1, the following relationship is true: P¯ ¼2 BoxððN 1 N 2 . . . N l Þ8 Þ where the value of the Box function is given by Eq. (6). 2. Finding the center of a given point. The relation between the coordinate of a point and its box index enables easy finding of the coordinates of the center for a given box index n at level l, simply converting this information in binary form, using Eq. (6) and then decomposing the result into the three coordinates:
Fig. 2. Cpu time comparison by using different algorithms for the calculation of the RHS.
n ¼ ða1 b1 c1 a2 b2 c2 . . . al bl cl Þ2 , x ¼ ð0:a1 a2 . . . al Þ2 , y ¼ ð0:b1 b2 . . . bl Þ2 , z ¼ ð0:c1 c2 . . . cl Þ2 . Combining this algorithm with the multipole calculation it is possible to reduce the computational effort due to the proposed formulation.
Fig. 3. Cpu time cost of all solution phases: int, computation of integral coeffcients; mat, matrix computation and assembling; sol, linear system solution; tnot, RHS computation.
ARTICLE IN PRESS G. Gruosso, M. Repetto / Physica B 403 (2008) 368–371
the first case appealing in the iterative use of magnetostatic solution. The computational cost of all solution phases is reported in Fig. 3. In the following activity the effect of using multilevel FMM on the computational cost of the procedure will be investigated. References [1] [2] [3] [4]
E. Tonti, IEEE Trans. Magn. 38 (2002) 333. G. Giuffrida, C. Gruosso, M. Repetto, IEEE Trans. Magn. 42 (2006) 1503. C. Giuffrida, C. Ragusa, M. Repetto, Physica B 372 (1–2) (2006) 299. J.D. Jackson, Classical Electrodynamics, Wiley, New York, 1998.
371
[5] L. Greengard, V. Rokhlin, J. Comput. Phys. 73 (1987) 325. [6] K. Nabors, J. White, IEEE Trans. Comput.-Aided Des. 10 (1991) 1447. [7] R. Beatson, L. Greengard, A short course on fast multipole methods, in: Wavelets, Multilevel Methods and Elliptic PDEs, Oxford Science Publications, Oxford, 1997. [8] A. Buchau, C. Huber, W. Rieger, W. Rucker, IEEE Trans. Magn. 36 (4) (2000) 680. [9] N.A. Gumerov, R. Duraiswami, E.A. Borovikov, Data structures, optimal choice of parameters, and complexity results for generalized multilevel fast multipole methods in d dimensions, Technical Report UMIACS-TR-2003-28, Institute for Advanced Computer Studies, University of Maryland, 2003.