A Reference Discretization Strategy for the Numerical Solution of Physical Field Problems

ADVANCES IN IMAGING AND ELECTRON PHYSICS, VOL. 121
A Reference Discretization Strategy for the Numerical Solution
of Physical Field Problems
CLAUDIO MATTIUSSI*
Clampco Sistemi-NIRLAB, AREA Science Park,
Padriciano 99, 34012 Trieste, Italy
I. I ntroduction . . .
. . . . . . .
144
.
II. Foundations . . . . . . . . . . . . . . . . . .
147
A. The Mathematical Structure of Physical Field Theories
147
B. Geometric Objects and Orientation . .
150
I. Space-Time Objects . . . . . . .
155
C. Physical Laws and Physical Quantities .
157
I. Local and Global Quantities . . .
158
2. Equations. . . . . . . . . . .
159
D. Classification of Physical Quantities .
163
I. Space-Time Viewpoint
165
E. Topological Laws . . . . . .
168
F. Constitutive Relations. . . . .
172
I. Constitutive Equations and Discretization E rror .
175
G. Boundary Conditions and Sources. .
176
H. The Scope of the Structural Approach
177
III. Representations. . .
183
A. Geometry. . . .
183
I. Cell Complexes
184
2. Primary and Secondary Mesh
186
3. Incidence Numbers.
188
4. Chains . . . . . . . .
190
5. The Boundary of a Chain
191
B. Fields . . . . .
193
I. Coc hains . . .
193
2. Limit Systems.
197
C. Topological Laws
199
I. The Coboundary Operator .
200
2. Properties of the Coboundary Operator
202
3. Discrete Topological Equations.
204
D. Constitutive Relations. . .
205
E. Continuous Representations
207
I. Differential Forms
210
2. Weighted Integrals . . .
211
'Current affiliation: E volutionary and Adaptive Systems Team, Institute of Robotic Systems
(ISR), Department of Micro-E ngineering (DMT), Swiss Federal Institute of Technology (EPFL),
CH-1015 Lausanne, Switzerland.

---

1 44
CLAUDIO MATTIUSSI
3. Differential Operators . . . .
214
4. Spread Cells . . . . . . . .
217
5. Weak Form of Topological Laws
220
IV. Methods. . . . . . . . . . . . .
222
A. The Reference Discretization Strategy
222
I. Domain Discretization . . . . .
223
2. Topological Time Stepping . . .
225
3. Strategies for Constitutive Relations Discretization
231
4. Edge Elements and Field Reconstruction . .
239
B. Finite Difference Methods . . . . . . . . .
246
I. The Finite Difference Time-Domain Method
246
2. The Support Operator Method
252
3. Beyond the FDTD Method . . . .
254
C. Finite Volume Methods . . . . . . .
255
I. The Discrete Surface Integral Method
256
2. The Finite Integration Theory Method .
260
D. Finite Element Methods . . . . . . . .
264
1. Time-Domain Finite Element Methods
267
2. Time-Domain Edge Element Method .
269
3. Time-Domain Error-Based FE Method
271
V. Conclusions
273
V I. Coda . . .
275
References .
276
I. INTRODUCTION
One of the fundamental concepts of mathematical physics is that offield; that is,
naively speaking, of a spatial distribution of some mathematical object repre
senting a physical quantity. The power of this idea lies in that it allows the mod
eling of a number of very important phenomena-for example, those grouped
under the labels "electromagnetism," "thermal conduction," "fluid dynamics,"
and "solid mechanics," to name a few-and of the combinations thereof.
When the concept of field is used, a set of "translation rules" is devised,
which transforms a physical problem belonging to one of the aforementioned
domains-a physical field problem-into a mathematical one. The properties
of this mathematical model of the physical problem-a model which usually
takes the form of a set of partial differential or integrodi fferential equations,
supplemented by a set of initial and boundary conditions-can then be sub
jected to analysis in order to establish if the mathematical problem is well
posed (Gustafsson et at., 1 995). If the result of this inquiry is judged satisfac
tory, it is possible to proceed to the actual derivation of the solution, usually
with the aid of a computer.
The recourse to a computer implies, however, a further step after the model
ing step described so far, namely, the reformulation of the problem in discrete

---

