35, NO. 3, MAY 1999. Some realizations of a discrete Hodge operator: A reinterpretation of ... Due to the apparent relationship between (4) and (5), we shall call ...
IEEE TRANSACTIONS ON MAGNETICS, VOL. 35, NO. 3, MAY 1999
1494
Some realizations of a discrete Hodge operator: A reinterpretation of finite element techniques Tim0 Tarhasaari, Lauri Kettunen Tampere University of Technology, Laboratory of Electromagnetics, P.O.Box 692, FIN-33101 Tampere, Finland
Alain Bossavit Electricit6 de fiance, 1 AV. du Gal de Gaulle, 92141 Clamart, France Abstract-In this paper some structures which underlie the numerical treatment of second-order boundary value problems are studied using magnetostatics as an example. We show that the construction of a discrete Hodge is a central problem. In this light we interpret finite element techniques as a realization of the discrete Hodge operator in the Whitney complex. This enables us to view the Galerkin method as a way to set up circuit equations, the metric of space being encoded in the values of branch impedances. I n d e x terms-Finite element methods, differential forms, Hodge operator
I. INTRODUCTION In electromagnetics, as a rule, one uses two different mathematical entities - for instance, the magnetic flux density B and the magnetic field strength H - to describe one and the same physical field. The two entities are needed to express dual properties of the same phenomenon. For simplicitly, without losing generality, let us consider magnetostatics as an example. In magnetostatics B and H fulfill Gauss and AmpBre’s law divB = 0 , curlH = J .
(1) (2)
These are the two basic properties of the magnetic phenomenon. And now, the constitutive law
B=pH
(3)
express how the two “sides” of the magnetic field are related. In terms of differential forms [l]instead of B and H one has magnetic flux b (a 2-form) and magnetic field h (a 1-form). The constitutive law can be written by b=*h,
(4)
where the definition of the *-operator is contrived in order t o embed the parameter 1.1. This reveals now the idea of duality: if one entity, b or h, describing the magnetic field is known, then the *-operator yields immediately the other, dual one. On the discrete level b and h correspond with arrays of degrees of freedom (DoF-arrays) & and ?L, which contain Manuscript received June 1, 1998. T. Tarhasaari, e-mail timo.tarhasaariOcc.tut.fi
the fluxes across some surfaces and the circulations along some paths, respectively. We shall first examine how the Gauss apd Ampkre’s law are fulfilled with the DoF-arrays 8 and h. After this we construct a mesh-dependent matrix * which relates the dual arrays 6 and h such that (5) Due to the apparent relationship between (4) and ( 5 ) , we shall call matrix * the discrete counterpart of the Hodge operator *. Our aim in this paper is to show that the construction of a discrete approximation of the *-operator is a central problem. Starting from the basic definitions we i n t r e duce a discrete Hodge and then show, that the Galerkin method imposes another one although both approaches result in the same solution for a second order boundary value problem. Furthermore, we show how the +operator can be exploited to separate the metric-dependent and metric-independent parts of boundary value problems. 11.
DUALENTITIES
IN DUAL GRIDS
In a n-dimensional space the Hodge (star) operator yields an isomorphism - a one-to-one correspondence between pforms and (n - p)-forms. So, the *-operator relates locally the two sides of the magnetic field t o each other. This is what we formally express with (4). On the discrete level, as is well known, one can’t r e p resent both the 2-form b and the 1-form h simultaneously in the same finite element mesh. However, the problem can be circumvented - this is a key point - if b and h are assigned on a pair of dual grids. (For a same kind of idea, see [21, [31 and [4.) A grid or a primal mesh is a cell complex, either rectangular boxes, simplices, polyhedra etc., and the dual grid or dual mesh to the primal one is another grid with onet e o n e correspondence between pcells and (n - p)-cells. For example, if the primal grid is a simplicial mesh and dimension n = 3, the dual counterpart of a node so, edge sl, facet s2 and tetrahedron s3 is a polyhedron S 3 , a poly onal face 2,a broken segment of line S1 and a vertex s’ , respectively. The dual counterpart of a simplicial mesh is called a burycentric or regular subdivision. The DelaunayVoronoi meshes (exploited e.g. in references [2] and [3]) are another example of this kind of duality. A 2D-example is shown in Fig. 1.l
8-
‘Be aware, the problem of Delaunay meshes with obtuse angles can be ruled out by requiring that the pcell support and its dual ( n - p)-cell support should effectively intersect. 0018-9464/99$10.00 0 1999 IEEE
1495
As the integration of differential forms does not need metric structures, there is no metric involved in (6) or (7): the de R h m maps are independent of met+ In consequence of this in the SI-system the units of b and h are in webers and in amphres, r-p~timly, instead of Wb/m2 and A/m. (The unit of length does not pop up.) IV. EXTERIOR OPERATOR On the left, a simplicial mesh and its dual complex. The dual counterpart of edge s t (p = 1) is the broken segment of line St ( n - p = 1). On the right, Delaunay-Voronoi meshes. Fig. 1.
On the discrete level the field problems, such as our magnetostatic example, can be interpreted as looking for a recipe or rule, how the two DoF-arrays representing the p and (n-p)-forms, one living on the primal and the other on the dual side, are related. As either b or h can be put on the primal side, we should always have two families of formulations, i.e the b and h-oriented methods. However, we shall deal only with &formulations in order to keep the exposition short and compact.
Let us now assume a simplicial mesh and its barycentric subdivision. The sets of simplices are denoted by SP and where dimenthe sets of their dual counterparts by ,!?--P, sion n = 3 and p = 0,. . . ,3. The space of p-forms is denoted by Fp. Obviously we need also spaces for the DoFarrays. So, let us call CP and C n - p the spaces spanned by the DoF-arrays on the primal and dual side, respectively. We shall call de Rham map, and denote by C, the machinery to map F P to C P on all levels p = 0 , . . . , 3 . If w E FP then the scalar which is the ith entry of array Cw is defined by
.
dCw = Cdw V w E F P .
(10)
As the Gauss's law for magnetics (1) can be expressed in the language of differential forms by db = 0, applying the de Rham map yields d C b = Cdb = 0 .
111. DE RHAMMAPS
(cw)i = J w
To impose the Gauss and Ampike's law we need differential operators. So, let us denote the discrete counterparts of grad, curl and div by matrices G , C and D, respectively, and call them with a generic n-9. This way matrix d is the discrete counterpart d e primal side of what is known as the exterior operator d in differential geometry. Between d and C holds
(6)
:s
where simplex s: E SP. (Such integrals make sense without specifying a length, area, etc., over the simplex. This is the point of working with differential forms.) Correspondingly, the dual de Rham map c" maps from FP to C P . If w belongs to Fp and d E @'is the dual counterpart of simplex sn-P E Sn-P,then operator c" is defined by
(11) Thence, in the discrete case the Gauss's law can be given bY db = Di = 0. (12) In other words, (12) states that the sum of the DoF's (i.e. the total flux) over a set of faces which is the boundary of an assembly of elements - i.e. over a so called bounding cycle [5] - is null. Thus, (12) expresses a property of b in a combinatorial manner. In other words the matrix d involves only -1, 1 or O's, which is in consequence of also d being metric independent. At the discrete level the relation between d and its adjoint corresponds to the relation between the matrix d and its transpose d T . Therefore between d and c" should hold dTCw = c"dw V W EF p . (13) And now, this enables us to impose also AmpBre's law in discrete spaces. On the continuous level it is given by dh = j. Thereby, according to (13) we have dTfh = f d h = f j ,
(14)
and now, as k = c" h, we end up with
dTK = CTh =
3.
(15)
Note that array 3 is associated with the dual mesh. A graphical illustration of the relationship between the d operator and the de Rham maps is shown in Fig. 2. Summing up, the discrete magnetostatic problem can be formally expressed as Find vectors b and k such that (5), (12) and (15) hold b = Cb. (8) simultaneously. Notice that (12) and (15) can be imposed ezactly as In other words, operator C yields the array E C 2repre- they involve no metric. (Metric-independent data can be senting b in the primal mesh. _In the same fashion when expressed with integers, and computers do execute integer h E F 1 is inserted, operator C gives the dual side array arithmetics exactly.) This fact emphasizes the central role iiE C l , of the discrete Hodge operator although we have not yet K =eh. (9) defined it explicitly. The de Rham maps are mathematical tools needed to construct a discrete representation of any kind of field, be it scalar, l-form, 2-form or 3-form. As an example consider the field b E F 2 . In this case we have
I
1496
Id
Id
The primal de Rham map maps spaces of pforms onto spaces spanned by DoF-arrays associated with the primal mesh. The dual de Rham does the same but now onto spaces of DoF-arrays associated with the dual mesh. The exterior operator and its discrete counterparts are mappings from level Fig. 2.
The property given in (17) raises a question whether a dual Whitney map corresponding with c' could be found. This is an open question and since we do not need the dual side Whitney elements in this paper, we shall bypass the problem. VI. CONSTRUCTION OF A
c
F P
*
c p
4
4
*I
HODGE
What we have now is enough t o suggest a construction of the +matrix. If (16) is imposed in a weaker sense for only those forms which are in the range of W, i.e.
*CW3 =
ptop+l.
As the *-operator maps pforms to (n-p)-forms, consequently the matrix * should map from C P onto C n - P , Fig. 3. Since there is a one-to-one correspondence between the simplicial mesJl and the barycentric subdivision, the DoFarrays 6 and h are of same dimension. This implies that the matrix * should be a non-singular square matrix.
DISCRETE
t*w3 v3 E c p ,
(18)
then, thanks to (17), that would imply G =
C*W3 Q 3E C p .
(19)
So, if we define
H = &W, (20) this square matrix is a fair approximation of the Hodge operator. There is one such matrix H at each level p = 0,. . . ,3, and using the one for p = 2, we may rewrite ( 5 ) now by
I*
Hz6 =
L.
(21)
According to (7) and (20) the entry i , j of matrix H2 is given by
H2(i,j) =
The relationship between the spaces of forms and the spaces of DoF-arrays. Fig. 3.
In the ideal case we should have *CW = &w
E Fp
.
(16) Informally, this means that if we skirt say from b and then find its discrete representation b, then *b should be precisely the same array as the discrete counterpart of h = *b. This property will guide us in setting up a definition for the discrete Hodge operator. QW
v. WHITNEY MAP To work towards the discrete Hodge operator we need now an operator which maps the other way compared t o C. For this purpose we select the Whitney forms (called also Whitney elements) [6] as the tool to construct pforms from the DoF-arrays in CP. As the background of Whitney elements is reported elsewhere - see for instance references [6] and [7] - we bypass the details here and just call the Whitney map, denoted by W , the metric independent machinery which maps arrays of DoF's onto forms on all levels p E 0,. .. ,3. As an example, let us say we know an array b c_oresponding with fluxes across the facets. In this case W b yields a field b E F2 by using facet elements (P = 2). Between the Whitney map and the de Rham map holds
cw3 = 3 ' d 3 € C P , (17) where p = 0 , . . . , 3 . However, the converse is not true, WC # I, which means there is a class of fields which have the same discrete representation.
/*Wj
(22)
9
i;
where wj is the Whitney 2-form, i.e. the facet element associated with the j t h facet s;. The right-hand side of (22) is the integral of the 2-form on the dual counterpart of facet i. So, in terms of classical vector analysis (22) is given by Hz(i,j) = / - t1. W j , (23) P at where Wj is the facet element associated with the j t h facet, and t is the tangent vector of the broken segment of line named 5:. (The dual counterpart of a primal facet is a broken segment of line in the dual mesh.) A graphical interpretation of the matrices Hp is shown in Fig. 4. F P
.
C
*
c p
*I (e:
p - P
c'
Cn-P
Fig. 4. Graphical illustration of matrix H,which is an approximation of the +operator.
Matrix H does not have precisely the same properties as the *-operator. Ideally, the discrete Hodge operator should retain the property of a local mapping. Matrix H can be inverted, and it has some locality, but it is not a
1497
diagonal matrix, and therefore it is not strictly speaking a local operator. This can be easily confirmed with Fig. 4. Start with a continous field w in the bottom left corn_er of Fig. 4. Next revert the order and construct WH-lCw t o get t o the upper left-corner. This is clearly not the same as *w, since WH-IC w is in the range of W .
A . Example We may now rephrase our discrete magnetostatic problem and state: Find 8 E C2 and ?L E C1 such that
DE = 0 , CTi = j ,
(24) (25)
H26 = 6, (26) hold simultaneously. To impose (24), let us choose b = da, where a is _the potential.2 On the discrete level the same becomes b = Czi, and equivalently to (24), (25) and (26) we may state: find zi E C1 such that CTH2Csi =
holds. The interpretation of this is: seek for an array ii such that: (i) the fluxes,i.e. the sums of the entries of 6 = Csi vanish across all possible assemblies of facets forming closed surfaces (of course, having b = Csi automatically satisfies this), and (ii) on the dual mesh the circulation of H2b fulfills Amphe’s law along all “loops” made of assemblies of the dual counterparts of the primal facets. But the problem has now the same structure as circuit equations. The first property corresponds with Kirchhoff’s law of closed current loops and the second with Kirchhoff’s law of zero voltage around 1 0 0 ~ s . In ~ this light the matrix H2 - including all metric data - is like the impedance matrix Z in circuit theory, with the exception that the metric of space is encoded into H2but not into Z .4 VII. THEGALERKIN METHOD The Galerkin approach for the discrete magnetostatic problem is given by: find si E C1such that CTM2Csi = j, (28) holds. Matrix M2 is a so called mass matrix [7] whose entries are given by M2(irj) = / i W i . W j ,
(29)
R 2From the mathematical point of view this is not necessarily always possible, but the mwmtic flux density b happens t o have enough properties which justifies this move. 31n the case of Amphre’s law, of course, the circulation of magnetic field does not vanish around cycles. However, the idea of Kirchhoff’s law is t o imply a certain property for voltages around loops, The fact that with circuits this voltage happens to be null is very much a secondary point here. 41t is possible t o find a similar kind of interpretation for all kinds of numerical approaches in electromagnetics. For example, see references [2], [8] and [9].
where the W’s axe the Whitney facet elements. Although the matrices H2 and M2 are not the same, the Galerkin method implies the same system matrix as in (27), i.e. CTM2C = CTHzC. (30) This means that the circuit equation interpretation ap._plies to the Galerkin method as well - see also reference ‘ [lo]- and what is more important, similarly as H also the mass matrix M is a discrete representation of the Hodge operator. Thus, the Galerkin method can be understood as a realization of a discrete Hodge operator. The reinterpretation of the Galerkin method makes immediately some strong suggestions of how the finite dement techniques should be linked with circuit eq&&ns. However, more important is that the difference be&een various numerical approaches boils down t o the chosen *-matrix simulating the Hodge operator on the discrete level. (See reference [ll].)In practice, this means that one should be able to write effectively a single software system which is able to find several kinds of numerical solutions for the same boundary value problem just by providing access to several Hodge operator subroutines. VIII. CONCLUSION In this paper it is shown how metric-dependent and metric-independent parts of second order boundary value problems can be separated. As the metric-independent equations are imposed exactly in the sense of circuit equations, the discrete solution is characterized by the metricdependent data which is all involved in the discrete Hodge operator. Therefore the construction of the discrete *operator is a central problem. This reinterpretation of finite element techniques suggests several practical aspects of how finite element software systems should be built. REFERENCES 111 H.Flanders, Diflerential Forms with Applications to the Physical Sciences. Dover Publications, 1989. [2] T. Weiland, “Time domain electromagnetic field computation with finite difference methods,” Int. J. Num. Modell., vol. 9, pp. 295-319, 1996. [3] E.Tonti, “Algebraic topology and computational electromagnetism,” in Proc. of 4th Int. Workshop on Elec. and M a p . Fields, (Marseille, France), pp. 285-294, 1998. [4] C. Mattiussi, “An analysis of finite volume, finite element, and finite difference methods using some concepts from algebraic topology,” J. of Comp. Physics, vol. 133,. pp. 289-309, 1997. [5] 3. Stillwell, Classical Topology and Combanatorial Group Theory. New York: Springer-Verlag, 1980. [6] A. Bossavit, “Whitney forms: a class of finite elements for three-dimensional computations in electromagnetism,” IEE Proc., vol. 135, Pt. A, no. 8, pp. 493-500, 1988. [7] A. Bossavit, Computational Electromagnetism, Variational Formulations, Edge Elements, Complementarity. Boston: Academic Press, 1998. I81 R. Albanese and G. Rubinacci. “Interrral formulation for 3D eddy-current computation using edggelements,” IEE Pruc., vol. 135, Pt. A, no. 7, pp. 457-462, 1988. [9] A. M. Winslow, “Numerical solution of the quasilinear Poisson equation in a nonuniform triangle mesh,” J. Comp. Phys., vol. 1. no. 2., DD. 149-172. 1967. [io] A. Bossavit, “HOWweak is the weak solution in finite element methods?,” IEEE Pans. Magn., vol. 34, no. 5, pp. 2429-2432, 1998. [ll] A. Bossavit and L. Kettunen, ‘‘Yedike schemes on a tetrahedral mesh with diagonal lumping,” Int. J. Numer. Modelling, 1999. In press.
..
..
