HOMOLOGY AND COHOMOLOGY COMPUTATION IN FINITEELEMENT MODELING
M. PELLIKKA† , S. SUURINIEMI† , L. KETTUNEN† , AND C. GEUZAINE‡
Abstract. A homology and cohomology solver for finite element meshes is represented. It is
an integrated part of the finite element mesh generator Gmsh. We demonstrate the exploitationof the cohomology computation results in a finite element solver, and use an induction heatingproblem as a working example. The homology and cohomology solver makes the use of a vectorscalar potential formulation straightforward. This gives better overall performance than a vector
potential formulation. Cohomology computation also clarifies the lumped parameter coupling of the
problem and enables the user to obtain useful post-processing data as a part of the finite element
solution.
Key words. homology computation, cohomology computation, finite element method, lumped
parameter coupling, electromagnetics
AMS subject classifications. 57R19, 58Z05, 65M60, 78M10

1. Introduction. We present a tool for the homology and cohomology computation of domains tessellated with finite element meshes. The tool is an integrated
part of the finite element mesh generator Gmsh [17]. Homology and cohomology computation can be exploited to exhaustively fix the so-called cohomology class of thesolution of a boundary value problem that is solved with the finite element method.
As a concrete application, we demonstrate how such computations greatly benefit themodeling of an induction heating machine.
In boundary value problems that involve the Hodge-Laplace operator, one often
needs to choose the cohomology class of the solution. Such problems are usual inelectromagnetics, which is why our working example is chosen from that field. The
cohomology classes of the problem are generated by the choice of the boundary conditions and by the homology of the problem domain. Informally, homology is about
the quantity and the quality of holes in an object, whether it has voids or tunnels
or both. Relative homology captures whether the object has holes when one “disregards” a part of the model. In the finite element method, the disregarded part is a
subdomain where the solution is fixed by a boundary condition. Cohomology can becharacterized by saying that it assigns quantities to these holes. In boundary value
problems such assignments fix the cohomology class of the solution. For the technicaldefinitions of homology and cohomology spaces, see appendix A.
In the finite element method, typically only a bounded portion of the device andthe surrounding space is modeled. The modeling domain may contain holes, andboundary conditions are assigned to confine the fields and couple them with external
phenomena outside the domain. Further, the domain can be split into many coupled
regions where different approximations and potential formulations are being employed.
These modeling aspects give rise to homology and cohomology and their relative forms
in numerical models, since one is required to assign source quantities to entities thatare absent from the model.
For electrical engineers, an evident manifestation of homology and cohomology
are Maxwell’s equations in their integral form. Classically, electromagnetic phenom2 Tampere University of Technology, Dept. of Electrical Engineering, Electromagnetics, Finland
3 University of Liege, Dept. of Electrical Engineering and Computer Science, Montefiore Institute,

Belgium1

2

M. Pellikka, S. Suuriniemi, L. Kettunen, and C. Geuzaine

ena are assumed to take place in the entire 3-dimensional Euclidean space that has
no holes. However, once boundaries or holes are introduced to the modeling domain,
the coupling with the external electromagnetic phenomena can often be formalizedwith cohomology. For a more detailed discussion, see [18]. To employ so-called magnetic field conforming method [4] in our working example, cohomology computation isneeded to make the use of an efficient vector-scalar potential formulation straightforward. Making it easy to use such a formulation instead of having to resort to a pure
vector potential formulation is the main benefit of the cohomology computation. As
other benefits, it helps to couple the problem with circuit models, provides efficientand accurate post-processing possibilities, and in general gives a structural insight intosuch problems. Methods that exploit homology and cohomology in electromagnetics
are also presented in [11, 10, 32]. As heat conduction is in many respects analogous toelectrostatics, the presented homology and cohomology solver can be exploited in heat
conduction boundary value problems as well. The role of homology and cohomology
have also been recognized in the field of linear elasticity [1, 37]. Thus, homology andcohomology computations might also have applications in that field.
Another related application for cohomology computations is the parametrizationof surface meshes. Good quality surface meshes of objects are desirable in computer
graphics and are a prerequisite for 3D-meshes of good computational quality. To
produce a good surface mesh, a surface parametrization is often needed [17, 30]. A
closed surface cannot be parametrized with a single coordinate chart, for a discontinuity is introduced at the coordinate patch boundary. To avoid this, cohomology
computation can be exploited to find a global conformal parametrization for surfaces
with an arbitrary genus [19]. For compact non-closed surfaces with boundary andholes, relative cohomology might be exploited in a similar manner.
The presented homology and cohomology solver is an integrated part of the finite
element mesh generator Gmsh [17]. It can also be accessed through the application
programming interface (API) of Gmsh. Thus, homology and cohomology computation
can be embedded as a part of the mesh generation workflow. The working example issolved using the finite element solver GetDP [12] which is able to fluently exploit theresults of the cohomology solver in its problem definition. However, any solver canbe programmed to exploit the results, as they are saved along with the finite element
mesh.
The design principles of the presented homology and cohomology solver are thefollowing. First, it operates on large, up to 3 dimensional finite element meshes
that are possibly non-simplicial. It is expected that the usual input is such thatthe reduction techniques are be able to tremendously reduce the size of the algebraic
problem. Second, it produces a basis representation of the homology and cohomology
spaces in addition to Betti numbers. Such output is needed in finite element solvers
as we will demonstrate. Examples of other software that perform homology andcohomology computation include CHomP [6], jPlex [33], and GAP homology [22]packages. The design objectives of these popular packages are slightly different fromours, with less emphasis on problem position in finite element modeling.
2. The homology and cohomology solver for up to 3D meshes. A finite element mesh is easily converted into a regular cell complex 1 [20] which is the
1 A mesh that contains only triangles or tetrahedra can be converted into a special case calleda simplicial complex, which is often used in the introductory expositions of homology computation.
The regular cell complex is required to handle more general finite element meshes, for instance,
containing hexahedra, prisms, and pyramids.

Homology and Cohomology Computation in FE Modeling3