NUMERICAL METHODS FOR PHYSICAL FIELD PROBLEMS
145
terms, as a tinite set of algebraic equations, which are more suitable than a
set of partial differential equations to the number-crunching capabilities of
present-day computing machines. If this discreti zation step is made by start
ing from the mathematical problem in terms of partial differential equations,
the resulting procedures can logically be called numerical methods for par tial
d ifferential equations. This is indeed how the tinite difference (FD), tinite ele
ment (FE), tinite volume (FV), and many other methods are often categorized.
Finally, the system of algebraic equations produced by the discretization step
is solved, and the result is interpreted from the point of view of the original
physical problem.
More than 30 years ago, while considering the impact of the digital computer
on mathematical activity, Bellman (1968) wrote
Much of the mathematical analysis that was developed over the eighteenth and
nineteenth centuries originated in attempts to circumvent arithmetic. With our
ability to do large-scale arithmetic . .. we can employ simple, direct methods
requiring much less old-fashioned mathematical training .. .. This situation by
no mean implies that the mathematician has been dispossessed in mathematical
physics. It does signify that he is urgently needed ... to transform the original
mathematical problems to the stage where a computer can be utilized profitably
by someone with a suitable scientific training .
. . . Good mathematics, like politics, is the an of the possible. Unfortunately,
people quickly forget the origins of a mathematical formulation with the result
that it soon acquires a life of its own. Its genealogy then protects it from scrutiny.
Because the digital computer has so greatly increased our ability to do arithmetic,
it is now imperative that we reexamine all the classical mathematical models of
mathematical physics from the standpoints of both physical significance and
feasibility of numerical solution. It may well tum out that more realistic descrip
tions are easier to handle conceptually and computationally with the aid of the
computer. (pp. 44-45)
In this spirit, the present work describes an alternative to the classical par
tial differential equations-based approach to the discretization of physical tield
problems. This alternative is based on a preliminary reformulation of the math
ematical model in a partially discrete form, which preserves as much as pos
sible the physical and geometric content of the original problem, and is made
possible by the existence and properties of a common mathematical structure
of physical tield theories (Tonti, 1 975). The goal is to maintain the focus,
both in the modeling step and in the discretization step, on the physics of
the problem, thinking in terms of numerical methods for physical field prob
lems . and not for a particular mathematical form (e.g., a partial differential
equation) into which the original physical problem happens to be translated
(Fig. I ).

---

146
CLAUDIO MATTIUSSI
physical
( field problem
"discrete,,
modeling
partially discrete
continuous
mathematical
mathematical model:
model
p.d.e.
) discretization
system of
algebraic
equations
numerical solution 1
discrete solution
I
approx. reconstruction I
...
continuous field
representation
FIGURE I. The alternative paths leading from a physical field problem to a system of alge
braic equations. p.d.e partial differential equation.
.*
The advantages of this approach are various. First, it provides a unifying
viewpoint for the discretization of physical field problems, which is valid for a
multiplicity of theories. Second, by basing the discretization of the problems
on the structural properties of the theory to which they belong, this approach
gives discrete form ulations which preserve many physically significant prop
erties of the original problem. Finally, being based on very intuitive geometric
and physical concepts, this approach facilitates both the analysis of existing
numerical methods and the development of new ones. The present work con
siders both these aspects, introducing first a reference discretization strategy
directly inspired by the res ults of the analysis of the structure of physical field
theories. Then, a n umber of popular numerical methods for partial 'differential
equations are considered, and their workings are compared with those of the
reference strategy. in order to ascertain to what extent these methods can be
interpreted as discretization methods for physical field problems.
The realization of this plan requires the preliminary introduction of the
basic ideas of the structural analysis of physical field theories. These ideas are
simple, but unfortunately they were formalized and given physically unintuitive
names at the time of their first application, within certain branches of advanced

---