usual input for homology or cohomology computation. In the presented homology
and cohomology solver, the extraction of the cell complex from the finite element
mesh is the first stage of the computation. From the cell complex, integer matrix
representations of so-called boundary operators can be obtained. They are basicallythe volumes-to-faces, faces-to-edges, and edges-to-vertices incidence matrices of themesh.
The general algorithm for homology and cohomology group2 computation is theSmith normal form integer matrix decomposition, but it suffers from high worst casecomputational complexity [36]. Thus, before the matrix decomposition is employed,
the problem size is reduced using a variety of methods that usually run in linearithmic
time [23, 34]. In the homology and the cohomology computation of the finite element
meshes, the problem size is reduced by removing “homologically irrelevant” parts ofthe mesh and the mesh elements in the relevant parts are combined. These operations
keep the homology of the mesh invariant [23]. Heuristically, the problem size is usually
reduced by an order of 103 to 105 depending on both the geometry and the topology ofthe mesh. In the represented homology and cohomology solver, reduction algorithms
described in Pellikka et al. [28] are employed. A short description on the essentials of
the solver is given in the section 2.2.
The input for the solver is a mesh that is generated by Gmsh or imported into
Gmsh. The solver only uses the connectivity information of the mesh elements, andignores the coordinates of the mesh vertices. Therefore, the solver can be utilized withany data set that can be realized as an at most 3-dimensional simplicial complex,
once it is converted to a Gmsh-readable mesh format. In section 2.1, we consider
an example where a point cloud in R3 is interpreted as a 2-dimensional simplicial
complex.
Results of the homology and cohomology solver are the bases for the relative homology spaces Hk (M, S) and the relative cohomology spaces H k (M, S) of the domain
M and its subdomain S for each k from 0 to 3. The bases are represented by aset of k-chains and k-cochains whose cosets 3 form the actual bases. With the finite
element mesh, the bases thus have matrix representations Hk = [z1 , z2 , . . . , zβk ] andHk = [z1 , z2 , . . . , zβk ], respectively, where for βk = dim Hk (M, S) = dim H k (M, S)
holds. 4 The sparse integer coefficient vectors zi and zi contain the information
which k-dimensional mesh cells are a part of the corresponding k-(co)chain zi or z i .
The homology and cohomology solver produces such matrices given a finite element mesh of a computational domain M with a designated subdomain S. The basis
representations are non-unique on two levels. First, the representative chains andcochains can be any from the corresponding coset. Second, the basis can be transformed to an another basis, with different basis cosets, by an unimodular transformation [36]. The basis representation issue is further discussed in Pellikka et al. [29].
2.1. Example: surface reconstruction from a point cloud. A techniqueto construct a surface mesh based on a point cloud in R3 is called alpha shapes [15].
They aim to approximate the “shape” of a point cloud, where the level of detail is
2 Inherently, homology and cohomology have the structure of an abelian group, which is also
solved by the presented homology and cohomology solver. However, in the context of the finite
element method, we are primarily interested on the homology and cohomology vector spaces instead.
The free subgroup of an abelian group is interpreted as the vector space of interest.
3 Two k-chains that are k-cycles belong to the same homology coset, if their difference is a kboundary. This equivalence relation generalizes to the relative k-cycles and k-cocycles.
4 If S is empty, β is called the k:th Betti number of M . Then, β is the number of connected0kcomponents of M , β1 is the number of tunnels through M , and β2 is the number of voids in M .

4

M. Pellikka, S. Suuriniemi, L. Kettunen, and C. Geuzaine

expressed with a parameter α. As extrema, if α = ∞ holds, the resulting shape is theconvex hull of the point cloud. If α = 0 holds, the resulting shape is empty.

Fig. 2.1: A point cloud of 6000 points in R3 whose shape is to be determined.
The usual problem is to choose such an α that a shape with desired characteristics
is recovered. The alpha shape is always a subcomplex of the Delaunay triangulation of
the point cloud. Therefore, the homology solver of Gmsh can be applied to computeits topological features.
In this example, we use a priori knowledge of the shape to find a suitable valueof α. Namely, we insist that the shape consists of two tori. As a single torus has
Betti numbers β0 = β2 = 1 and β1 = 2, we can use the homology solver to find such
α that the Betti number criteria for the shape is met. In Fig 2.1 we have depicted
an example point cloud. Fig 2.2 depicts the shape the point cloud determined withdifferent values for α along with the Betti numbers of the shape. We used the program
Hull [7] for the alpha shape computation.
2.2. Homology and cohomology computation process. The homology andcohomology solver in Gmsh relies on reduction of chain complexes [23] before thecomputation of the actual homology or cohomology groups. It is due to the reduction
algorithms that the computed bases have somewhat small support. See Fig. 7.3 forexample. While the details of the employed reduction algorithms are presented in
Pellikka et al. [28], here we describe the overall strategy how the reduction algorithms
are combined together to make up the solver in Gmsh.
The first step of homology and cohomology computation is to construct a chain
complex [26, 20] from the finite element mesh. The chain groups Ck are formal sumsof the mesh cells of each dimension k. That is, to construct the chain complex froman n-dimensional mesh, first all the n-dimensional mesh elements are enumerated.
Then, the set of their unique n − 1-faces is constructed and the faces are oriented.
Thereafter, the n − 2-faces of the n − 1-cells are processed similarly and so on, until
the 0-dimensional cells are reached. To construct the boundary operators of the chain
complex, the incidence relations of the cells are stored along the process. The chain
complex construction has computational complexity bound O(n log n), where n is thenumber of most abundant k-cells.
Once the chain complex has been constructed, the chain complex reduction takes
place. The reduction algorithms for homology and cohomology computation have sub-

Homology and Cohomology Computation in FE Modeling(a) α = ∞, β0 = 1, β1 = 0, β2 = 1(b) α = 20e3, β0 = 1, β1 = 0, β2 = 1(c) α = 10e3, β0 = 1, β1 = 15, β2 = 1(d) α = 9e3, β0 = 2, β1 = 4, β2 = 2(e) α = 2e3, β0 = 2, β1 = 127, β2 = 1430

(f) α = 0.5e3, β0 = 296, β1 = 99, β2 = 487

5Fig. 2.2: Alpha shapes with different values for α and their Betti numbers.

tle differences, but they all rely on chain equivalences [26]. A chain equivalence mapsa chain complex to a chain complex with equivalent homology and cohomology. The
reduction algorithms apply chain equivalences that reduce the number of chain groupbasis elements computationally efficiently. Effectively, this makes the mesh incidence
matrices smaller in the last stage of the homology and cohomology computation.
The chain equivalences employed in Gmsh have geometric interpretations, thatare depicted in Fig. 2.3. They appear in the following works: [23, 25]. Their computational complexities range from O(n log n) to O(n2 ), where n is the size of a basisof a chain group that is being modified. The algorithms are applied in succession sothat most of the work is done by algorithms with lowest computational complexities.
The reduction techniques are able to reduce the size of the problem by a large

6

M. Pellikka, S. Suuriniemi, L. Kettunen, and C. Geuzaine

factor when the input is a typical finite element mesh. That is, when the Betti
numbers of the chain complex are much smaller than the size of the bases of the chaingroups, i.e. the number of mesh elements. For example, with the mesh in Fig. 2.2f,
which no longer resembles a typical finite element mesh, the reduction factor is only
94.4893 %. Compare that to the reduction factor of 99.9981 % in the computation of
the results in Fig. 7.3. In conclusion, the reduction techniques are a good heuristic
for homology computations within the finite element method.
The reduction algorithms are responsible for the somewhat geometrically smoothappearance of the computed homology and cohomology bases. They remove mesh cells
from the same vicinity instead of in a random manner. This forces the actual homology
and cohomology computation to operate on a chain complex that is geometrically
localized. That is, the employed reduction algorithms are also a good heuristic forthe geometric appeal of the results. Except for local deformation [29] of the homology
basis representatives, we do not apply more advanced methods to attain appealingbases with respect to some criteria. However, there is a lot of research on such methods[8, 38, 9, 16].

Fig. 2.3: Geometric interpretations of two chain equivalences used in the homology
computation. In the upper one, a cell on the boundary of the domain is “pushedin.” In the lower one, two cells are combined. The chain equivalences used in thecohomology computation have similar interpretations, but in the dual mesh.

Once the chain complex has been reduced, its homology or cohomology groups
are computed using the Smith normal form integer matrix decomposition. The computation concentrates on the short exact sequence of abelian groupsijkk
0 → Bk −→Zk −→ Zk /Bk → 0.

(2.1)

The group Bk is the codomain of the k + 1-cells to k-cells incidence matrix of themesh, i.e. the group of k-boundaries. The group Zk is the kernel of the k-cells tok − 1-cells incidence matrix of the mesh, i.e. the group of k-cycles. The quotient group
Zk /Bk is by definition the k:th homology group Hk of the chain complex: the groupof equivalence classes of non-bounding k-cycles.

Homology and Cohomology Computation in FE Modeling

7

The Smith normal form is then used in three steps. First, to find out the kernels
and the codomains of the incidence matrices. Second, to construct the matrix representation of the map ik . Third, the matrix representation of ik is decomposed to solve
the structure of the k:th homology group Hk = Zk /Bk ≃ Zk / im ik , where the lastisomorphism results from the algebra of groups [26]. The procedure for cohomology
computation is analogous.
3. Working example: induction heating with magnetic field oriented
formulation. To introduce and motivate the homology and cohomology computations in the finite element method, we take an eddy current problem that models an
induction heating machine as an example. We use magnetic field oriented formulationand consider both current- and voltage-driven variants.
We start the problem definition with the problem domain. Our primary interest
in modeling sense are the eddy currents in the workpiece, while the model also needs
to take the inductor coil configuration in the surrounding air into account. Due toa large non-conducting air region, a magnetic field oriented formulation that treats
most of the air with a scalar potential only is likely more efficient than a magnetic flux
density oriented one, where a vector potential is employed throughout the domain [13].
The topological aspects of the 3D modeling domain M ⊂ R3 in the Euclidean
space5 are the main challenge of using a h-oriented formulation. The domain is divided
into conducting and non-conducting subdomains Mc and Ma , respectively, depicted
in Fig. 3.1, left.
ζe

Sj

Mc1

ζj

Ma

Σ

Mc2

Se

Fig. 3.1: Topology of the induction heating problem.

Topologically interesting aspects of the domain M are described as follows. The
conducting domain Mc = Mc1 ∪ Mc2 consists of two parts, the inductor Mc1 and theworkpiece Mc2 . The workpiece has a tunnel through it, and the inductor pierces the
non-conducting domain. Therefore, the non-conducting domain Ma has two tunnels
through it. These features are captured by the homology and cohomology of thedomains.
For clear and uniform notation, we use differential forms in this treatment. See
the appendix B for a crash course on differential forms on the Euclidean space R3 . Let
F k (M ) denote the space of differential k-forms on M . Moreover, let h, e ∈ F 1 (M ) andj, b ∈ F 2 (M ) denote the magnetic field, electric field, current density, and magnetic
flux density, respectively. The tensor fields µ and ρ are maps F 1 (M ) → F 1 (M ) called
5 More generally, M is a compact oriented 3-dimensional Riemannian manifold (M, g) with boundary that can be covered with a single chart. That is, the notation and the treatment in this paper
remain valid under the change of the metric tensor g and under the change of charts of M .

8

M. Pellikka, S. Suuriniemi, L. Kettunen, and C. Geuzaine

permeability and resistivity, respectively. The operators d, t, and ⋆ are the exterior
derivative, trace, and Hodge operators, respectively; see the appendix B for theirdefinitions.
Our problem statement on M reads: find the fields h, e ∈ F 1 (M ) and j, b ∈2
F (M ) that satisfy the magnetoquasistatic approximation to the Maxwell’s equations
and material relations
de = −∂t b,
dh = j,db = 0,

(3.1)
(3.2)

b = ⋆µh,e = ρ ⋆ j.

(3.3)

In addition, we require the following local constraints and boundary conditions to
hold:dh = 0 on Ma ,tj = 0 on (∂M ∪ ∂Mc )/Sj ,(3.4)

tb = 0 on ∂M.

te = 0 on Sj ,(3.5)

where Sj = ∂M ∩ ∂Mc and Se = ∂M ∩ ∂Ma hold, that is ∂M = Se ∪ Sj holds. With
such boundary conditions, the problem qualifies as a circuit element [35].
To couple the problem with a driving circuit, we need to pose a non-local constraint to take the source current or source voltage into account. Such constraints fix
the cohomology class of the problem. For the source current Is driven problem werequireZ
h = Is , [ζj ] ∈ H1 (Se )
(3.6)
ζj

It follows from the long exact homology sequence [26] that ζj is a boundary of any
surface Σ in M that isolates the two terminals on Sj . The coset [Σ] of such surface
belongs to the homology space
H2 (M,
R in Fig. 3.1, right, ∂Σ qualifies
R example
R Se ). ForR
as such ζj . Then, Is = ζj h = ∂Σ h = Σ dh = Σ j holds by Stokes’ theorem.
Alternatively, we can drive the problem by a source voltage Vs and requireZe = Vs , [ζ] ∈ H1 (Se , ∂Se )
(3.7)
ζe

That is, ζe is any path along Se between the two inductor terminals on Sj , for example
the one in Fig. 3.1, right.
Remark: By the Lefschetz duality theorem the spaces H1 (Se ) and H1 (Se , ∂Se )are isomorphic [26, 18]. This induces a duality between current driven and voltage
driven problems.
4. Vector-scalar potential formulation of the eddy current problem.
In boundary value problems in 3-dimensional domains, the unknown field is often avector field. The computational burden of such problem is relieved if the unknown
vector field can be expressed as a gradient of a scalar potential. Such a potential is
known to exists, but only locally, if the vector field is curl-free [3]. In the terminology
of differential forms, one seeks for a potential 0-form of a closed 1-form. In the eddy
current problem, the magnetic field 1-form h is closed in the non-conducting domain
Ma . Thus, h can be locally expressed as an exterior derivative of a potential 0-form
in Ma . However, in the conducting domain Mc , one still needs to solve for a 1-form.

Homology and Cohomology Computation in FE Modeling9

Such formulation in eddy current problems is called T − Ω formulation, which we will
now employ.
The T -Ω potential formulation in magnetodynamics draws its name from themixed use of the electric current vector potential T and the magnetic scalar potentialΩ [31]. For a better insight, instead of writing the formulation in terms of T andΩ, we first write it in terms of magnetic field h. This prevents the structure of the
problem from being obscured by the excess of notation. Later, we describe a basis of
the function space from which an approximate solution for h is sought. The basis will
have terms resembling the potentials T and Ω, and terms arising from the homologyof the non-conducting subdomain Ma .
In our example problem, the function space from which we look for the solutionh ish ∈ H(M ) = { h ∈ F 1 (M )|dh = 0 on Ma , tdh = 0 on (∂M ∪ ∂Mc )/Sj ,t(⋆µh) = 0 on ∂M, t(ρ ⋆ dh) = 0 on Sj } (4.1)

together with (3.6) or (3.7) depending on the source quantity. The following weak formulation for the equations (3.1) to find h ∈ H(M ) is obtained by using the equations(3.2) and (3.3) and by taking the b-side constraints (3.5) in (4.1) into account:ZZZt(e ∧ h′ ) = 0 ∀ h′ ∈ H(M ),(4.2)
ρ ⋆ dh ∧ dh′ +∂t ⋆ µh ∧ h′ +
Se
Mc
M
Z
⋆µh ∧ h′ = 0 ∀ h′ ∈ H(M ).
(4.3)
−
M