NUMERICAL METHODS FOR PHYSICAL FIELD PROBLEMS
1 47
mathematics. Therefore, in applying them to other fields, one is faced with
the dilemma of inventing fo r these concepts new and, one would hope, more
meaningful names, o r maintaining the names inherited from mathematical
tradition. After some hesitation, I chose to keep the o riginal names, to avoid
a p roliferation of typically ephemeral new definitions and in consideration of
the fact that there can be difficult concepts, not difficult names; we must try to
clarify the former, not avoid the latter (Do\Cher, 1978).
The intended audi ence fo r this article is wide. On the one hand, novices to
the field of numerical methods for physical field p roblems will find herein a
framewo rk which will help them to intuitively g rasp the common concepts hid
den under the surface of a variety of methods and thus smooth the path to thei r
mastery. On the other hand, the ideas p resented should also p rove helpful to the
experienced numerical p ractitioner and to the researcher as additional tools that
can be applied to the evaluation of existing methods and the development of
new ones.
Finally, it is worth remembering that the result of the discretization must be
subjected to analysis also, in o rder to establish its p roperties as a new mathe
matical p roblem, and to measure th e effects of the discretization on the solution
when it is compared with that of nondiscrete mathematical models. This fur
ther analysis will not be dealt with he re, th e emphasis being on the unveiling of
the common discretization substratum for existing methods, the convergence,
stability, consistency, and error analyses of which abound in th e literature.
II. FOUNDATIONS
A. The Mathematical Structure of Physical Field Theories
It was mentioned in the Introduction that the approach to the disc retization
that will be p resented in this work is based on the observation that physical
field theories possess a common structure. Let us, therefo re, start by explaining
what we mean when we tal k of th e structure of a physical theo ry.
It is a common exp eri ence that exposure to more than one physical field
theo ry (e.g., thermal conduction and electrostatics) aids the comp rehension o f
each single one and facilitates th e quic k g rasping o f new ones. This occurs be
cause there are easily recognizable similarities in the mathematical fo rmulation
o f theories describing diffe rent phenomena, which p ermit the transfer of intu
ition and image ries developed fo r more familiar cases to unfamiliar realms.*
Building in a systematic way on these similarities, one can fill a correspondence
*One may say that this is the essence of explanation (i.e., the mapping of the unexplained on
something that is considered obvious),

---