As dh = 0 holds in Ma , we wished to express h exhaustively in terms of a magnetic
scalar potential φ in Ma , so that h = dφ holds. This would be computationallyefficient. However, we cannot, which we know from de Rham’s first theorem [18]:Z(dh = 0 =⇒ h = dφ) ⇐⇒h = 0 ∀ z ∈ H1 (Ma ).
(4.4)
z

In our case the integral condition fails to hold at two places. First, integrals of h over
1-cycles around the inductor should equal to the net current in the inductor. Second,
integrals of h over 1-cycles that go through the tunnel in Mc2 should equal the net
current in the workpiece around the tunnel, see Fig. 5.1, left. Two cosets of such
1-cycles are a basis for the homology space H1 (Ma ).
With the aid of the homology and cohomology solver, one can automaticallyconstruct the missing part of an approximate subspace of H(M ) in Ma : the edgebased cohomology basis functions that are often called “thick-cuts” [24, 21]. The
approximate subspace of H(M ) is still mainly spanned by the nodal finite element
shape functions in the non-conducting domain Ma , but the two cohomology basis
functions enables one to have non-zero circulations of the magnetic field h aroundthe inductor and the workpiece. Typically, the computed cohomology basis functions
have a small support, and thus do not excessively populate the system matrix.
Let us have a finite element mesh on the domain M . Let N (M ) denote the set ofmesh nodes on M , and let E(Ma ) denote set of mesh edges on Ma and let E(Mc /Ma )denote the set of mesh edges on Mc that are not on Ma . For each node and edge
we associate a Whitney 0-form ni and Whitney 1-form ei , respectively. These form abasis for Whitney spaces W 0 (M ) and W 1 (M ) [4].

10

M. Pellikka, S. Suuriniemi, L. Kettunen, and C. Geuzaine

Then, a field that will suffice to represent the approximate solution is of the form
X
X
X
X
h=
φi dni +ti e i + I 1
z1i ei + I2
z2i ei ∈ W 1 (M ) ∩ H(M ),(4.5)
i∈N (M )i∈E(Mc /Ma )i∈E(Ma )i∈E(Ma )where the coefficients I1 and I2 of the cohomology basis functions correspond to thevalues of the integrals of h over two linearly independent 1-cycle cosets in H1 (Ma ),and W 1 (M ) ∩ H(M ) is the approximate subspace of H(M ). The sparse coefficient
vectors z1 and z2 represent the 1-cochains that represent a computed basis {z 1 , z 2 }of the space H 1 (Ma ). Note that the interpretation of the unknown coefficients I1 andI2 depends on the basis of H 1 (Ma ).
Remark: The h-oriented constraints in (4.1) are satisfied by the choice of basis for the function space H(M ). The b-oriented constraints are satisfied “weakly”
by the absence of boundary integrals over Sj and ∂M in equations (4.2) and (4.3),
respectively [5].
5. Choosing a useful basis for the cohomology space H 1 (Ma ). In theexample problem, we would like that the coefficients I1 and I2 have an evident interpretation for the engineer. For example, in order make it easy to fix the source
current or the source voltage, the coefficient I1 should correspond to the net current
through the inductor Mc1 . Also, if the coefficient I2 would correspond to the net
current through the workpiece, one could read off useful post-processing informationright from the solution, without any actual post-processing! For example, if our apparatus was a transformer, the coefficient I2 would indicate the net short-circuit current
in the secondary winding. In order to make such connotations, the engineer needs to
be involved in the choice of the cohomology basis.
However, in practice it is actually easier to choose the homology basis of therepresentative 1-chains and then adjust the cohomology basis accordingly. Let thematrix H1 = z1 z2 . . . zβ1 , β1 = dim H1 (Ma ), represent a computed basis forH1 (Ma ), where zi are coefficient vectors of the representatives of the basis 1-chains zi .
For any unimodular matrix U, the matrix H1 U also represents a basis for the space
H1 (Ma ). The intuitive basis representation for H1 (Ma ) would be the one where therepresentatives zi only loop around a single tunnel in Ma once, and loop around noother tunnels. Then, one can directly relate the circulation of the magnetic fieldh over the representative zi to the net current Ii in the tunnel. That is, the basis
would correspond to the source current Is and net current Ii in the workpiece. Sucha basis is represented in the Fig. 5.1, left. Unfortunately, the homology solver cannot
distinguish such loops without help from the user: Effectively, the user needs to feed
in a unimodular matrix U to produce such a desired basis from the computed one.
For example, if the homology solver produced the basis representation on the righton the left can be attained by the (2 × 2)-matrix
in Fig.
 5.1, the basis representationU = 1 −1; 0 −1 .
Since H2 (M ) = H1 (M ) = 0 holds in our example problem, an automated method
to obtain such clear basis for the space H1 (Ma ) does exist. We can infer from thelong exact homology sequence [26]
j∗

i

∂

∗
∗
H1 (M ) = 0
H1 (Ma ) −→0 = H2 (M ) −→ H2 (M, Ma ) −→(5.1)

that the spaces H2 (M, Ma ) and H1 (Ma ) are isomorphic6 . Therefore, we obtaina basis representation for the space H1 (Ma ) by applying the boundary operator
6 Furthermore,

H2 (M, Ma )≃H2 (Mc , ∂Mc

∩

∂Ma )

holds

and

moreoverHomology and Cohomology Computation in FE Modeling

11

z1

z1

z2

z2

Fig. 5.1: Two basis representations for the homology space H1 (Ma ).

to the computedof the space H2 (M, Ma ). That is, representative basis 2-chains
H1 = D2 H2 = D2 z1 D2 z2 . . . D2 zβ2 holds, where H2 is the basis representation matrix for the space H2 (M, Ma ) produced by the homology solver, and D2 isthe faces-to-edges incidence matrix of the mesh. Such a basis is guaranteed to beequivalent to the one in Fig. 5.1, left.
Once we have obtained a desired basis representation H̃1 := H1 U for the space
H1 (Ma ), it is possible to havespace basis in the following a compatible cohomology
sense. Let the matrix H1 = z1 z2 . . . zβ1 be a basis representation for the space
H 1 (Ma ) produced by the cohomology solver. Then one can automatically obtain abasis representation H̃1 = H1 V for the space H 1 (Ma ) that satisfies
I = H̃T1 H̃1 = H̃T1 H1 V ⇐⇒ V = (H̃T1 H1 )−1 = ((H1 U)T H1 )−1(5.2)

once the β1×β1 -matrix H̃T1 H1 is known. Then, H̃1 is the cobasis of the basis H̃1 , sinceZ
zTi zj =[z j ] = δij , 1 ≤ i, j ≤ β1
(5.3)
[zi ]