148
CLAUDIO MATIIUSSI
table that relates physical quantities and laws playing a similar role within dif
ferent theories. Usually we say that there are analog ies between these theories.
These analogies are often reported as a trivial, albeit useful curiosity, but some
scholars have devoted considerable efforts to unveiling their origin and mean
ing. In these scholars' quest, they have discovered that these similarities can
be traced to the common geometric background upon which the "physics" is
built. In the book that, building on a long tradition, took these enquiries almost
to their present state, Tonti (1975) emphasized the following:
*
The existence within physical theories of a natural association of many
physical quantities, with geometric objects in space and space-time*
*
The necessity to consider as oriented the geometric objects to which phys
ical quantities are associated
*
The existence of two kinds of orientation for these geometric objects
*
The primacy and priority, in the foundation of each theory, of global phys
ical quantities associated with geometric objects, over the corresponding
densities
From this set of observations there follows naturally a classification of phys
ical quantities, based on the type and kind of orientation of the geometric
object with which they are associated. The next step is the consideration of
the relations held between physical quantities within each theory. Let us call
them generically the physical laws. From our point of view, the fundamental
observation in this context relates to
*
The existence within each theory of a set of intrinsically discrete physical
laws
These observations can be given a graphical representation as follows. A clas
s ificat ion d iagram for physical quantities is devised, with a series of "slots" for
the housing of physical quantities, each slot corresponding to a different kind
<?f oriented geometric object (see Figs. 7 and 8). The slots of this diagram can
be filled for a number of different theories. Physical laws will be represented
in this diagram as links between the slots housing the physical quantities (see
Fig. 17). The classification diagram of physical quantities, complemented by
the links representing physical laws, will be called the factor izat ion d iagram
of the physical field problem, to emphasize its role in singling out the terms
in the governing equations of a problem, according to their mathematical and
physical properties.
The classification and factorization diagrams will be used extensively in
this work. They seem to have been first introduced by Roth (see the discussion
*For the time being, we give the concept of oriented geometric object an intuitive meaning
(points, and sufficiently regular lines, surfaces, volumes, and hypervolumes, along with time
instants and time intervals).

---

NUMERICAL METHODS FOR PHYSICAL FIELD PROBLEMS
149
in Bowden, 1 990, who calls them Roth's diagrams). B ranin ( 1 966) used a
modified version of Roth 's diagrams, calling them transformation diagrams.
Tonti ( 1 975, 1 976a, 1976b, 1 998) refined and used these diagrams-which
he called classification schemes-as the basic rep resentational tool for the
analysis of the fo rmal structure of physical theories. We will refer here to
this last version of the diagrams, which were subsequently adopted by many
authors with slight g raphical variations and under various names (Baldomir
and Hammond, 1 996; Bossavit, 1 998a; Palmer and Shapiro, 1 993; Oden and
Reddy, 1 983) and for which the name Tonti diagrams was suggested. *
The Tonti classification and factorization diagrams are an ideal starting
point for the discretization of a field problem. The association of physical
quantities with geometric objects gives a rationale for the construction of the
discretization meshes and the association of the variables to the constituents of
the meshes, whereas singling out in the diagram the intrinsically discrete terms
of the field equation permits us both to pu rsue the direct discrete rendering of
these terms and to focus on the disc retization effort with the remaining terms.
Having found this common starting point for the discretization of field p rob
lems, one might be tempted to adopt a ve ry abstract viewpoint, based on a
generic field theory, with a corresponding generic terminology and factoriza
tion diagram. However, although many problems share the same structure of
the diagram, there are classes of theo ries whose diagrams differ markedly and
consequently a generic diagram would be either too simple to encompass all
the cases o r too complicated to work with. For this reason we are going to
p roceed in conc rete terms, selecting a model field theory and referring mainly
to it, in the belief that this could aid intuition, even if the reader's main in
terest is in a different field. Considering the focus of the series in which this
article appears, electromagnetism was selected as the model theory. Readers
having another background can easily translate what follows by comparing
the facto rization diagram for electromagnetism with that of the theory they are
inte rested in. To give a feeling of what is required for the development of the
factorization diagram for other theories, we discuss the case of heat transfer,
thought of as rep resentative of a class of scalar transpo rt equations.
It must be said that there are still issues that wait to be clarified in relation to
the factorization diagrams and the mathematical structure of physical theories.
This is true in particular for some issues concerning the position of energy
quantities within the diagrams and the role of o rientation with reference to
'In fact, the diagrams used in this work (and in Mattiussi, 1997) differ from those originally
conceived by Tonti in their admitting only cochains within the slots, whereas the latter had chains
in some slots and cochains in others (depending on the kind of orientation of the subjacent
geometric object). This difference reflects our advocating the use of the chain--cochain pair to
distinguish the discrete representation of the geometry (which is always made in terms of chains)
from that of the fields (which is always based on cochains).

---

1 50
CLAUDIO MATTIUSSI
time. Luckily this touches only marginally on the application of the theory to
the discretization of physical problems finalized to their numerical solution.
B. Geometric Objects and Orientation
The concept of geometric object is ubiquitous in physical field theories. For
example, in the theory of thermal conduction the heat balance equation links the
difference between the amount of heat contained inside a volume Vat the initial
and final time instants Ti and Tjof a time interval /, to the heat flowing through
the surface S, which is the boundary of V, and to the heat produced or absorbed
within the volume during the time interval. In this case, V and S are geometric
objects in space, whereas /, Ti, and Tj are geometric objects in time. The
combination of a space and a time object (e.g., the surface S considered during
the time interval /, or the volume Vat the time instant Ti, or Tj) gives a space
time geometric object. These examples show that by "geometric object" we
mean the points and the sufficiently well-behaved lines, surfaces, volumes, and
hypervolumes contained in the domain of the problem, and their combination
with time instants and time intervals. This somewhat vague definition will be
substituted later by the more detailed concept of the p-d imens ional cell .
The preceding example also shows that each mention of an object comes
with a reference to its or ientat ion . To write the heat balance equation, we must
specify if the heat flowing out of a volume or that flowing into it is to be
considered positive. This corresponds to the selection of a preferred direc
tion through the surface. * Once this direction is chosen, the surface is said to
have been given e xte rnal or ientat ion, where the qualifier "external" hints at
the fact that the orientation is specified by means of an arrow that does not
lie on the surface. Correspondingly, we will call inte rnal or ientat ion of a sur
face that which is specified by an arrow that lies on the surface and that specifies
a sense of rotation on it (Fig. 2). Note that the idea of internal orientation for
surfaces is seldom mentioned in physics but is very common in everyday ob
jects and in mathematics (Schutz, 1980). For example, a knob that must be
rotated counterclockwise to ensure a certain effect is usually designed with a
suitable curved arrow drawn on its surface, and in plane affine geometry, the
ordering of the coordinate axes corresponds to the choice of a sense of rotation
on the plane and defines the orientation of the space.
'Of course it must be possible to assign such a direction consistently, which is true if the
geometric object is orientable (Schutz, 1 980), as we will always suppose to be the case. Once
the selection is made, the object acquires a new status. As pointed out by Mac Lane (1986):
"A plane with orientation is really not the same object as one without. The plane with an orientation
has more structure-namely, the choice of the orientation" (p. 84).

---

NUMERICAL METHODS FOR PHYSICAL FIELD PROBLEMS
151
(a)
(b)
FIGURE 2. (a) External and (b) internal orientations for surfaces.
In fact, all geometric objects can be endowed with two kinds of orientations
but, for historical reasons, almost no mention of this distinction survives in
physics. * Since both kinds of orientation are needed in physics, we will show
how to build the complete orientation apparatus . We will start with internal
orientation, using the preceding affine geometry example as inspiration. An
n-dimensional affine space is oriented by fixing an order of the coordinate axes:
this, in the three-dimensional case, corresponds to the choice of a screw-sense,
or that of a vortex; in the two-dimensional case, to the choice of a sense of
rotation on the plane; and in the one-dimensional case, to the choice of a sense
(an arrow) along the line. These images can be extended to geometric objects.
Therefore, the internal orientation of a volume is given by a screw-sense; that
of a surface, by a sense of rotation on it; and that of a line, by a sense along it
(see Fig. 5).
Before we proceed further, it is instructive to consider an example of a
physical quantity that, contrary to common belief, is associated with internally
oriented surfaces: the magnetic flux <p. This association is a consequence of
the invariance requirement of Maxwell's equations for improper coordinate
transformations; that is, those that invert the orientation of space, transforming
a right-handed reference system into a left-handed one. Imagine an experi
mental setup to probe Faraday's law, for example, verifying the link between
the magnetic flux <p "through" a disk S and the circulation U of the electric
field intensity E around the loop r which is the border of S. If we suppose, as
is usually the case, that the sign of <p is determined by a direction through the
disk, and that of U by the choice of a sense around the loop, a mirror reflection
through a plane parallel to the disk axis changes the sign of U but not that of
<p. Usually the incongruence is avoided by using the right-hand rule to define
B and invoking for it the status of axial vector (Jackson, 1 975). In other words,
we are told that for space reflections, the sense of the "arrow" of the B vector
'However, for example, Maxwell ( 1 87 1 ) was well aware of the necessity within the context
of electromagnetism of at least four kinds of mathematical entities for the correct representation
of the electromagnetic field (entities referred to lines or to surfaces and endowed with internal or
with external orientation).

---

1 52
CLAUDIO MATTIUSSI
(a)
(b)
FIGURE 3. Orientational issues in Faraday's law. The intervention of the right-hand rule,
required in the classical version (a), can be avoided by endowing both geometric objects rand
S with the same kind of orientation (b).
does not count; only the right-hand rule does. It is, however, apparent that for
the invariance of Faraday's law to hold true without such tricks, all we have to
do is either to associate with internally oriented surfaces and U with inter
nally oriented lines, or to associate with externally oriented surfaces and U
with lines oriented by a sense of rotation around them (i.e., externally oriented
lines, as will soon be clear). Since the effects of an electric field act along the
field lines and not around them, the first option seems preferable (Schouten,
1 989; Fig. 3).
This example shows that the need for the right -hand rule is a consequence of
our disregarding the existence of two kinds of orientation. This attitude seems
reasonable in physics as we have become accustomed to it in the course of
our education, but consider that if it were applied systematically to everyday
objects, we would be forced to glue an arrow pointing outward from the afore
mentioned knob, and to accompany it with a description of the right-hand
rule. Note also that the difficulties in the classical formulation of Faraday's
law stem from the impossibility of comparing directly the orientation of the
surface with that of its boundary, when the surface is externally oriented and
the bounding line is internally oriented. In this case, "directly" means "without
recourse to the right-hand rule" or similar tricks. The possibility of making
this direct comparison is fundamental for the correct statement of many phys
ical laws. This comparison is based on the idea of an orientation induced by
an object on its boundary. For example, the sense of rotation that internally
orients a surface induces a sense of rotation on its bounding curve, which can
be compared with the sense of rotation which orients the surface internally.
The same is true for the internal orientation of volumes and of their bounding
surfaces. The reader can check that the direct comparison is indeed possible
if the object and its boundary are both endowed with internal orientation as
defined previously for volumes, surfaces, and lines. However, this raises an
interesting issue, since our list of internally oriented objects does not so far
include points, which nevertheless form the boundary of a line. To make inner
orientation a coherent system, we must, therefore, define internal orientations
for points (as in algebra we extend the definition of the nth power of a number
to include the case n = 0). This can be done by means of a pair of symbols

---