holds for the bases. In Gmsh, one can request a cohomology cobasis for a computed
homology basis at the post-processing stage of the (co)homology computation.
6. Realization of the non-local constraints and the linear system. So
far, we have left the boundary term in equation (4.2) unspecified. It turns out thatit can be related to the non-local constraints of the example problem. Consequently,
we can write down the linear system that is either constrained by the source voltage
or by the source current. The other will be a part of the solution to the system.
Let {z1 , z2 } represent the basis of the space H1 (Ma ) in Fig. 5.1, left, and let
1{z1 , z2 } representof the space
P H (Ma2 ).i Then, the cohomology1 basis
P the cobasis
1 i2
1functions E = i∈E(Ma ) zi e and E = i∈E(Ma ) zi e of the function space W (M )∩H(M ) correspond to the source current Is and the net current Ii in the workpiece,
respectively.
Let us take a closer look to the boundary term in (4.2). On the boundary ∂M ,the field te can be expressed by a scalar potential. As dte = −∂t tb = 0 holds by
H2 (Mc , ∂Mc ∩ ∂Ma ) ≃ H2 (Mc1 , ∂Mc1 ∩ ∂Ma ) ⊕ H2 (Mc2 , ∂Mc2 ∩ ∂Ma ) holds.
The first isomorphism holds by the excision theorem [26], and the second since Mc1 and Mc2 are notconnected.

12

M. Pellikka, S. Suuriniemi, L. Kettunen, and C. Geuzaine

the boundary condition (3.5) and dim H1 (∂M ) = 0 holds, it follows from de Rham’s
first theorem (4.4) that dϕ = te holds for a scalar field ϕ ∈ F 0 (∂M ). Then fore, h′ ∈ F 1 (M ),ZZ′
dϕ ∧ th′(6.1)
t(e ∧ h ) =
Se
Se
Z
Z
=
ϕ ∧ tdh′ −d(ϕ ∧ th′ )
(6.2)
Se

Se

holds. The first term vanishes, since for any h′ ∈ H(M ) it holds that tdh′ = tj = 0 on
Se . The voltage ϕ is a constant on each connected component of Sj , with a voltage
difference V1 = ϕ1 − ϕ2 .
Let the boundary ∂Sj = ∂Se be representedRby a 1-chain
R c1 − c2 , where c1 andc2 are homologous in the subdomain Ma , and c1 h = − c2 h′ = I1 hold for any
h′ ∈ H(M ). That is, c1 and c2 belong to the same homology class of H1 (Ma ),represented by z1 . Because zT1 z1 = 1 and zT1 z2 = 0 hold, we obtain as the boundary
term of the equation (4.2)
(
Z
Z
Z

 Z
−V1 when h′ = E1h′ =h ′ − ϕ2t(e ∧ h′ ) = −ϕ ∧ th′ = − ϕ1.
0when h′ 6= E1
c2
c1
Se
∂Se
(6.3)
Similar developments without the connection to the (co)homology were made in [14,13].
Now, by plugging in the approximation in equation (4.5) to the equations (4.2)
and (4.3) and treating the problem as time harmonic, we obtain the following block
linear system:   
0Adn,dn Ae,dn AE1 ,dn AE2 ,dn
φ0 Adn,e1 ,e2 ,e   t 
A
A
A
e,e
E
E

   (6.4)
Adn,E1 Ae,E1 AE1 ,E1 AE2 ,E1  I1  = V1  ,0I2Adn,E2 Ae,E2 AE1 ,E2 AE2 ,E2

where the blocks of the complex matrix correspond to the different kinds of basis
functions of the space W 1 (M ) ∩ H(M ).
Since z1 is homologous to ζj , which was the boundary of a surface Σj that isolates
the terminals on Sj , the current-driven problem is obtained by fixing I1 equal to the
known source current Is . Then, the equation involving V1 plays no role in the solutionof the problem. However, once the problem has been solved, one can obtain the valueof V1 by substitution. Analogously, since the 0-chain ∂ζe lies on the distinct terminals
on Sj , the voltage-driven problem is obtained by fixing V1 equal to the known source
voltage Vs . Then, the net current I1 will be a part of the solution. Note that inboth current- and voltage-driven problems, the net current I2 through the workpieceis obtained as a part of the solution.
Remark : The linear system (6.4) is actually an affine relation between the current
I1 and the voltage V1 .
Remark : If the workpiece is interpreted as the secondary winding of a transformer,
one could also drive the problem with either net current I2 through the winding or
voltage V2 across the winding in addition to either I1 or V1 . This is done by insistingthat there are two electrodes very close to each other in the workpiece Mc2 , that

Homology and Cohomology Computation in FE Modeling13

creates a void in the domain M . The electrodes would belong to Sj , while the sidesof the void would belong to Se , see Fig. 6.1. Then again since zTi zj = δij holds, theboundary term of the equation (4.2) would be

′
1

Z
Z
−V1 when h = E′′t(e ∧ h ) = −ϕ ∧ th = −V2 when h′ = E2
(6.5)

Se
∂Se

0
when h′ 6= E1 and h′ 6= E2

However, then the current leads from the electrodes to the circuit are excluded fromthe model. Thus, the model assumes that negligible electromotive force is induced
to the current leads and the contribution to the magnetic field by the current in theleads is negligible [35].

Sj

Mc1
Ma
Mc2

Se

Se

Sj

Fig. 6.1: A void in the workpiece Mc2 that is excluded from the actual model. The
electrodes are part of Sj while the insulating boundary of the void is a part of Se .

7. Numerical example. We consider an example problem whose geometry is
depicted in Fig. 7.1. For such a geometry, to imagine a suitable basis representation
for the space H 1 (Ma ) with mind’s eye is a difficult task. Thus, having a provenalgorithm for the task is a considerable advantage. Note that from a topological pointof view, the problem is still equivalent to the Fig. 3.1.
In our numerical example, the workpiece is made of aluminum, and the conductor
is copper. AC current of 1 A at 50 Hz frequency is led to the conductor. The thicknessof the workpiece is a bit larger than the skin depth at that frequency.
This example is available at a website [27] together with Gmsh and GetDP script
files. They demonstrate how to perform the required homology and cohomology
computations in Gmsh and how to input the approximation in the equation (4.5) of
the function space H(M ) and the equations (4.2) and (4.3) to the finite element solver
GetDP.
In Figures 7.2 and 7.3 are computed representations for the basis {z1 , z2 } of thespace H1 (Ma ) and for the cobasis {z 1 , z 2 } of the space H 1 (Ma ). The bases are suchthat z1 and z 1 correspond to the net current in the inductor, and z2 and z 2 correspondto the net current through the workpiece. Therefore, the coefficients I1 and I2 in theequation (4.5) will be equal to them. Note how the coefficients of z 1 equal 6 inside
the workpiece to ensure that the incidence matrix of the bases is an identity matrix.
In Table 7.1 we compare the homology and the cohomology computation timefor various meshes of the geometry. The computation time overhead added by thehomology and cohomology computations are of the order of magnitude of the mesh

14

M. Pellikka, S. Suuriniemi, L. Kettunen, and C. Geuzaine

Fig. 7.1: Induction heating geometry and surface mesh of Mc .

Fig. 7.2: A representation of the basis of the homology space H1 (Ma ).

Fig. 7.3: A representation of the basis of the cohomology space H 1 (Ma ).

15

Homology and Cohomology Computation in FE Modeling

Table 7.1: Comparing computation times with different meshes. Number of tetrahedra in the mesh; Time to generate the mesh; Number of tetrahedra in Ma and Mc
which are the cohomology and homology computation domains, respectively; Time tocompute the cohomology H 1 (Ma ) and the homology H2 (Mc , ∂Mc ∩ ∂Ma ); Degrees of
freedom in the T − Ω-formulation and in the A − V -formulation for the same problem;
Solution time for the problems using a direct solver; Impedance Z = U1 /I1 of theT − Ω solution.
#Tet. M
Mesh time#Tet. MaH 1 time#Tet. McH2 timeT −Ω#DoF’stime

A−V#DoF’s
time

Re(Z)
Im(Z)

Mesh 1106 141
2.2 s

88 8495.9 s

17 2922.0 s

61 568301 s

257 148
3 823 s0.532 mΩ1.01 mΩMesh 2148 962
3.4 s

123 977
8.2 s

24 9853.2 s

88 406666 s

366 252–0.485 mΩ0.962 mΩMesh 3268 238
8.1 s

232 22717.9 s36 0119.2 s

179 514
2 126 s649 332–

0.494 mΩ0.932 mΩMesh 4392 218
13.8 s299 328
26.7 s

84 534
15.3

279 490
4 920 s949 848–0.494 mΩ0.918 mΩgeneration. Note that the cohomology computation domain is not the full mesh,but just the air subdomain Ma , and the homology computation domain is just the
conducting subdomain Mc .
We also compare the number of degrees of freedom between the T −Ω-formulationand the A − V -formulation [13] of the entire problem. In the A − V -formulation of theeddy current problem, one seeks for a vector field in the whole domain M : There’s
one DoF for each edge in the mesh of M , and one DoF for each node in the mesh of
the conducting subdomain Mc . When the A − V -formulation is used, (co)homology
computation typically isn’t needed.7
As a conclusion of the results in Table 7.1, the benefits of being able to use
T − Ω-formulation outweigh the additional cost of the homology and cohomology
computation.
In Fig. 7.4 we demonstrate
the effect of the cohomology basis functions E1 andP
2
E to the scalar part i∈N φi ni of the function space. The values of the scalar part
have a jump equal to I1 z1i and I2 z2i across the edges i ∈ E(Ma ). In Fig. 7.5 is thecomputed current density j = dh of the problem, where h ∈ W 1 (M ) ∩ H(M ).
8. Conclusions. A tool for the homology and the cohomology computation for
finite element meshes and exploitation of its results in computational electromagnetics
was presented. The solver is an integrated part of the finite element mesh generator
Gmsh. As such, the homology and cohomology computation step can be seamlessly
incorporated to the finite element modeling workflow. The performance of the solver
is found to be good when the input is a typical finite element mesh. The output of
7 However, if we interpreted the workpiece as the secondary winding of a transformer, one could
use cohomology basis functions arising from the space H 1 (Mc , Sj ) to drive it with voltage or net
current.

16

M. Pellikka, S. Suuriniemi, L. Kettunen, and C. Geuzaine

Fig. 7.4: Real and imaginary parts of the computed magnetic field h ∈ H(M ).

Fig. 7.5: Real and imaginary parts of the computed current density j = dh, whereh ∈ H(M ).

the solver is associated with the input mesh, and thus can be easily exploited within
a finite element solver.
The cohomology solver can be used to produce cohomology basis functions for thefinite element method. Such basis functions make the use of vector-scalar formulationsof boundary value problems straightforward. The degrees of freedom associated withthe cohomology basis functions provide useful post-processing information, as they
can be related to the lumped parameter model of the problem.
In this paper, the exploitation of homology and cohomology computation wasdemonstrated by a T − Ω-formulation of an eddy current problem, which provides aconsiderable performance increase when compared to the A − V -formulation of the
problem. Without the cohomology computation, the usage of the T − Ω-formulation
is more cumbersome in many cases. Homology and cohomology computation also
clarified the circuit coupling of the problem.

Homology and Cohomology Computation in FE Modeling17

Acknowledgments. This work was supported in part by the Belgian Science
Policy (IAP grant P7/02) and the Walloon region (WIST3 grant “DOMHEX”).
Appendix A. Homology and cohomology spaces. Briefly, the homology
and the cohomology of a finite element mesh is described as follows. Based on a meshon a domain M , vector spaces Ck (M ) of k-chains ck areP
defined to be formal sums of
k-dimensional mesh cells σ k . For example 1-chain c1 = σ1 ∈M ci σi1 , where ci ∈ R, isia formal sum of the edges in the mesh of M . Thus, on a given mesh a k-chain ck canTbe represented by its coefficient vector c = c1 c2 . . . cNk . Boundary operator ∂kmaps k-chains to (k − 1)-chains. They canPbe represented as the incidence matrices
Dk of the mesh. Then for example ∂1 c1 = σ0 ∈M (D1 c)i σi0 is a sum of nodes σi0 thaticonstitute the endpoints of the 1-chain c1 .
A k-chain whose boundary is the zero element of Ck−1 (M ) is called a k-cycle andoften denoted by zk . A k-chain that is a boundary of a k + 1-chain is called a kboundary. That is, k-cycles constitute the space ker(∂k ) and k-boundaries constitutethe space im(∂k+1 ). Two k-cycles are called homologous if their difference is a kboundary. Such relation on k-cycles is an equivalence relation and thus can be used
to construct a quotient space whose elements are cosets: equivalence classes of kcycles.
Cochains are linear maps from chains to real (or complex) numbers. That is, they
assign a scalar quantity to the given chain. Integration of a field over a domain is an
example of a cochain when the field is fixed, but the integration domain varies.
A k-cochain ck can be represented as a complex or real coefficient vector ck suchthat for a given k-chain ck , ck (ck ) = (ck )T ck holds. The coboundary operator δkmaps k-cochains to k + 1-cochains and it is defined by δk ck (ck+1 ) = ck (∂k+1 ck+1 ).
Therefore, their matrix representations are DTk+1 . The elements of the spaces ker(δk )and im(δk−1 ) are called k-cocycles and k-coboundaries, respectively.
Homology and cohomology spaces Hk (M ) and H k (M ) are the quotient spaces
Hk (M ) = ker(∂k )/im(∂k+1 ) = {ck + im(∂k+1 ) | ck ∈ ker(∂k )},kkk

H (M ) = ker(δk )/im(δk−1 ) = {c + im(δk−1 ) | c ∈ ker(δk )}.

(A.1)
(A.2)

For example, the elements of Hk (M ) are represented by k-chains that are nonbounding k-cycles, such as loops around tunnels, or surfaces that enclose voids. One
can add to such k-chain any bounding k-chain, and it still represents the same elementof Hk (M ).
Cosets Ck (M, S) = {c + Ck (S) | c ∈ Ck (M )} are called k-chains relative to a
subdomain S ⊂ M . Usually, a choice of representative is made which zeroes thepart of the chain that lies in S. 8 With such choice, the representations of therelative boundary operator matrices are obtained by removing rows and columns of theboundary operator matrices Dk that correspond to mesh cells lying on the subdomain
S. From the definition of relative k-chains, spaces of relative k-cochains, relativehomology spaces, and relative cohomology spaces can be induced. For example, a
relative non-bounding 1-cycle would draw from S to S so that it wouldn’t enclose asurface, part of whose boundary might lie in S.
Appendix B. Differential forms on M ⊂ R3 . In the Euclidean space R3 ,the algebra and the analysis of the differential forms is very similar to the classical8 Such choice disregards torsion in the relative chain groups C (M, S). In the finite elementkcontext, this is not a loss since boundary conditions are not typically assigned to a multiple of aboundary.

18

M. Pellikka, S. Suuriniemi, L. Kettunen, and C. Geuzaine

vector analysis. The main difference is the different point of view: Differential forms
are functionals on vectors rather than vectors. However, once these functionals are
given a component representation, the operations on the component vectors resemble
the operations on the vector fields in the vector analysis. Due to this similarity, it
may seem that the differential forms are just glorified vector fields. However, evenin the Euclidean 3-space, differential forms have some virtues: While the 0-, 1-, 2-,
and 3-forms have clearly different geometric nature, the operators acting on them
are uniform and have a uniform notation. Therefore, differential forms make a goodcompanion to the concepts of the homology and the cohomology.
A differential 0-form ω 0 ∈ F 0 (M ) on a subset M of the Euclidean space R3 is adifferentiable map ω 0 : M → F to a scalar field F of real or complex numbers.
A tangent vector vp of M ⊂ R3 is a vector v ∈ R3 = Tp M that is “bound” toa point p ∈ M . The vector space Tp M is actually just R3 , but denoted in such a
way to emphasize that it is a vector space
to a point p of M . Such
P3 that is bounditangent vector is often denoted vp =
v
∂/∂x
,
where
{∂/∂x1 , ∂/∂x2 , ∂/∂x3 }
i=1 i
3
1
2
is the standard basis of R with Cartesian coordinates x , x , and x3 . Informally, atangent vector can be considered as a “line element” at a point p of a curve lying in R3 .
A pair of tangent vectors up , vp can be considered to form a “surface element” thathas the coefficients of a vector normal to the surface they span: n = u × v. Similarly,
three tangent vectors up , vp , wp form a “volume element” that is represented by ascalar volume V = u × v · w of their span.
A differential 1-form ω 1 ∈ F 1 (M ) is a differentiable linear map ω 1 : Tp M → F.
That is, it takes a tangent vector as an argument and produces a scalar. A differential
1-form can be defined by providing three scalar functions ω 1i : M → F that are theP3“components” the 1-form: ω 1 = i=1 ω 1i (p)dxi at p, where informally, the notation
dxi denotes that
P3 1-forms are integrands of a line integral. The evaluation is then givenby ω 1 (vp ) = i=1 ω 1i (p)vi = ω 1 (p) · v ∈ F. Here and in the following paragraphs, we
use the symbols · and × merely to express the computation rules of the componentsof the differential forms. That is, they do not denote any binary operators that haverepresentation independent meaning.
A differential 2- and 3-forms are similar maps but they take as arguments twoand three tangent vectors, respectively. Or informally, the arguments are surface andvolume elements, respectively. A 2-form has three components: ω 2 = ω 21 (p)dx2 dx3 +ω 2 (p)dx3 dx1 +ω 23 (p)dx1 dx2 . It is evaluated by feeding in a surface element: ω 2 (up , vp ) =
P23
2
2
3
3
1
2
3
i=1 ω i (p)ni = ω (p) · n. A 3-form has just one component ω = ω (p)dx dx dx .
3
3
It is evaluated for a volume element: ω (up , vp , wp ) = ω 1 (p)V .
The wedge product ∧ of a k-form and an l-form produces a k + l-form whenk + l ≤ 3, otherwise it produces zero. The components of ω 0 ∧ η k are ω 0 (p)η(p) for
any k, for 1-forms the components of ω 1 ∧ η 1 are ω 1 (p) × η 1 (p), for a 1-form and a
2-form the component of ω 1 ∧ η 2 is ω 1 (p) · η 2 (p).
The exterior derivative d maps a k-form to a k + 1-form. The components of dω k
for k = 0, 1, 2 are ∇ω 0 (p), ∇ × ω 1 (p), and ∇ · ω 2 (p), respectively.
The Hodge operator ⋆ maps a 1-form ω 1 to a 2-form η 2 = ⋆ω 1 . The components
2of η are η 2 (p) = A(p)ω 1 (p), where A(p) is the identity 3 × 3-matrix in the Euclidean
3-space. However, if the space was not Euclidean but a more general Riemannian
manifold (M, g), the expression for the positive definite matrix A(p) would be obtained
from the metric tensor g.
The trace operator t “restricts” a differential form to a subdomain of the Euclidean
space. Given a set of one or two tangent vectors at a point p ∈ M , the trace operator

Homology and Cohomology Computation in FE Modeling19

picks the tangential part of a 1-form with respect to the line element or the surface
element spanned by the tangent vectors. That is, the components of tω 1 are (ω 1 (p) ·û)û or ω 1 (p) − (ω 1 (p) · n̂)n̂, respectively, where the vectors û and n̂ are normalizedwith respect to the Euclidean norm. For 2-forms, the trace operator picks the normalpart with respect to the surface element spanned by two tangent vectors: tω 2 has
the components (ω 2 (p) · n̂)n̂. For a 0-form, tω 0 is the ordinary restriction of a scalarfunction.
the integration
of aP
differential
k-form ω k is defined over a k-chain ck =R
R
P Finally,kkkki ci σik tω , where each k-dimensional cell σi is a ki ci σi . It is given by ck ω =dimensional
set in the Euclidean space R3 . The Stokes’ theorem is stated as
R boundedRkkdω = ∂ck ω . This reveals that the exterior derivative d can be considered as the
ck
coboundary operator δ. Indeed, the cohomology where the cochains are differential
forms is called de Rham cohomology. Together with the Hodge operator and the wedge
product, integration provides an inner product space structure for the function
spacesR
F k (M ) of differential k-forms. The inner product is defined by hω k , η k i = M ⋆ω k ∧η k .
Thus, one can construct various Hilbert spaces of differential forms needed for arigorous treatment of differential forms in the finite element method [2].
A tensor field C that maps a 1-form ω 1 to a 1-form η 1 can be represented as a
(3 × 3)-matrix C(p). The components of η 1 are η 1 (p) = C(p)ω 1 (p). In the possible
non-linear case, the elements of the matrix C(p) depend on the components of ω 1 (p).
Then, the matrix C(p, ω 1 (p)) represents a linear approximation to the non-linear map
C : F 1 → F 1 near ω 1 (p) (in the sense of the inner product of the space F 1 ). To use
such a (possibly non-linear) map in the finite element method, the map C is required
to be strictly monotonous and symmetric with respect to the inner product of thedifferential 1-forms [18].
REFERENCES[1] D. N. Arnold, R. S. Falk, and R. Winther, Finite element exterior calculus, homological techniques, and applications, Acta Numer., 15 (2006), pp. 1–155.
doi:10.1017/S0962492906210018.
[2], Finite element exterior calculus: from Hodge theory to numerical stability, Bull. Amer.
Math. Soc. (N.S.), 47 (2010), pp. 281–354. doi:10.1090/S0273-0979-10-01278-4.
[3] A. Bossavit, Magnetostatic problems in multiply connected regions: some properties of thecurl operator, Physical Science, Measurement and Instrumentation, Management and Education - Reviews, IEE Proceedings A, 135 (1988), pp. 179 –187.
[4] A. Bossavit, Computational Electromagnetism, Academic Press, 1998.
[5] A. Bossavit, How weak is the “weak solution” in finite element methods?, Magnetics, IEEE
Transactions on, 34 (1998), pp. 2429 –2432. doi:10.1109/20.717558.
[6] CHomP Group, CHomP. http://chomp.rutgers.edu/software/, 2012. Accessed: 25.6.2013.
[7] K. Clarkson, A program for convex hulls. http://www.netlib.org/voronoi/hull.html, 1995.
Accessed: 25.6.2013.
[8] T. K. Dey, A. N. Hirani, and B. Krishnamoorthy, Optimal homologous cycles, total unimodularity, and linear programming, in Proceedings of the 42nd ACM symposium on Theoryof computing, STOC ’10, ACM, 2010, pp. 221–230. doi:10.1145/1806689.1806721.
[9] T. K. Dey, K. Li, and J. Sun, Computing geometry-aware handle and tunnel loops in 3dmodels, ACM Trans. Graph, (2008), pp. 1–9.
[10] P. Dlotko and R. Specogna, Efficient cohomology computation for electromagnetic modeling,
CMES: Computer Modeling in Engineering & Sciences, 60 (2010), pp. 247–278.
[11] P. Dlotko, R. Specogna, and F. Trevisan, Automatic generation of cuts on large-sized
meshes for the T − Ω geometric eddy-current formulation, Computational Methods in
Applied Mechanical Engineering, 198 (2009), pp. 3765–3781.
[12] P. Dular, C. Geuzaine, F. Henrotte, and W. Legros, A general environment for thetreatment of discrete problems and its application to the finite element method, Magnetics,
IEEE Transactions on, 34 (1998), pp. 3395 –3398. doi:10.1109/20.717799.

20

M. Pellikka, S. Suuriniemi, L. Kettunen, and C. Geuzaine[13] P. Dular, P. Kuo-Peng, C. Geuzaine, N. Sadowski, and J.P.A. Bastos, Dual magnetodynamic formulations and their source fields associated with massive and stranded inductors,
Magnetics, IEEE Transactions on, 36 (2000), pp. 1293 –1299. doi:10.1109/20.877677.
[14] P. Dular, W. Legros, and A. Nicolet, Coupling of local and global quantities in various finite element formulations and its application to electrostatics, magnetostatics
and magnetodynamics, Magnetics, IEEE Transactions on, 34 (1998), pp. 3078 –3081.
doi:10.1109/20.717720.
[15] H. Edelsbrunner and E. P. Mücke, Three-dimensional alpha shapes, in Proceedings of the
1992 Workshop on Volume Visualization, ACM, 1992, pp. 75–82.
[16] J. Erickson and K. Whittlesey, Greedy optimal homotopy and homology generators, in Proceedings of the Sixteenth Annual ACM-SIAM Symposium on Discrete Algorithms, SODA
’05, Society for Industrial and Applied Mathematics, 2005, pp. 1038–1046.
[17] C. Geuzaine and J.-F. Remacle, Gmsh: a three-dimensional finite element mesh generator with built-in pre- and post-processing facilities, International Journal for Numerical
Methods in Engineering, 79 (2009), pp. 1309–1331. doi:10.1002/nme.2579.
[18] P. W. Gross and P. R. Kotiuga, Electromagnetic Theory and Computation, Cambridge
University Press, 2004.
[19] X. Gu and S.-T. Yau, Global conformal surface parameterization, in Proceedings of the 2003
Eurographics/ACM SIGGRAPH symposium on Geometry processing, SGP ’03, Eurographics Association, 2003, pp. 127–137.
[20] A. Hatcher, Algebraic Topology, Cambridge University Press, 2002.
[21] F. Henrotte and K. Hameyer, An algorithm to construct the discrete cohomology basis functions required for magnetic scalar potential formulations without cuts, Magnetics, IEEE
Transactions on, 39 (2003), pp. 1167 – 1170. doi:10.1109/TMAG.2003.810404.
[22] J.-G. Dumas, F. Heckenbach, B. D. Saunders and V. Welker, GAP Homology.
http://www.eecis.udel.edu/~dumas/Homology/, 2011. Accessed: 25.6.2013.
[23] T. Kaczynski, M. Mrozek, and M. Ślusarek, Homology Computation by Reduction of Chain
Complexes, Computers Math. Applic., 35 (1998), pp. 59 – 70.
[24] L. Kettunen, K. Forsman, and A. Bossavit, Formulation of the eddy current problem in
multiply connected regions in terms of h, Int. J. Numer. Meth. Engng., 41 (1998), pp. 935– 954. doi:10.1002/(SICI)1097-0207(19980315)41:5<935::AID-NME321>3.0.CO;2-F.
[25] M. Mrozek and B. Batko, Coreduction homology algorithm, Discrete & Computational Geometry, 41 (2009), pp. 96–118. doi:10.1007/s00454-008-9073-y.
[26] J. R. Munkres, Elements of Algebraic Topology, Perseus Books, 1984.
[27] M.
Pellikka,
Magnetodynamics
withcohomology
conditions.
http://onelab.info/wiki/Magnetodynamics_with_cohomology_conditions.
Accessed:
25.6.2013.
[28] M. Pellikka, S. Suuriniemi, and L. Kettunen, Homology in Electromagnetic Boundary
Value Problems, Boundary Value Problems, 2010 (2010). Article ID 381953.
[29] M. Pellikka, S. Suuriniemi, and L. Kettunen, Powerful heuristics and basis selection
bring computational homology to engineers, Magnetics, IEEE Transactions on, 47 (2011),
pp. 1226 –1229. doi:10.1109/TMAG.2010.2091109.
[30] J.-F. Remacle, C. Geuzaine, Compère, G., and E. Marchandise, High-quality surface
remeshing using harmonic maps, Int. J. Numer. Meth. Engng., 83 (2010), pp. 403–425.
doi:10.1002/nme.282.
[31] Z. Ren, t − ω -formulation for eddy-current problems in multiply connected regions, Magnetics,
IEEE Transactions on, 38 (2002), pp. 557 –560. doi:10.1109/20.996146.
[32] R. Specogna, Complementary geometric formulations for electrostatics, International Journal
for Numerical Methods in Engineering (IJNME), 86 (2011), pp. 1041–1068.
[33] Stanford University, JPlex. http://comptop.stanford.edu/u/programs/jplex/, 2008. Accessed: 25.6.2013.
[34] S. Suuriniemi, Homological Computations in Electromagnetic Modeling, PhD thesis, Dept.
Elect. Eng., Tampere Univ. Technology, Tampere, Finland, 2004.
[35] S. Suuriniemi, J. Kangas, and L. Kettunen, Driving a coupled field-circuit problem, COMPEL: The International Journal for Computation and Mathematics in Electrical and Electronic Engineering, 26 (2007), pp. 899 – 909.
[36] C. K. Yap, Fundamental Problems of Algorithmic Algebra, Oxford University Press, 2000.
[37] A. Yavari, On geometric discretization of elasticity, Journal of Mathematical Physics, 49(2008). doi:10.1063/1.2830977.
[38] A. Zomorodian and G. Carlsson, Localized homology, Computational Geometry, 41 (2008),
pp. 126 – 148. doi:http://dx.doi.org/10.1016/j.comgeo.2008.02.003.