You do not have permission to edit this page, for the following reason:

You are not allowed to execute the action you have requested.


You can view and copy the source of this page.

x
 
1
==Abstract==
2
3
We revisit the theory of Discrete Exterior Calculus (DEC) in 2D for general triangulations, relying only on Vector Calculus and Matrix Algebra.  We present DEC numerical solutions of the Poisson equation and compare them against those  found using the Finite Element Method with linear elements (FEML).
4
5
'''Keywords''': Discrete element method, Poisson equation
6
7
==1 Introduction==
8
9
The purpose of this paper is to introduce the theory of Discrete Exterior Calculus (DEC)  to the widest possible audience and, therefore, we will  rely mainly on Vector Calculus and Matrix Algebra. Discrete Exterior Calculus is a relatively new method for solving partial differential equations <span id='citeF-8'></span>[[#cite-8|[8]]] based on the idea of discretizing the mathematical theory of Exterior Differential Calculus, a theory that goes back to E. Cartan <span id='citeF-3'></span>[[#cite-3|[3]]] and is fundamental  in the areas of Differential Geometry and Differential Topology. Although Exterior Differential Calculus is an abstract mathematical theory, it has been introduced in various fields such as in digital geometry processing <span id='citeF-4'></span>[[#cite-4|[4]]],  numerical schemes for partial differential equations <span id='citeF-8'></span><span id='citeF-1'></span>[[#cite-8|[8,1]]], etc.
10
11
In his PhD thesis <span id='citeF-8'></span>[[#cite-8|[8]]], Hirani laid down the fundamental concepts  of Discrete Exterior Calculus (DEC), using discrete combinatorial and geometric  operations on simplicial complexes (in any dimension), proposing discrete equivalents for differential forms, vector fields, differential and geometric operators, etc. Perhaps the first numerical application of DEC to PDE was given in <span id='citeF-9'></span>[[#cite-9|[9]]] in order to solve Darcy flow and Poisson's equation.  In <span id='citeF-7'></span>[[#cite-7|[7]]], the authors develop a modification of DEC and show that in simple cases (e.g. flat geometry and regular meshes), the equations resulting from DEC are equivalent to classical numerical schemes such as finite difference or finite volume discretizations. In <span id='citeF-11'></span>[[#cite-11|[11]]], the authors used DEC to solve the Navier-Stokes equations and,  in <span id='citeF-5'></span>[[#cite-5|[5]]] DEC was used with a discrete lattice model to simulate  elasticity, plasticity and failure of isotropic materials.
12
13
In this expository paper, we review the various operators of Exterior Differential Calculus in 2D in terms of ordinary vector calculus, and introduce only the geometrical ideas that are essential  to the formulation.  Among those ideas is that of duality between the differentiation operator (on vector fields) and the boundary operator (on the domain) contained in Green's theorem.  This duality is one of the key ideas of the method, which justifies taking the discretized derivative matrix as the transpose of the boundary operator matrix on the given mesh.   Another important ingredient is the Hodge star operator, which is hidden in the notation of Vector Calculus. In order to show the necessity of the Hodge star operator,  we carry out some simple calculations. In particular, we will introduce the notion of wedge product of vectors which, roughly speaking, helps us assign algebraic objects to parallelograms and carry out algebraic manipulations with them. We present  DEC in the simplest terms possible using easy examples. We also review the formulation of DEC for arbitrary meshes, which was first considered in <span id='citeF-10'></span>[[#cite-10|[10]]]. Performance of the method is tested on the Poisson equation and compared with the Finite Element Method with linear elements (FEML).
14
15
The paper is organized as follows.  In Section [[#2 2D Exterior Differential Calculus as Vector Calculus|2]], we introduce the wedge product of vectors and the ''geometric'' Hodge star operator, and rewrite Green's theorem appropriately in order to display the duality between the differentiation and the boundary operators. In Section [[#3 Discrete Exterior Calculus|3]], we present the operators of DEC (mesh, dual mesh, discrete derivation, discrete Hodge star operator), showing simple examples throughout. In Section [[#4 DEC for general triangulations|4]], we present  the formulation of DEC on arbitrary triangulations. In Section [[#5 Numerical examples|5]], we present the numerical solution of a Poisson equation with DEC and FEML, in order to compare their performance. In Section [[#6 Conclusions|6]], we present our conclusions.
16
17
==2 2D Exterior Differential Calculus as Vector Calculus==
18
19
In this section we introduce two geometric operators (the wedge product and the Hodge star) and explain how to use them together with the gradient operator in order to obtain the Laplacian.
20
21
===2.1 <span id='lb-2.1'></span>Wedge product for vectors in R²===
22
23
Let <math display="inline">a,b</math> be vectors in <math display="inline">\mathbb{R}^2</math>. We can assign to these vector the parallelogram they span, and to such a parallelogram its area. The latter is equal to the determinant of the transformation matrix  sending <math display="inline">e_1,e_2</math> to <math display="inline">a,b</math> respectively. In Exterior Calculus, such a parallelogram is regarded as an algebraic object <math display="inline">a\wedge b</math>, a bivector, and the set of bivectors is equipped with a vector space structure to form the one-dimensional vector space
24
25
{| class="formulaSCP" style="width: 100%; text-align: left;" 
26
|-
27
| 
28
{| style="text-align: left; margin:auto;width: 100%;" 
29
|-
30
| style="text-align: center;" | <math>\textstyle \wedge ^2\mathbb{R}^2</math>
31
|}
32
|}
33
34
Thus, we have the following spaces
35
36
{| class="formulaSCP" style="width: 100%; text-align: left;" 
37
|-
38
| 
39
{| style="text-align: left; margin:auto;width: 100%;" 
40
|-
41
| style="text-align: center;" | <math>\mathbb{R}, \quad \quad \mathbb{R}^2, \quad \quad  \wedge^2\mathbb{R}^2</math>
42
|}
43
|}
44
45
of scalars, vectors and bivectors, respectively. For instance, the square formed by <math display="inline">e_1</math> and <math display="inline">e_2</math> is represented by the symbol
46
47
{| class="formulaSCP" style="width: 100%; text-align: left;" 
48
|-
49
| 
50
{| style="text-align: left; margin:auto;width: 100%;" 
51
|-
52
| style="text-align: center;" | <math>e_1\wedge e_2</math>
53
|}
54
|}
55
56
which is read as "<math display="inline">e_1</math> wedge <math display="inline">e_2</math>" (see Figure [[#lb-2.1|2.1]])
57
58
{| class="floating_imageSCP" style="text-align: center; "
59
|-
60
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
61
62
<span style="text-align: center; font-size: 75%;">'' [[Image:Review_239195285364_6178_text-Square00-eps-converted-to.png|180px|Wedge product of two vectors.]] ''</span></div>
63
|}
64
65
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
66
<span style="text-align: center; font-size: 75%;">'''Figure 2.1:''' Wedge product of two vectors.</span></div>
67
68
69
In <math display="inline">\mathbb{R}^2</math>, this represents an “element” of unit area.
70
71
Note that if we list the vectors in the opposite order, we have a different orientation and, therefore, the algebraic objects must satisfy <math display="inline">e_2\wedge e_1=-e_1\wedge e_2</math> (see Figure [[#lb-2.1|2.2]])
72
73
{| class="floating_imageSCP" style="text-align: center; "
74
|-
75
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
76
77
<span style="text-align: center; font-size: 75%;">''[[Image:Review_239195285364_6648_text-Square01-eps-converted-to.png|250px|Change of orientation of the parallelogram implies anticommutativity in the wedge product of two vectors.]]''</span></div>
78
|}
79
80
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
81
<span style="text-align: center; font-size: 75%;">'''Figure 2.2:''' Change of orientation of the parallelogram implies anticommutativity in the wedge product of two vectors.</span></div>
82
83
84
85
More generally, given two vectors <math display="inline">a,b\in \mathbb{R}^2</math>, their wedge product <math display="inline">a\wedge b</math> looks as follows (see Figure [[#lb-2.1|2.3]])
86
87
{| class="floating_imageSCP" style="text-align: center; "
88
|-
89
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
90
91
<span style="text-align: center; font-size: 75%;">''[[Image:Review_239195285364_4688_text-Parallelogram-eps-converted-to.png|220px|The wedge product of the vectors <math>a</math> and <math>b</math>.]]''</span></div>
92
|}
93
94
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
95
<span style="text-align: center; font-size: 75%;">'''Figure 2.3:''' The wedge product of the vectors <math>a</math> and <math>b</math>.</span></div>
96
97
98
99
The properties of the wedge product are
100
101
* it is anticommutative: <math display="inline">a\wedge b = - b\wedge a</math> (see Figure [[#lb-2.1|2.4]])
102
103
{| class="floating_imageSCP" style="text-align: center; "
104
|-
105
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
106
107
<span style="text-align: center; font-size: 75%;">[[Image:Review_239195285364_4197_text-ParalPos-eps-converted-to.png|220px|]] [[File:Review_239195285364_6087_text-ParalNeg-eps-converted-to.png|220px|]]</span></div>
108
|}
109
110
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
111
<span style="text-align: center; font-size: 75%;">'''Figure 2.4:''' Anticommutativity of the wedge product.</span></div>
112
113
114
* <math display="inline">a\wedge a = 0</math> since it is a parallelogram with area zero (see Figure [[#lb-2.1|2.5]])
115
116
{| class="floating_imageSCP" style="text-align: center; "
117
|-
118
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
119
<span style="text-align: center; font-size: 75%;">[[Image:Review_239195285364_7320_text-ParalZero-eps-converted-to.png|180px|Parallelogram with very small area, depicting what happens when <math>b</math> tends to <math>a</math>.]]</span></div>
120
|} 
121
122
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
123
<span style="text-align: center; font-size: 75%;">'''Figure 2.5:''' Parallelogram with very small area, depicting what happens when <math>b</math> tends to <math>a</math>.</span></div>
124
125
126
* it is distributive
127
128
{| class="formulaSCP" style="width: 100%; text-align: left;" 
129
|-
130
| 
131
{| style="text-align: left; margin:auto;width: 100%;" 
132
|-
133
| style="text-align: center;" | <math>(a + b )\wedge c = a\wedge c + b\wedge c</math>
134
|}
135
|}
136
* it is associative
137
138
{| class="formulaSCP" style="width: 100%; text-align: left;" 
139
|-
140
| 
141
{| style="text-align: left; margin:auto;width: 100%;" 
142
|-
143
| style="text-align: center;" | <math>(a \wedge b )\wedge c = a\wedge (b \wedge  c)</math>
144
|}
145
|}
146
147
For example, let
148
149
{| class="formulaSCP" style="width: 100%; text-align: left;" 
150
|-
151
| 
152
{| style="text-align: left; margin:auto;width: 100%;" 
153
|-
154
| style="text-align: center;" | <math>a= (a_1,a_2)\,\, = a_1 e_1 + a_2 e_2</math>
155
|-
156
| style="text-align: center;" | <math>  b= (b_1,b_2)\,\, \,= b_1 e_1 + b_2 e_2 </math>
157
|}
158
|}
159
160
Then
161
162
{| class="formulaSCP" style="width: 100%; text-align: left;" 
163
|-
164
| 
165
{| style="text-align: left; margin:auto;width: 100%;" 
166
|-
167
| style="text-align: center;" | <math>a\wedge b      =(a_1 e_1 + a_2 e_2)\wedge (b_1 e_1 + b_2 e_2)</math>
168
|-
169
| style="text-align: center;" | <math>     =a_1b_1\, e_1\wedge e_1 + a_1b_2\, e_1\wedge e_2 + a_2b_1\, e_2\wedge e_1 + a_2b_2\, e_2\wedge e_2 </math>
170
|-
171
| style="text-align: center;" | <math>     =a_1b_2\, e_1\wedge e_2 - a_2b_1\, e_1\wedge e_2  </math>
172
|-
173
| style="text-align: center;" | <math>     =(a_1b_2 - a_2b_1)\, e_1\wedge e_2  </math>
174
|}
175
|}
176
177
i.e. the determinant
178
179
{| class="formulaSCP" style="width: 100%; text-align: left;" 
180
|-
181
| 
182
{| style="text-align: left; margin:auto;width: 100%;" 
183
|-
184
| style="text-align: center;" | <math>\det \left(\begin{array}{cc} a_1 & a_2 \\  b_1 & b_2 \end{array} \right)</math>
185
|}
186
|}
187
188
times the canonical bivector/square <math display="inline">e_1\wedge e_2</math>.
189
190
===2.2 <span id='lb-2.2'></span>Hodge star operator for vectors in R²===
191
192
Now consider the following situation: given the vector <math display="inline">e_1</math>, find another vector <math display="inline">v</math> such that the parallelogram that they form has area <math display="inline">1</math>. It is readily seen that, for instance, <math display="inline">v=e_2, -e_2,e_1+e_2</math> are all solutions. Requiring orthogonality and standard orientation, we see that <math display="inline">e_2</math> is the unique solution. This process is summarized in the Hodge star operator, which basically says that the complementary vector for <math display="inline">e_1</math> is <math display="inline">e_2</math>, and the one for <math display="inline">e_2</math> is <math display="inline">-e_1</math>,
193
194
{| class="formulaSCP" style="width: 100%; text-align: left;" 
195
|-
196
| 
197
{| style="text-align: left; margin:auto;width: 100%;" 
198
|-
199
| style="text-align: center;" | <math>\star e_1 = e_2</math>
200
|-
201
| style="text-align: center;" | <math> {}\star e_2 = -e_1 </math>
202
|}
203
|}
204
205
In general, the equation that defines the Hodge star operator for any given vector <math display="inline">v\in \mathbb{R}^2</math> is the following
206
207
{| class="formulaSCP" style="width: 100%; text-align: left;" 
208
|-
209
| 
210
{| style="text-align: left; margin:auto;width: 100%;" 
211
|-
212
| style="text-align: center;" | <math>w\wedge (\star v) = (w\cdot v) \, e_1\wedge e_2</math>
213
|}
214
|}
215
216
for every <math display="inline">w\in \mathbb{R}^2</math>. In particular, if we take <math display="inline">v=w</math>,
217
218
{| class="formulaSCP" style="width: 100%; text-align: left;" 
219
|-
220
| 
221
{| style="text-align: left; margin:auto;width: 100%;" 
222
|-
223
| style="text-align: center;" | <math>v\wedge (\star v) = |v|^2 \, e_1\wedge e_2</math>
224
|}
225
|}
226
227
which means that <math display="inline">v</math> and <math display="inline">\star v</math> form a square of area <math display="inline">|v|^2</math>. Thus, the Hodge star operator on a vector <math display="inline">v\in \mathbb{R}^2</math> can be thought of as finding the vector that makes with <math display="inline">v</math> a positively oriented square of area <math display="inline">\vert v\vert ^2</math>.
228
229
It is somewhat less intuitive to work out the Hodge star of a bivector. First of all, we have to treat bivectors as vectors of a different space, namely the space of bivectors <math display="inline">\textstyle \wedge ^2\mathbb{R}^2</math>. Secondly, the length of a bivector <math display="inline">v\wedge w</math> is its area
230
231
{| class="formulaSCP" style="width: 100%; text-align: left;" 
232
|-
233
| 
234
{| style="text-align: left; margin:auto;width: 100%;" 
235
|-
236
| style="text-align: center;" | <math>{length}(v\wedge w):= {Area}(v\wedge w)</math>
237
|}
238
|}
239
240
Thus, the defining equation of the Hodge star applied to <math display="inline">(v\wedge w)</math> and <math display="inline">\star (v\wedge w)</math> reads as follows
241
242
{| class="formulaSCP" style="width: 100%; text-align: left;" 
243
|-
244
| 
245
{| style="text-align: left; margin:auto;width: 100%;" 
246
|-
247
| style="text-align: center;" | <math>(v\wedge w)\wedge \star (v\wedge w)= \left\langle v\wedge w, v\wedge w\right\rangle e_1\wedge e_2 </math>
248
|}
249
|}
250
251
which means
252
253
{| class="formulaSCP" style="width: 100%; text-align: left;" 
254
|-
255
| 
256
{| style="text-align: left; margin:auto;width: 100%;" 
257
|-
258
| style="text-align: center;" | <math>(v\wedge w)\wedge \star (v\wedge w) = {Area}(v\wedge w)^2 e_1\wedge e_2 </math>
259
|}
260
|}
261
262
Since <math display="inline">v\wedge w={Area}(v\wedge w) e_1\wedge e_2</math> is already a bivector, <math display="inline">\star (v\wedge w)</math> must be a scalar, i.e.
263
264
{| class="formulaSCP" style="width: 100%; text-align: left;" 
265
|-
266
| 
267
{| style="text-align: left; margin:auto;width: 100%;" 
268
|-
269
| style="text-align: center;" | <math>\star (v\wedge w)={Area}(v\wedge w)</math>
270
|}
271
|}
272
273
When <math display="inline">v=e_1</math> and <math display="inline">w=e_2</math> we have
274
275
{| class="formulaSCP" style="width: 100%; text-align: left;" 
276
|-
277
| 
278
{| style="text-align: left; margin:auto;width: 100%;" 
279
|-
280
| style="text-align: center;" | <math>\star (e_1\wedge e_2) = 1</math>
281
|}
282
|}
283
284
Finally, the Hodge star of a number <math display="inline">\lambda </math> is a bivector, i.e.
285
286
{| class="formulaSCP" style="width: 100%; text-align: left;" 
287
|-
288
| 
289
{| style="text-align: left; margin:auto;width: 100%;" 
290
|-
291
| style="text-align: center;" | <math>\star \lambda =\lambda \, e_1\wedge e_2</math>
292
|}
293
|}
294
295
===2.3 The Laplacian===
296
297
Let <math display="inline">f:\mathbb{R}^2\rightarrow \mathbb{R}</math> and consider the gradient
298
299
{| class="formulaSCP" style="width: 100%; text-align: left;" 
300
|-
301
| 
302
{| style="text-align: left; margin:auto;width: 100%;" 
303
|-
304
| style="text-align: center;" | <math>\nabla f = {\partial f\over \partial x}e_1+{\partial f\over \partial y}e_2</math>
305
|}
306
|}
307
308
Apply the Hodge star operator to it
309
310
{| class="formulaSCP" style="width: 100%; text-align: left;" 
311
|-
312
| 
313
{| style="text-align: left; margin:auto;width: 100%;" 
314
|-
315
| style="text-align: center;" | <math>\star \nabla f     = {\partial f\over \partial x} \star e_1+{\partial f\over \partial y} \star e_2</math>
316
|-
317
| style="text-align: center;" | <math>    = {\partial f\over \partial x}e_2-{\partial f\over \partial y}e_1</math>
318
|}
319
|}
320
321
Now take the gradient of each coefficient function together with wedge product
322
323
{| class="formulaSCP" style="width: 100%; text-align: left;" 
324
|-
325
| 
326
{| style="text-align: left; margin:auto;width: 100%;" 
327
|-
328
| style="text-align: center;" | <math>\nabla^\wedge \star \nabla f     := \nabla \left({\partial f\over \partial x}\right)\wedge e_2-\nabla \left({\partial f\over \partial y}\right)\wedge e_1</math>
329
|-
330
| style="text-align: center;" | <math>    = \left({\partial ^2 f\over \partial x^2}e_1+{\partial ^2 f\over \partial y\partial x}e_2\right)\wedge e_2    -\left({\partial ^2 f\over \partial x\partial y}e_1+{\partial ^2  f\over \partial ^2 y}e_2\right)\wedge e_1</math>
331
|-
332
| style="text-align: center;" | <math>    = {\partial ^2 f\over \partial x^2}e_1\wedge e_2+{\partial ^2 f\over \partial y\partial x}e_2\wedge e_2    -{\partial ^2 f\over \partial x\partial y}e_1\wedge e_1-{\partial ^2  f\over \partial ^2 y}e_2\wedge e_1</math>
333
|-
334
| style="text-align: center;" | <math>    = \left({\partial ^2 f\over \partial x^2}+{\partial ^2  f\over \partial ^2 y}\right)e_1\wedge e_2 </math>
335
|}
336
|}
337
338
By taking the Hodge star of this last expression we get the ordinary Laplacian of <math display="inline">f</math>
339
340
{| class="formulaSCP" style="width: 100%; text-align: left;" 
341
|-
342
| 
343
{| style="text-align: left; margin:auto;width: 100%;" 
344
|-
345
| style="text-align: center;" | <math>\star \nabla ^\wedge \star \nabla (f)     = {\partial ^2 f\over \partial x^2}+{\partial ^2  f\over \partial ^2 y} </math>
346
|}
347
|}
348
349
'''Remarks. ''' 
350
351
(i) This rather convoluted looking way of computing the Laplacian of a function is based on Exterior Differential Calculus, a theory that generalizes  the operators of vector calculus (gradient, curl and divergence) to arbitrary dimensions, and is the basis for the differential topological theory of deRham cohomology.
352
353
(ii) We would like to emphasize the  necessity of using the Hodge star operator <math display="inline">\star</math> in order to make the combination of differentiation  and wedge product produce the correct answer.
354
355
===2.4 Duality in Green's theorem===
356
357
Green's theorem states that for a vector field <math display="inline">(L,M)</math> defined on a region <math display="inline">D\subset  \mathbb{R}^2</math>,
358
359
{| class="formulaSCP" style="width: 100%; text-align: left;" 
360
|-
361
| 
362
{| style="text-align: left; margin:auto;width: 100%;" 
363
|-
364
| style="text-align: center;" | <math>\int _D\left({\partial M\over \partial x}-{\partial L\over \partial y}\right)dxdy=\int _{C=\partial D} (Ldx+Mdy) .</math>
365
|}
366
|}
367
368
In this section, we will explain how this identity encodes a duality between the operator of differentiation and that of taking the boundary of the domain of integration.
369
370
====2.4.1 Rewriting Green's theorem====
371
372
Note that if <math display="inline">F=(L,M)</math> is a vector field, we can write it as
373
374
{| class="formulaSCP" style="width: 100%; text-align: left;" 
375
|-
376
| 
377
{| style="text-align: left; margin:auto;width: 100%;" 
378
|-
379
| style="text-align: center;" | <math>F= Le_1+Me_2</math>
380
|}
381
|}
382
383
Then, we can apply the gradient operator together with wedge product in the following fashion
384
385
{| class="formulaSCP" style="width: 100%; text-align: left;" 
386
|-
387
| 
388
{| style="text-align: left; margin:auto;width: 100%;" 
389
|-
390
| style="text-align: center;" | <math>\nabla ^\wedge (L e_1 +  M  e_2)    :=\nabla L\wedge e_1 + \nabla M \wedge e_2</math>
391
|-
392
| style="text-align: center;" | <math>    =\left({\partial L\over \partial x}e_1+{\partial L\over \partial y}e_2\right)\wedge e_1   +\left({\partial M\over \partial x}e_1+{\partial M\over \partial y}e_2\right)\wedge e_2</math>
393
|-
394
| style="text-align: center;" | <math>    ={\partial L\over \partial x}e_1\wedge e_1+{\partial L\over \partial y}e_2\wedge e_1   +{\partial M\over \partial x}e_1\wedge e_2+{\partial M\over \partial y}e_2\wedge e_2</math>
395
|-
396
| style="text-align: center;" | <math>      ={\partial L\over \partial y}e_2\wedge e_1   +{\partial M\over \partial x}e_1\wedge e_2</math>
397
|-
398
| style="text-align: center;" | <math>      =\left({\partial M\over \partial x}   -{\partial L\over \partial y}\right)e_1\wedge e_2 </math>
399
|}
400
|}
401
402
Note that we have defined a new operator <math display="inline">\nabla ^\wedge </math> which combines differentiation and wedge product. Applying the Hodge star operator we obtain
403
404
{| class="formulaSCP" style="width: 100%; text-align: left;" 
405
|-
406
| 
407
{| style="text-align: left; margin:auto;width: 100%;" 
408
|-
409
| style="text-align: center;" | <math>\star \nabla ^\wedge (L e_1 +  M e_2)    =\left({\partial M\over \partial x}   -{\partial L\over \partial y}\right) </math>
410
|}
411
|}
412
413
i.e. the integrand of the left hand side of the identity of integrals in Green's theorem. Thus, as a first step, Green's theorem can be rewritten as follows:
414
415
{| class="formulaSCP" style="width: 100%; text-align: left;" 
416
|-
417
| 
418
{| style="text-align: left; margin:auto;width: 100%;" 
419
|-
420
| style="text-align: center;" | <math>\int _D\nabla ^\wedge (L, M)dxdy=\int _{\partial D} (L ,  M )\cdot (dx,dy)</math>
421
|}
422
|}
423
424
====2.4.2 The duality in Green's theorem====
425
426
Let us recall the following fact from Linear Algebra. Given a linear transformation <math display="inline">A</math> of Euclidean space <math display="inline">\mathbb{R}^n</math>, the transpose <math display="inline">A^T</math> satisfies
427
428
{| class="formulaSCP" style="width: 100%; text-align: left;" 
429
|-
430
| 
431
{| style="text-align: left; margin:auto;width: 100%;" 
432
|-
433
| style="text-align: center;" | <math>\left\langle A(v), w \right\rangle= \left\langle v, A^T(w)\right\rangle</math>
434
|}
435
|}
436
437
for any two vectors <math display="inline">v,w\in \mathbb{R}^n</math>, where <math display="inline">\left\langle\cdot , \cdot \right\rangle</math> denotes the standard inner/dot product.  In fact, such an identity characterizes the transpose <math display="inline">A^T</math>.
438
439
Now, let us do the following notational trick: substitute the integration symbols in Green's theorem by <math display="inline">\ll \cdot ,\cdot \gg </math> as follows:
440
441
{| class="formulaSCP" style="width: 100%; text-align: left;" 
442
|-
443
| 
444
{| style="text-align: left; margin:auto;width: 100%;" 
445
|-
446
| style="text-align: center;" | <math>\int _C (L,M)\cdot (dx,dy) = \ll (L,M), C\gg </math>
447
|-
448
| style="text-align: center;" | <math> \int _D\nabla ^\wedge (L, M)dxdy = \ll \nabla ^\wedge (L,M), D\gg  </math>
449
|}
450
|}
451
452
where <math display="inline">C=\partial D</math> is the boundary of the region. Using this notational change, Green's theorem reads as follows
453
454
{| class="formulaSCP" style="width: 100%; text-align: left;" 
455
|-
456
| 
457
{| style="text-align: left; margin:auto;width: 100%;" 
458
|-
459
| style="text-align: center;" | <math>\ll \nabla ^\wedge (L,M), D\gg \,\,=\,\,\ll (L,M), \partial D\gg</math>
460
|}
461
|}
462
463
Roughly speaking, this means that the differential operator <math display="inline">\nabla ^\wedge </math> is the transpose of the boundary operator <math display="inline">\partial </math> by means of the product <math display="inline">\ll \cdot ,\cdot \gg </math>.
464
465
'''Remark'''. The previous observation is fundamental in the development of DEC, since the boundary operator is well understood and easy to calculate on meshes.
466
467
==3 Discrete Exterior Calculus==
468
469
Now we will discretize the differentiation operator <math display="inline">\nabla ^\wedge </math> presented above. We will start by describing the discrete version of the boundary operator on simplices/triangles. Afterwards, we will treat the differentiation operator as the transpose of the boundary operator.
470
471
We are interested in using certain geometric subsets of a given triangular mesh of a 2D region. Such subsets include vertices/nodes, edges/sides and faces/triangles. We will describe each one of them by means of the ordered list of vertices whose convex closure constitutes the subset of interest. For instance, consider the triangular mesh of the planar hexagonal region in Figure [[#3 Discrete Exterior Calculus|3.1]],
472
473
{| class="floating_imageSCP" style="text-align: center; "
474
|-
475
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
476
<span style="text-align: center; font-size: 75%;">'' [[Image:Review_239195285364_2873_text-Region-01-eps-converted-to.png|180px|]]''</span></div>
477
|}
478
479
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
480
<span style="text-align: center; font-size: 75%;">'''Figure 3.1:''' Triangular mesh of a planar hexagonal region.</span></div>
481
482
483
where the shaded triangle will be denoted by <math display="inline">[v_0,v_1,v_6]</math>, and its edge joining the vertices <math display="inline">v_0</math> and <math display="inline">v_1</math> will be denoted by <math display="inline">[v_0,v_1]</math>. For the sake of notational consistency, we will denote the vertices also enclosed in brackets, e.g. <math display="inline">[v_0]</math>.
484
485
===3.1 Boundary operator and discrete derivative===
486
487
There is a well known boundary operator <math display="inline">\partial </math> for oriented triangles, edges and points:
488
489
* For points/vertices:
490
491
{| class="floating_imageSCP" style="text-align: center; "
492
|-
493
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
494
<span style="text-align: center; font-size: 75%;">''[[Image:Review_239195285364_9903_text-boundary00-eps-converted-to.png|80px|]]''</span></div>
495
|}
496
497
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
498
<span style="text-align: center; font-size: 75%;">'''Figure 3.2:''' Boundary of a vertex: <math>\partial [v_0] = 0</math>.</span></div>
499
500
501
* For sides/edges:
502
503
{| class="floating_imageSCP" style="text-align: center; "
504
|-
505
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
506
<span style="text-align: center; font-size: 75%;">''[[Image:Review_239195285364_1507_text-boundary01-eps-converted-to.png|200px|]]''</span></div>
507
|}
508
509
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
510
<span style="text-align: center; font-size: 75%;">'''Figure 3.3:''' Boundary of an edge: <math>\partial [v_0,v_1] = [v_1] - [v_0]</math>.</span></div>
511
512
513
* For faces/triangles:
514
515
{| class="floating_imageSCP" style="text-align: center; "
516
|-
517
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
518
<span style="text-align: center; font-size: 75%;">''[[Image:Review_239195285364_1837_text-boundary02-eps-converted-to.png|200px|]]''</span></div>
519
|}
520
521
522
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
523
<span style="text-align: center; font-size: 75%;">'''Figure 3.4:''' Boundary of a face: <math>\partial [v_0,v_1,v_2] = [v_1,v_2] - [v_0,v_2]+[v_0,v_1]</math>.</span></div>
524
525
526
'''Example'''. Let us consider again the mesh of the planar hexagonal (with oriented triangles) in Figure [[#3.1 Boundary operator|3.5]]
527
528
{| class="floating_imageSCP" style="text-align: center; "
529
|-
530
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
531
<span style="text-align: center; font-size: 75%;">''[[Image:Review_239195285364_6355_text-boundary04-eps-converted-to.png|250px|]]''</span></div>
532
|}
533
 
534
535
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
536
<span style="text-align: center; font-size: 75%;">'''Figure 3.5:''' Oriented triangular mesh of a planar hexagonal region.</span></div>
537
538
539
We will denote a triangle by the list of its vertices listed in order according to the orientation of the tringle. Thus, we have the following ordered lists:
540
541
* list of faces
542
543
{| class="formulaSCP" style="width: 100%; text-align: left;" 
544
|-
545
| 
546
{| style="text-align: left; margin:auto;width: 100%;" 
547
|-
548
| style="text-align: center;" | <math> \{[v_0,v_1,v_6],[v_1,v_2,v_6],[v_2,v_3,v_6], [v_3,v_4,v_6],[v_4,v_5,v_6],[v_5,v_0,v_6]\}  </math>
549
|}
550
|}
551
* list of edges
552
553
{| class="formulaSCP" style="width: 100%; text-align: left;" 
554
|-
555
| 
556
{| style="text-align: left; margin:auto;width: 100%;" 
557
|-
558
| style="text-align: center;" | <math> \{  [v_0,v_6], [v_1,v_6], [v_2,v_6], [v_3,v_6], [v_4,v_6], [v_5,v_6], [v_0,v_1], [v_1,v_2], [v_2,v_3], [v_3,v_4], [v_4,v_5], [v_5,v_0] \} </math>
559
|}
560
|}
561
* list of vertices
562
563
{| class="formulaSCP" style="width: 100%; text-align: left;" 
564
|-
565
| 
566
{| style="text-align: left; margin:auto;width: 100%;" 
567
|-
568
| style="text-align: center;" | <math> \{  [v_0],[v_1],[v_2],[v_3],[v_4],[v_5],[v_6] \} </math>
569
|}
570
|}
571
572
A key idea in DEC is to consider each face as an element of a basis of a vector space. Namely,  column coordinate vectors are associated to faces as follows:
573
574
{| class="formulaSCP" style="width: 100%; text-align: left;" 
575
|-
576
| 
577
{| style="text-align: left; margin:auto;width: 100%;" 
578
|-
579
| style="text-align: center;" | <math>[v_0,v_1,v_6] \longleftrightarrow  (1,0,0,0,0,0)^T</math>
580
|-
581
| style="text-align: center;" | <math> {}[v_1,v_2,v_6] \longleftrightarrow  (0,1,0,0,0,0)^T</math>
582
|-
583
| style="text-align: center;" | <math> {}[v_2,v_3,v_6] \longleftrightarrow  (0,0,1,0,0,0)^T</math>
584
|-
585
| style="text-align: center;" | <math> {}[v_3,v_4,v_6] \longleftrightarrow  (0,0,0,1,0,0)^T</math>
586
|-
587
| style="text-align: center;" | <math> {}[v_4,v_5,v_6] \longleftrightarrow  (0,0,0,0,1,0)^T</math>
588
|-
589
| style="text-align: center;" | <math> {}[v_5,v_0,v_6] \longleftrightarrow  (0,0,0,0,0,1)^T</math>
590
|}
591
|}
592
593
Similarly,  coordinate vectors are associated to the edges
594
595
{| class="formulaSCP" style="width: 100%; text-align: left;" 
596
|-
597
| 
598
{| style="text-align: left; margin:auto;width: 100%;" 
599
|-
600
| style="text-align: center;" | <math>{} [v_0,v_6] \longleftrightarrow  (1,0,0,0,0,0,0,0,0,0,0,0)^T</math>
601
|-
602
| style="text-align: center;" | <math> {} [v_1,v_6] \longleftrightarrow  (0,1,0,0,0,0,0,0,0,0,0,0)^T</math>
603
|-
604
| style="text-align: center;" | <math> {} [v_2,v_6] \longleftrightarrow  (0,0,1,0,0,0,0,0,0,0,0,0)^T</math>
605
|-
606
| style="text-align: center;" | <math> {} [v_3,v_6] \longleftrightarrow  (0,0,0,1,0,0,0,0,0,0,0,0)^T</math>
607
|-
608
| style="text-align: center;" | <math> {} [v_4,v_6] \longleftrightarrow  (0,0,0,0,1,0,0,0,0,0,0,0)^T</math>
609
|-
610
| style="text-align: center;" | <math> {} [v_5,v_6] \longleftrightarrow  (0,0,0,0,0,1,0,0,0,0,0,0)^T</math>
611
|-
612
| style="text-align: center;" | <math> {} [v_0,v_1] \longleftrightarrow  (0,0,0,0,0,0,1,0,0,0,0,0)^T</math>
613
|-
614
| style="text-align: center;" | <math> {} [v_1,v_2] \longleftrightarrow  (0,0,0,0,0,0,0,1,0,0,0,0)^T</math>
615
|-
616
| style="text-align: center;" | <math> {} [v_2,v_3] \longleftrightarrow  (0,0,0,0,0,0,0,0,1,0,0,0)^T</math>
617
|-
618
| style="text-align: center;" | <math> {} [v_3,v_4] \longleftrightarrow  (0,0,0,0,0,0,0,0,0,1,0,0)^T</math>
619
|-
620
| style="text-align: center;" | <math> {} [v_4,v_5] \longleftrightarrow  (0,0,0,0,0,0,0,0,0,0,1,0)^T</math>
621
|-
622
| style="text-align: center;" | <math> {} [v_5,v_0] \longleftrightarrow  (0,0,0,0,0,0,0,0,0,0,0,1)^T</math>
623
|}
624
|}
625
626
Finally, we do the same with the vertices
627
628
{| class="formulaSCP" style="width: 100%; text-align: left;" 
629
|-
630
| 
631
{| style="text-align: left; margin:auto;width: 100%;" 
632
|-
633
| style="text-align: center;" | <math>{} [v_0] \longleftrightarrow  (1,0,0,0,0,0,0)^T</math>
634
|-
635
| style="text-align: center;" | <math> {} [v_1] \longleftrightarrow  (0,1,0,0,0,0,0)^T</math>
636
|-
637
| style="text-align: center;" | <math> {} [v_2] \longleftrightarrow  (0,0,1,0,0,0,0)^T</math>
638
|-
639
| style="text-align: center;" | <math> {} [v_3] \longleftrightarrow  (0,0,0,1,0,0,0)^T</math>
640
|-
641
| style="text-align: center;" | <math> {} [v_4] \longleftrightarrow  (0,0,0,0,1,0,0)^T</math>
642
|-
643
| style="text-align: center;" | <math> {} [v_5] \longleftrightarrow  (0,0,0,0,0,1,0)^T</math>
644
|-
645
| style="text-align: center;" | <math> {} [v_6] \longleftrightarrow  (0,0,0,0,0,0,1)^T</math>
646
|}
647
|}
648
649
Now, if we take the boundary of each face, we have
650
651
{| class="formulaSCP" style="width: 100%; text-align: left;" 
652
|-
653
| 
654
{| style="text-align: left; margin:auto;width: 100%;" 
655
|-
656
| style="text-align: center;" | <math>\partial [v_0,v_1,v_6] = [v_1,v_6] - [v_0,v_6]+[v_0,v_1]</math>
657
|-
658
| style="text-align: center;" | <math> \partial [v_1,v_2,v_6] = [v_2,v_6] - [v_1,v_6]+[v_1,v_2]</math>
659
|-
660
| style="text-align: center;" | <math> \partial [v_2,v_3,v_6] = [v_3,v_6] - [v_2,v_6]+[v_2,v_3]</math>
661
|-
662
| style="text-align: center;" | <math> \partial [v_3,v_4,v_6] = [v_4,v_6] - [v_3,v_6]+[v_3,v_4]</math>
663
|-
664
| style="text-align: center;" | <math> \partial [v_4,v_5,v_6] = [v_5,v_6] - [v_4,v_6]+[v_4,v_5]</math>
665
|-
666
| style="text-align: center;" | <math> \partial [v_5,v_0,v_6] = [v_0,v_6] - [v_5,v_6]+[v_5,v_0] </math>
667
|}
668
|}
669
670
which, under the previous assignments of coordinate vectors, corresponds to the linear transformation given by the following matrix which will multiply the column coordinate vectors on the left
671
672
{| class="formulaSCP" style="width: 100%; text-align: left;" 
673
|-
674
| 
675
{| style="text-align: left; margin:auto;width: 100%;" 
676
|-
677
| style="text-align: center;" | <math>\partial _{2,1}=\left(\begin{array}{rrrrrr} -1 & 0 & 0 & 0 & 0 & 1\\ 1  & -1 & 0 & 0 &0 & 0\\  0  & 1 & -1 & 0 & 0 & 0\\  0  & 0 & 1 & -1 & 0 & 0\\  0  &0  &0  & 1 & -1 & 0\\  0  & 0 & 0 & 0 & 1 & -1\\ 1  & 0 & 0 & 0 & 0 & 0\\ 0   & 1 & 0 & 0 & 0 & 0\\  0  & 0 & 1 &0  &0  &0 \\   0 & 0 & 0 & 1 & 0 & 0\\  0  & 0 & 0 & 0 & 1 & 0\\  0  & 0 & 0 & 0 & 0 & 1         \end{array} \right)</math>
678
|}
679
|}
680
681
where the subindices in <math display="inline">\partial _{2,1}</math> indicate that we are taking the boundary of 2-dimensional elements and obtaining 1-dimensional ones. Similarly, taking the boundaries of all the edges gives
682
683
{| class="formulaSCP" style="width: 100%; text-align: left;" 
684
|-
685
| 
686
{| style="text-align: left; margin:auto;width: 100%;" 
687
|-
688
| style="text-align: center;" | <math>\partial [v_0,v_6] = [v_6] - [v_0]</math>
689
|-
690
| style="text-align: center;" | <math> \partial [v_1,v_6] = [v_6] - [v_1]</math>
691
|-
692
| style="text-align: center;" | <math> \partial [v_2,v_6] = [v_6] - [v_2]</math>
693
|-
694
| style="text-align: center;" | <math> \partial [v_3,v_6] = [v_6] - [v_3]</math>
695
|-
696
| style="text-align: center;" | <math> \partial [v_4,v_6] = [v_6] - [v_4]</math>
697
|-
698
| style="text-align: center;" | <math> \partial [v_5,v_6] = [v_6] - [v_5]</math>
699
|-
700
| style="text-align: center;" | <math> \partial [v_0,v_1] = [v_1] - [v_0]</math>
701
|-
702
| style="text-align: center;" | <math> \partial [v_1,v_2] = [v_2] - [v_1]</math>
703
|-
704
| style="text-align: center;" | <math> \partial [v_2,v_3] = [v_3] - [v_2]</math>
705
|-
706
| style="text-align: center;" | <math> \partial [v_3,v_4] = [v_4] - [v_3]</math>
707
|-
708
| style="text-align: center;" | <math> \partial [v_4,v_5] = [v_5] - [v_4]</math>
709
|-
710
| style="text-align: center;" | <math> \partial [v_5,v_0] = [v_0] - [v_5] </math>
711
|}
712
|}
713
714
which, under the previous assignments of coordinate vectors, corresponds to the linear transformation given by the following matrix
715
716
{| class="formulaSCP" style="width: 100%; text-align: left;" 
717
|-
718
| 
719
{| style="text-align: left; margin:auto;width: 100%;" 
720
|-
721
| style="text-align: center;" | <math>\partial _{1,0}=\left( \begin{array}{rrrrrr|rrrrrr} -1 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & 0 & 1\\ 0 & -1 & 0 & 0 & 0 & 0 & 1 & -1 & 0 & 0 & 0 & 0\\ 0 & 0 & -1 & 0 & 0 & 0 & 0 & 1 & -1 & 0 & 0 & 0\\ 0 & 0 & 0 & -1 & 0 & 0 & 0 & 0 & 1 & -1 & 0 & 0\\ 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & 0 & 1 & -1 & 0\\ 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & 0 & 1 & -1\\ 1 & 1 & 1 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \end{array} \right)</math>
722
|}
723
|}
724
725
'''Remark'''. Note how these matrices encode different levels of connectivity with orientations, such as who are the edges of which oriented triangle, or which are the  end points of a given oriented edge.
726
727
Due to the duality between <math display="inline">\partial </math> and <math display="inline">\nabla ^\wedge </math>, we can define the discretiztion of <math display="inline">\nabla ^\wedge </math>  by
728
729
{| class="formulaSCP" style="width: 100%; text-align: left;" 
730
|-
731
| 
732
{| style="text-align: left; margin:auto;width: 100%;" 
733
|-
734
| style="text-align: center;" | <math>\nabla ^\wedge := (\partial )^T</math>
735
|}
736
|}
737
738
For instance, we see that the operator <math display="inline">\nabla ^\wedge _{0,1}</math> reads as follows
739
740
{| class="formulaSCP" style="width: 100%; text-align: left;" 
741
|-
742
| 
743
{| style="text-align: left; margin:auto;width: 100%;" 
744
|-
745
| style="text-align: center;" | <math>\nabla ^\wedge _{0,1}=\left(\begin{array}{rrrrrrr} -1 & 0 & 0 & 0 & 0 & 0 & 1\\ 0 & -1 & 0 & 0 & 0 & 0 & 1\\ 0 & 0 & -1 & 0 & 0 & 0 & 1\\ 0 & 0 & 0 & -1 & 0 & 0 & 1\\ 0 & 0 & 0 & 0 & -1 & 0 & 1\\ 0 & 0 & 0 & 0 & 0 & -1 & 1\\ -1 & 1 & 0 & 0 & 0 & 0 & 0\\ 0 & -1 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & -1 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & -1 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & -1 & 1& 0 \\ 1 & 0 & 0 & 0 & 0 & -1& 0            \end{array} \right)</math>
746
|}
747
|}
748
749
===3.2 Dual mesh===
750
751
In order to discretize the Hodge star operator, we must first introduce  the notion of the dual mesh of a triangular mesh.
752
753
Consider the triangular mesh <math display="inline">K</math> in Figure [[#3.2 Dual mesh|3.6]](a). The construction of the dual mesh <math display="inline">K^*</math> is carried out as follows:
754
755
* The vertices of the dual mesh <math display="inline">K^*</math> are the circumcenters of the faces/triangles of the original mesh  (the blue dots in Figures   [[#3.2 Dual mesh|3.6]](b) and [[#3.2 Dual mesh|3.6]](c)).  For instance, the dual of the face <math display="inline">[v_0,v_1,v_6]</math> will be denoted by <math display="inline">[v_0,v_1,v_6]^*</math>.
756
* The edges of <math display="inline">K^*</math> are the straight line segments joining the circumcenters of two adjacent triangles (those which share an edge). Note that the resulting line segments   are orthogonal to one of the original edges (the blue straight line segments in   in Figures   [[#3.2 Dual mesh|3.6]](b) and [[#3.2 Dual mesh|3.6]](c)). For instance, the dual of the edge <math display="inline">[v_0,v_6]</math> will be denoted by <math display="inline">[v_0,v_6]^*</math>.
757
* The faces or cells of <math display="inline">K^*</math> are the areas enclosed by the new polygons determined by the new edges. For instance, the dual of the vertex <math display="inline">[v_6]</math> is the inner blue hexagon in  in Figures   [[#3.2 Dual mesh|3.6]](b) and [[#3.2 Dual mesh|3.6]](c)  and will be denoted by <math display="inline">[v_6]^*</math>.
758
759
760
761
{| class="floating_imageSCP" style="text-align: center;"
762
|-
763
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
764
<span style="text-align: center; font-size: 75%;">'' 
765
(a)  [[Image:Review_239195285364_6808_text-boundary04-eps-converted-to.png|150px|]]<math>\qquad</math>
766
(b)  [[Image:Review_239195285364_2587_text-Dual00-eps-converted-to.png|150px|]]<math>\qquad</math>
767
(c)  [[File:Review_239195285364_6445_Dual01-eps-converted-to.jpg|120px|]]
768
''</span></div>
769
|}
770
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
771
<span style="text-align: center; font-size: 75%;">'''Figure 3.6:''' Dual mesh construction: (a) triangular mesh <math>K</math>; (b) dual mesh <math>K^*</math> superimposed on the mesh <math>K</math>;      (c) dual mesh <math>K^*</math>.</span></div>
772
773
774
775
The orientation of the dual edges is given by the following recipe. If we have two adjacent triangles oriented as in  Figure [[#3.2 Dual mesh|3.7]](a), the dual edge crossing the edge of adjacency is oriented as in Figure [[#3.2 Dual mesh|3.7]](b).
776
777
778
{| class="floating_imageSCP" style="text-align: center; "
779
|-
780
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
781
<span style="text-align: center; font-size: 75%;">''(a)  [[Image:Review_239195285364_3618_text-TwoTriangles01-eps-converted-to.png|120px|]]<math>\qquad</math>(b) [[Image:Review_239195285364_1851_text-TwoTriangles02-eps-converted-to.png|120px|]]''</span></div>
782
|}
783
784
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
785
<span style="text-align: center; font-size: 75%;">'''Figure 3.7:''' (a) Two adjacent oriented triangles. (b) Compatibly oriented dual edge..</span></div>
786
787
===3.3 Boundary operator on the dual mesh===
788
789
Consider the dual mesh in Figure [[#3.3 Boundary operator on the dual mesh|3.8]] with the given labels and orientations.
790
791
{| class="floating_imageSCP" style="text-align: center; "
792
|-
793
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
794
<span style="text-align: center; font-size: 75%;">''[[Image:Review_239195285364_5980_text-Dual01a-eps-converted-to.png|160px|]]<math>\qquad</math>[[Image:Review_239195285364_7382_text-Dual01b-eps-converted-to.png|160px|n]]''</span></div>
795
|}
796
797
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
798
<span style="text-align: center; font-size: 75%;">'''Figure 3.8:''' Oriented dual mesh.</span></div>
799
800
801
The boundary operator is applied in a similar fashion as it was applied to triangles. In this case we have
802
803
{| class="formulaSCP" style="width: 100%; text-align: left;" 
804
|-
805
| 
806
{| style="text-align: left; margin:auto;width: 100%;" 
807
|-
808
| style="text-align: center;" | <math>\partial ^{dual}_{2,1} F_1 =  E_{1}  + E_{7}  -  E_{12}</math>
809
|-
810
| style="text-align: center;" | <math>  \partial ^{dual}_{2,1} F_2 =  E_{2}  - E_{7}  +  E_{8}</math>
811
|-
812
| style="text-align: center;" | <math>  \partial ^{dual}_{2,1} F_3 =  E_{3}  - E_{8}  +  E_{9} </math>
813
|-
814
| style="text-align: center;" | <math>  \partial ^{dual}_{2,1} F_4 =  E_{4}  - E_{9}  +  E_{10}</math>
815
|-
816
| style="text-align: center;" | <math>  \partial ^{dual}_{2,1} F_5 =  E_{5}  - E_{10} +  E_{11}</math>
817
|-
818
| style="text-align: center;" | <math>  \partial ^{dual}_{2,1} F_6 =  E_{6}  - E_{11} +  E_{12}</math>
819
|-
820
| style="text-align: center;" | <math>  \partial ^{dual}_{2,1} F_7 =  -E_{1}- E_{2}- E_{3}- E_{4}- E_{5}- E_{6} </math>
821
|}
822
|}
823
824
If we assign coordinate vectors to the dual faces and the dual edges as before
825
826
{| class="formulaSCP" style="width: 100%; text-align: left;" 
827
|-
828
| 
829
{| style="text-align: left; margin:auto;width: 100%;" 
830
|-
831
| style="text-align: center;" | <math>F_1 \longleftrightarrow  (1,0,0,0,0,0,0)^T</math>
832
|-
833
| style="text-align: center;" | <math>  F_2 \longleftrightarrow  (0,1,0,0,0,0,0)^T</math>
834
|-
835
| style="text-align: center;" | <math>  F_3 \longleftrightarrow  (0,0,1,0,0,0,0)^T</math>
836
|-
837
| style="text-align: center;" | <math>  F_4 \longleftrightarrow  (0,0,0,1,0,0,0)^T</math>
838
|-
839
| style="text-align: center;" | <math>  F_5 \longleftrightarrow  (0,0,0,0,1,0,0)^T</math>
840
|-
841
| style="text-align: center;" | <math>  F_6 \longleftrightarrow  (0,0,0,0,0,1,0)^T</math>
842
|-
843
| style="text-align: center;" | <math>  F_7 \longleftrightarrow  (0,0,0,0,0,0,1)^T</math>
844
|}
845
|}
846
847
and
848
849
{| class="formulaSCP" style="width: 100%; text-align: left;" 
850
|-
851
| 
852
{| style="text-align: left; margin:auto;width: 100%;" 
853
|-
854
| style="text-align: center;" | <math>E_1 \longleftrightarrow  (1,0,0,0,0,0,0,0,0,0,0,0)^T</math>
855
|-
856
| style="text-align: center;" | <math>  E_2 \longleftrightarrow  (0,1,0,0,0,0,0,0,0,0,0,0)^T</math>
857
|-
858
| style="text-align: center;" | <math>  E_3 \longleftrightarrow  (0,0,1,0,0,0,0,0,0,0,0,0)^T</math>
859
|-
860
| style="text-align: center;" | <math>  E_4 \longleftrightarrow  (0,0,0,1,0,0,0,0,0,0,0,0)^T</math>
861
|-
862
| style="text-align: center;" | <math>  E_5 \longleftrightarrow  (0,0,0,0,1,0,0,0,0,0,0,0)^T</math>
863
|-
864
| style="text-align: center;" | <math>  E_6 \longleftrightarrow  (0,0,0,0,0,1,0,0,0,0,0,0)^T</math>
865
|-
866
| style="text-align: center;" | <math>  E_7 \longleftrightarrow  (0,0,0,0,0,0,1,0,0,0,0,0)^T</math>
867
|-
868
| style="text-align: center;" | <math>  E_8 \longleftrightarrow  (0,0,0,0,0,0,0,1,0,0,0,0)^T</math>
869
|-
870
| style="text-align: center;" | <math>  E_9 \longleftrightarrow  (0,0,0,0,0,0,0,0,1,0,0,0)^T</math>
871
|-
872
| style="text-align: center;" | <math>  E_{10} \longleftrightarrow  (0,0,0,0,0,0,0,0,0,1,0,0)^T</math>
873
|-
874
| style="text-align: center;" | <math>  E_{11} \longleftrightarrow  (0,0,0,0,0,0,0,0,0,0,1,0)^T</math>
875
|-
876
| style="text-align: center;" | <math>  E_{12} \longleftrightarrow  (0,0,0,0,0,0,0,0,0,0,0,1)^T</math>
877
|}
878
|}
879
880
we have the associated matrix
881
882
{| class="formulaSCP" style="width: 100%; text-align: left;" 
883
|-
884
| 
885
{| style="text-align: left; margin:auto;width: 100%;" 
886
|-
887
| style="text-align: center;" | <math>\partial _{2,1}^{dual}=\left(\begin{array}{rrrrrrr} 1 & 0 & 0 & 0 & 0 & 0 & -1\\ 0 & 1 & 0 & 0 & 0 & 0 & -1\\ 0 & 0 & 1 & 0 & 0 & 0 & -1\\ 0 & 0 & 0 & 1 & 0 & 0 & -1\\ 0 & 0 & 0 & 0 & 1 & 0 & -1\\ 0 & 0 & 0 & 0 & 0 & 1 & -1\\ 1 & -1 & 0 & 0 & 0 & 0 & 0\\ 0 & 1 & -1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & -1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & -1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & -1& 0 \\ -1 & 0 & 0 & 0 & 0 & 1& 0            \end{array} \right)</math>
888
|}
889
|}
890
891
Notice that
892
893
{| class="formulaSCP" style="width: 100%; text-align: left;" 
894
|-
895
| 
896
{| style="text-align: left; margin:auto;width: 100%;" 
897
|-
898
| style="text-align: center;" | <math>\partial _{2,1}^{dual}=-\left(\partial _{1,0}\right)^T</math>
899
|}
900
|}
901
902
which arises from the duality between the two meshes <span id='citeF-6'></span>[[#cite-6|[6]]]. In general, the discrete differential to be applied is
903
904
{| class="formulaSCP" style="width: 100%; text-align: left;" 
905
|-
906
| 
907
{| style="text-align: left; margin:auto;width: 100%;" 
908
|-
909
| style="text-align: center;" | <math>\nabla _{1,2}^{\wedge ,dual}=\left(\partial _{2,1}^{dual}\right)^T=-\partial _{1,0}=-\left(\nabla ^\wedge _{0,1}\right)^T</math>
910
|}
911
|}
912
913
===3.4 Discrete Hodge star===
914
915
The discretization of the Hodge star <math display="inline">\star </math> uses the    geometrical ideas described in Section [[#lb-2.2|2.2]] and the dual mesh. More precisely, the 2D Hodge star operator rotates a vector <math display="inline">90^\circ </math> counterclockwise.  For the sake of clarity, let us focus on the edge <math display="inline">[v_0,v_6]</math>, its mesh dual <math display="inline">[v_0,v_6]^*</math> and its Hodge star image <math display="inline">\star [v_0,v_6]</math>. They are represented in Figure [[#3.4 Discrete Hodge star|3.9]].
916
917
{| class="floating_imageSCP" style="text-align: center;"
918
|-
919
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
920
<span style="text-align: center; font-size: 75%;">''[[Image:Review_239195285364_9611_text-Dual00a-eps-converted-to.png|160px]]<math>\qquad</math>[[Image:Review_239195285364_4015_text-Dual00aa-eps-converted-to.png|160px]]''</span></div>
921
|}
922
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
923
<span style="text-align: center; font-size: 75%;">'''Figure 3.9'''</span></div>
924
925
Since
926
927
{| class="formulaSCP" style="width: 100%; text-align: left;" 
928
|-
929
| 
930
{| style="text-align: left; margin:auto;width: 100%;" 
931
|-
932
| style="text-align: center;" | <math>{length}(\star [v_0,v_6])={length}([v_0,v_6])</math>
933
|}
934
|}
935
936
we see that the relationship between the dual edge <math display="inline">[v_0,v_6]^*</math> and the geometric <math display="inline">\star [v_0,v_6]</math> is the following
937
938
{| class="formulaSCP" style="width: 100%; text-align: left;" 
939
|-
940
| 
941
{| style="text-align: left; margin:auto;width: 100%;" 
942
|-
943
| style="text-align: center;" | <math>{1\over {length}([v_0,v_6]^*)}[v_0,v_6]^*= {1\over {length}([v_0,v_6])}\star [v_0,v_6]</math>
944
|}
945
|}
946
947
As can be seen, when applying the geometric Hodge star to the edges of the mesh,  we do not end up in the dual mesh but in multiples of the elements of the dual mesh. Thus, <math display="inline">\star [v_0,v_6]</math>  must be scaled to match <math display="inline">[v_0,v_6]^*</math>. If we do this to all the Hodge star images, we get the discrete Hodge star matrix
948
949
{| class="formulaSCP" style="width: 100%; text-align: left;" 
950
|-
951
| 
952
{| style="text-align: left; margin:auto;width: 100%;" 
953
|-
954
| style="text-align: center;" | <math>M_{1,1}=\left(\begin{array}{ccccc} {{length}([v_0,v_6]^*)\over {length}([v_0,v_6])} & 0 & 0 & \dots & 0 \\ 0 & {{length}([v_1,v_6]^*)\over {length}([v_1,v_6])} & 0 & \dots & 0 \\ 0 & 0 & {{length}([v_2,v_6]^*)\over {length}([v_2,v_6])} & \dots & 0\\ \vdots & \vdots & \vdots &    \ddots &\vdots \\ 0 & 0 & 0 & \dots &  {{length}([v_5,v_0]^*)\over {length}([v_5,v_0])}           \end{array} \right)</math>
955
|}
956
|}
957
958
where the subindices in <math display="inline">M_{1,1}</math> indicate that we are sending 1-dimensional elements of the original mesh to 1-dimensional elements of the dual mesh.
959
960
Somewhat less intuitive is the meaning of the geometric Hodge star operator on nodes of the original mesh. As we saw in Subsection [[#lb-2.2|2.2]],
961
962
{| class="formulaSCP" style="width: 100%; text-align: left;" 
963
|-
964
| 
965
{| style="text-align: left; margin:auto;width: 100%;" 
966
|-
967
| style="text-align: center;" | <math>\star 1=e_1\wedge e_2</math>
968
|}
969
|}
970
971
which geometrically means that the Hodge star <math display="inline">\star [v_6]</math> must be a polygon with area equal to 1 (classically, it is a parallelogram, but can also be a hexahedron of area 1 as in this example).  
972
973
{| class="floating_imageSCP" style="text-align: center; "
974
|-
975
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
976
<span style="text-align: center; font-size: 75%;">''[[Image:Review_239195285364_3989_text-HodgeStarPoint01-eps-converted-to.png|160px|]]''</span></div>
977
|}
978
979
However, in the dual mesh we have the polygon <math display="inline">[v_6]^*</math>       
980
981
{| class="floating_imageSCP" style="text-align: center; "
982
|-
983
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
984
<span style="text-align: center; font-size: 75%;">''[[Image:Review_239195285364_2797_text-HodgeStarPoint00-eps-converted-to.png|160px|]]''</span></div>
985
|}
986
987
so that we need to resize <math display="inline">\star [v_6]</math> as follows
988
989
{| class="formulaSCP" style="width: 100%; text-align: left;" 
990
|-
991
| 
992
{| style="text-align: left; margin:auto;width: 100%;" 
993
|-
994
| style="text-align: center;" | <math>[v_6]^*={Area}([v_6]^*)\, \star [v_6]</math>
995
|}
996
|}
997
998
If we do that for all the vertices, we obtain the discrete Hodge star matrix
999
1000
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1001
|-
1002
| 
1003
{| style="text-align: left; margin:auto;width: 100%;" 
1004
|-
1005
| style="text-align: center;" | <math> M_{0,2}=\left( \begin{array}{cccc} {Area}([v_0]^*) & 0 & \dots & 0\\ 0 & {Area}([v_1]^*) & \dots & 0\\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & {Area}([v_6]^*) \end{array} \right)</math>
1006
|}
1007
|}
1008
1009
The inverse matrix will deal with the case when we take the Hodge star of the 2D polygons in the dual mesh to obtain points (with weight 1) in the original mesh.
1010
1011
In summary, the various matrices <math display="inline">M</math> representing the discrete Hodge star operator send elements of the original mesh to elements of the dual mesh.
1012
1013
===3.5 DEC applied to 2D Poisson's equation===
1014
1015
Consider the 2D Poisson's equation
1016
1017
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1018
|-
1019
| 
1020
{| style="text-align: left; margin:auto;width: 100%;" 
1021
|-
1022
| style="text-align: center;" | <math>\kappa \Delta f= q</math>
1023
|}
1024
|}
1025
1026
As we have seen, this can be rewritten as
1027
1028
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1029
|-
1030
| 
1031
{| style="text-align: left; margin:auto;width: 100%;" 
1032
|-
1033
| style="text-align: center;" | <math>\kappa \,\,\star \nabla ^\wedge \star \nabla (f)     = q </math>
1034
|}
1035
|}
1036
1037
Suppose that we wish to solve the equation using the mesh <math display="inline">K</math>. Let us denote by <math display="inline">[f]</math> and <math display="inline">[q]</math> the column vector discretizations of the functions <math display="inline">K</math> and <math display="inline">q</math> at the nodes. According to DEC, the first derivation <math display="inline">\nabla(f)</math> is achieved left-multiplying <math display="inline">[f]</math> by the matrix <math display="inline">\nabla ^\wedge _{0,1}</math> of Subsection 3.1. 
1038
The resulting column vector is a collection of (unknown) directional derivative values of <math display="inline">f</math> assigned to the ordered set of edges. At
1039
this point we need to transfer such a set of values to the ordered set of dual edges of the dual mesh with the appropriate geometrical scaling, which is achieved multiplying on the left by the matrix <math display="inline">M_{1,1}</math> of Subsection 3.4. Now, the discrete version of the second derivation <math display="inline">\nabla^\wedge</math> is given by left multiplication with the matrix <math display="inline">\nabla^\wedge_{1,2}</math> of Subsection 3.3, which gives a column vector of (unknown) values assigned to the ordered set of 2-dimensional cells of the dual mesh. Finally, the task of mapping such values (with the appropriate scaling) to the ordered set of vertices of the original mesh is achieved by left multiplication with the matrix <math display="inline">M_{2,0}</math> of Subsection 3.4. Thus, the Poisson equation is discretized as the matrix equation
1040
1041
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1042
|-
1043
| 
1044
{| style="text-align: left; margin:auto;width: 100%;" 
1045
|-
1046
| style="text-align: center;" | <math>\kappa \,\,M_{2,0}\,\, \nabla _{1,2}^{\wedge ,dual} \,\, M_{1,1}\,\, \nabla ^\wedge _{0,1}\,\, [f]=[q]</math>
1047
|}
1048
|}
1049
1050
Later on, it will be convenient to work with the equivalent system
1051
1052
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1053
|-
1054
| 
1055
{| style="text-align: left; margin:auto;width: 100%;" 
1056
|-
1057
| style="text-align: center;" | <math>\kappa \,\,(\nabla ^\wedge _{0,1})^T\,\,M_{1,1}\,\,\nabla ^\wedge _{0,1}[f]= M_{0,2}[q]</math>
1058
|}
1059
|}
1060
1061
Note that the DEC-discretization of <math display="inline">\nabla(f)</math> is not a (discretized) flow vector. It is a collection of values
1062
of directional derivatives of <math display="inline">f</math>, one such value for each oriented edge.
1063
1064
==4 DEC for general triangulations==
1065
1066
Since the boundary operator is really concerned with the connectivity of the mesh and does not change under deformation of the mesh,  the change in the setup of DEC for a deformed mesh must be contained in the discrete Hodge star matrices. Since  such matrices are computed in terms of lengths and areas of oriented elements of the mesh, we will now examine how those ingredients  transform under deformation, a problem that was first considered in <span id='citeF-10'></span>[[#cite-10|[10]]].
1067
1068
===4.1 Dual mesh of an arbitrary triangle===
1069
1070
In order to explain how to implement DEC for general triangulations, let us consider first a mesh consisting of a single well-centered triangle,  as well as its dual mesh (see Figure [[#4.1 Dual mesh of an arbitrary triangle|4.1]]).
1071
1072
{| class="floating_imageSCP" style="text-align: center; "
1073
|-
1074
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
1075
<span style="text-align: center; font-size: 75%;">''[[Image:Review_239195285364_2549_NewCircumscribedTriangle00.png|160px|]]<math>\qquad</math>[[Image:Review_239195285364_2438_NewCircumscribedTriangle00b.png|180px|]]''</span></div>
1076
|}
1077
1078
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
1079
<span style="text-align: center; font-size: 75%;">'''Figure 4.1:''' Well-centered triangle and its dual mesh.</span></div>
1080
1081
1082
The dual cells are given as follows:
1083
1084
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1085
|-
1086
| 
1087
{| style="text-align: left; margin:auto;width: 100%;" 
1088
|-
1089
| style="text-align: center;" | <math>{} [v_1,v_2,v_3]^*= [c]</math>
1090
|-
1091
| style="text-align: center;" | <math> {} [v_1,v_2]^*= [p_1,c]</math>
1092
|-
1093
| style="text-align: center;" | <math> {} [v_2,v_3]^*= [p_2,c]</math>
1094
|-
1095
| style="text-align: center;" | <math> {} [v_3,v_1]^*= [p_3,c]</math>
1096
|-
1097
| style="text-align: center;" | <math> {} [v_1]^* = [v_1,p_1,c,p_3] </math>
1098
|-
1099
| style="text-align: center;" | <math> {} [v_2]^* = [v_2,p_2,c,p_1] </math>
1100
|-
1101
| style="text-align: center;" | <math> {} [v_3]^* = [v_3,p_3,c,p_2]  </math>
1102
|}
1103
|}
1104
1105
Now consider the cell <math display="inline">[v_3]^* = [v_3,p_3,c,p_2]</math> subdivided as in  Figure [[#4.1 Dual mesh of an arbitrary triangle|4.2]].
1106
1107
{| class="floating_imageSCP" style="text-align: center; "
1108
|-
1109
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
1110
<span style="text-align: center; font-size: 75%;">''[[Image:Review_239195285364_1856_NewCircumscribedTriangle00a.png|160px|]]''</span></div>
1111
|}
1112
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
1113
<span style="text-align: center; font-size: 75%;">'''Figure 4.2:''' The subdivision of a 2-dimensional dual cell of a well-centered triangle.</span></div>
1114
1115
1116
If we deform continuously the triangle <math display="inline">[v_1,v_2,v_3]</math> to become an obtuse triangle as in  Figure [[#4.1 Dual mesh of an arbitrary triangle|4.3]], we see in Figures [[#4.1 Dual mesh of an arbitrary triangle|4.1]](a) and [[#4.1 Dual mesh of an arbitrary triangle|4.1]](b) that the area of the blue subtriangle <math display="inline">[v_3,c,p_3]</math> decreases to 0 and in Figure [[#4.1 Dual mesh of an arbitrary triangle|4.1]](c) that it is completely outside of the triangle and, therefore, must be assigned a negative sign. The same can be said about the 1-dimensional cell <math display="inline">[p_3,c]</math>, which originally is completely contained in the triangle <math display="inline">[v_1,v_2,v_3]</math>, its size reduces to zero as the triangle is deformed (Figures [[#4.1 Dual mesh of an arbitrary triangle|4.1]](a) and [[#4.1 Dual mesh of an arbitrary triangle|4.1]](b)),  and eventually it is completely outside the triangle <math display="inline">[v_1,v_2,v_3]</math>  (Figure [[#4.1 Dual mesh of an arbitrary triangle|4.1]](c))  and a negative sign must be assigned to it. On the other hand, part of the red subtriangle <math display="inline">[v_3,c,p_2]</math> still intersects the interior of the triangle <math display="inline">[v_1,v_2,v_3]</math> and,  therefore, no assignment of sign is needed. Similarly for the segment <math display="inline">[p_2,c]</math>. In terms of numerical simulations, an implementation in terms of determinants contains intrinscally the aforementioned change of signs.
1117
1118
{| class="floating_imageSCP" style="text-align: center; "
1119
|-
1120
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
1121
<span style="text-align: center; font-size: 75%;">''[[Image:Review_239195285364_9664_NewCircumscribedTriangle01.png|150px|]](a)<math>\;</math>[[Image:Review_239195285364_6104_NewCircumscribedTriangle02.png|150px|]](b)<math>\;</math>[[Image:Review_239195285364_5608_NewCircumscribedTriangle03.png|150px|]](c)''</span></div>
1122
|}
1123
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
1124
<span style="text-align: center; font-size: 75%;">'''Figure 4.3:''' The subdivision of a 2-dimensional dual cell of a well-centered triangle.</span></div>
1125
1126
===4.2 Dual mesh of a general triangulation===
1127
1128
Now consider the well-centered mesh and its dual in Figure [[#4.2 Dual mesh of a general triangulation|4.4]].
1129
1130
1131
{| class="floating_imageSCP" style="text-align: center; margin: 1em auto; width: 100%;max-width: 100%;"
1132
|-
1133
|[[File:Review_239195285364_3465_NewHexagonalMesh00.png|160px|]]
1134
|[[File:Review_239195285364_6131_NewHexagonalMesh01.png|160px|]]
1135
|- style="text-align: center; font-size: 75%;"
1136
| colspan="2" | '''Figure 4.4:''' (a) Well-centered triangular mesh of hexagon. (b) Dual mesh.
1137
|}
1138
1139
Observe the deformation of the blue-colored dual cell <math display="inline">[v_6]^*</math>  in Figure [[#4.2 Dual mesh of a general triangulation|4.5]](a) as the vertex <math display="inline">v_0</math> is moved to make the triangle <math display="inline">[v_0,v_6,v_5]</math> non-well-centered in Figure [[#4.2 Dual mesh of a general triangulation|4.5]](b) and the vertex <math display="inline">v_4</math> is moved to make the triangle <math display="inline">[v_4,v_5,v_6]</math>  non-well-centered in Figure [[#4.2 Dual mesh of a general triangulation|4.5]](c).
1140
1141
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1142
|-
1143
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
1144
<span style="text-align: center; font-size: 75%;">''(a)[[File:Review_239195285364_7165_NewHexagonalMesh02.png|150px|]]<math>\;</math>(b)[[File:Review_239195285364_6070_NewHexagonalMesh03.png|150px|]]<math>\;</math>(c)[[File:Review_239195285364_6698_NewHexagonalMesh04.png|150px|]]''</span></div>
1145
|}
1146
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
1147
<span style="text-align: center; font-size: 75%;">'''Figure 4.5:''' Deformation of a 2-dimensional dual cell as the hexagonal region is deformed.</span></div>
1148
 
1149
1150
We have colored in red the part of the dual cell <math display="inline">[v_6]^*</math> in Figure [[#4.2 Dual mesh of a general triangulation|4.2]](c) that must have the opposite orientation of the blue-colored part. Similar considerations apply to the edges.
1151
1152
===4.3. Brief digression in 2D.===
1153
In 2D, there is an intuitively clear way of comparing Discrete Exterior Calculus with the Finite Element Method (FEM) and the Finite Volume Method (FVM), both with linear interpolation functions, by considering the circumcenters of the triangles.
1154
*  Since the circumcenter of an obtuse triangle lies outside the triangle, the circumcenter cannot be used as control points of an irregular mesh for FVM.
1155
*  DEC is capable of using the circumcenter of irregular meshes. While the rigidity matrix generated by DEC works in a similar fashion to that of FEM, the discretized source term generally differs.
1156
*  When the elements of the mesh are all regular (the circumcenter of every triangle is in the triangle), DEC, FEM and FVM are very similar.
1157
1158
==5 Numerical examples==
1159
1160
===5.1 First example===
1161
1162
Let us solve the Poisson equation in a circle of radius one under the following conditions (see Figure [[#5.1 First example|5.1]]):
1163
1164
* heat difussion constant <math display="inline">k=1</math>;
1165
* source term <math display="inline">q= -1</math>;
1166
* Dirichlet boundary condition <math display="inline">u=10</math>.
1167
1168
The exact solution is
1169
1170
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1171
|-
1172
| 
1173
{| style="text-align: left; margin:auto;width: 100%;" 
1174
|-
1175
| style="text-align: center;" | <math> u(x,y)=\frac{1}{4}(1-x^2-y^2)+10 </math>
1176
|}
1177
|}
1178
1179
{| class="floating_imageSCP" style="text-align: center; "
1180
|-
1181
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
1182
<span style="text-align: center; font-size: 75%;">''[[Image:Review_239195285364_5542_Fig_6_circle_prb.png|160px|]]''</span></div>
1183
|}
1184
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
1185
<span style="text-align: center; font-size: 75%;">'''Figure 5.1:''' Disk of radius one with difussion constant <math>k=1</math>, subject to a heat source <math>q=-1</math>.</span></div>
1186
1187
1188
The meshes used in this example are shown in Figure [[#5.1 First example|5.2]] and vary from very coarse to very fine.
1189
1190
{| class="floating_imageSCP" style="text-align: center;"
1191
|-
1192
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
1193
<span style="text-align: center; font-size: 75%;">''(a)[[Image:Review_239195285364_2714_Fig_7_circle_msh1.png|150px|]]<math>\;</math>(b)[[Image:Review_239195285364_3494_Fig_7_circle_msh2.png|150px|]]<math>\;</math>(c)[[Image:Review_239195285364_1920_Fig_7_circle_msh3.png|150px|]]''</span></div>
1194
1195
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
1196
<span style="text-align: center; font-size: 75%;">''(d)[[Image:Review_239195285364_1099_Fig_7_circle_msh4.png|150px|]]<math>\;</math>(e)|[[Image:Review_239195285364_9301_Fig_7_circle_msh5.png|150px|]]''</span></div>
1197
|}
1198
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
1199
<span style="text-align: center; font-size: 75%;">'''Figure 5.2:''' Meshes for unit disk.</span></div>
1200
1201
1202
The numerical results for the maximum temperature value (<math display="inline">u(0,0)=10.25</math>) are exemplified in Table [[#5.1 First example|5.1]] where  a comparison with the Finite Element Method with linear interpolation functions (FEML) is also shown.  The FEML methodology that we used in the comparison can be consulted <span id='citeF-13'></span><span id='citeF-12'></span><span id='citeF-2'></span>[[#cite-13|[13,12,2]]]. For the sake of completeness in our comparison, here we compute the flux vectors in the same way as in FEML.
1203
1204
1205
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1206
|-
1207
| 
1208
{| style="text-align: left; margin:auto;width: 100%;" 
1209
|-
1210
| style="text-align: center;" | <math> \begin{array}{|c|c|c|c|c|c|c|} {Mesh}& \# {nodes}& \# {elements} & {Max. Temp. Value} & {Max. Flux Magnitude}  \\ &&&{DEC\quad FEML}&{DEC\quad FEML}  \\  {Figure \,\,\, (a)}& 9& 8& {10.250\quad 10.285}& {0.270\quad 0.307}\\ {Figure \,\,\,(b)}& 17& 20& {10.250\quad 10.237}& {0.388\quad 0.405}\\ {Figure \,\,\,(c)}& 41& 56& {10.250\quad 10.246}& {0.449\quad 0.453}\\ {Figure \,\,\,(d)}& 713& 1304& {10.250\quad 10.250}& {0.491\quad 0.492}\\ {Figure \,\,\,(e)}& 42298& 83346& {10.250\quad 10.250}& {0.496\quad 0.496}\\ \end{array}     </math>
1211
|}
1212
|- style="text-align: center; font-size: 75%;"
1213
| colspan="3" | '''Table 5.1:''' Maximum temperature and Flux magnitude values in the numerical simulation.
1214
|}
1215
1216
1217
The temperature distribution and Flux magnitude fields for the finest mesh are shown in Figure [[#5.1 First example|5.3]].
1218
1219
{| class="floating_imageSCP" style="text-align: center; "
1220
|-
1221
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
1222
<span style="text-align: center; font-size: 75%;">''[[Image:Review_239195285364_7846_Fig_8_circle_contour_temp.png|200px|]]<math>\;</math>[[Image:Review_239195285364_1490_Fig_8_circle_contour_flux.png|200px|]]''</span></div>
1223
|}
1224
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
1225
<span style="text-align: center; font-size: 75%;">''''Figure 5.3:''' Temperature distribution and Flux magnitude fields for the finest mesh.</span></div>
1226
1227
1228
Figures [[#5.1 First example|5.4]](a), [[#5.1 First example|5.4]](b) and [[#5.1 First example|5.4]](c)  show the graphs of the temperature and flux magnitude values along a diameter of the circle for the different meshes of Figures [[#5.1 First example|5.4]](b), [[#5.1 First example|5.4]](c) and [[#5.1 First example|5.4]](d) respectively.
1229
1230
1231
{| class="floating_imageSCP" style="text-align: center; "
1232
|-
1233
| (a)
1234
|[[Image:Review_239195285364_4472_Fig23a_circle.png|260px|]]
1235
|[[Image:Review_239195285364_6864_Fig23a_circle-flux.png|260px|]]
1236
|-
1237
| (b)
1238
|[[Image:Review_239195285364_4411_Fig23b_circle.png|260px|]]
1239
|[[Image:Review_239195285364_6460_Fig23b_circle-flux.png|260px|]]
1240
|-
1241
| (c)
1242
|[[Image:Review_239195285364_8672_Fig23c_circle.png|260px|]]
1243
|[[Image:Review_239195285364_5450_Fig23c_circle-flux.png|260px|]]
1244
|- style="text-align: center; font-size: 75%;"
1245
| colspan="3" style="padding-top:10px;"| '''Figure 5.4:''' Temperature and Flux magnitude graphs along a diameter of the circle for different meshes:      (a) Graphs for the Mesh in Figure [[#5.1 First example|5.4]](b);      (b) Graphs for the Mesh in Figure [[#5.1 First example|5.4]](c);      (c) Graphs for the Mesh in Figure [[#5.1 First example|5.4]](d);
1246
|}
1247
1248
1249
As can be seen from Table [[#5.1 First example|5.1]] and Figure [[#5.1 First example|5.1]], DEC behaves very well on coarse meshes. Note that the maximum temperature values calculated with DEC matches the exact  theoretical value even on coarse meshes. As expected, the results of DEC and FEML are similar for fine meshes.  We would also like to point out the the computational costs of DEC and FEML are very similar.
1250
1251
===5.2 Second example===
1252
1253
In this example, we consider a region in the plane whose boundary consists of segments of a straight line, a circle, a parabola, a cubic and an ellipse (see Figure [[#5.2 Second example|5.5]]).
1254
1255
1256
{| class="floating_imageSCP" style="text-align: center; margin: 1em auto; width: 100%;max-width: 100%;"
1257
|-
1258
|[[Image:Review_239195285364_6878_Fig_10a_patata_prb.png|260px|]]
1259
|[[Image:Review_239195285364_1159_Fig_10b_patata_geom.png|260px|]]
1260
|- style="text-align: center; font-size: 75%;"
1261
| colspan="2" | '''Figure 5.5:''' Region with linear, quadratic and cubic boundary segments, together with boundary conditions.
1262
|}
1263
1264
1265
The meshes used in this example are shown in Figure [[#5.2 Second example|5.6]] and vary from coarse to very fine.
1266
1267
{| class="floating_imageSCP" style="text-align: center; margin: 1em auto; width: 100%;max-width: 100%;"
1268
|-
1269
|[[Image:Review_239195285364_6629_Fig_11c_patata_m2.png|260px|]]
1270
|[[Image:Review_239195285364_5402_Fig_11d_patata_m3.png|260px|]]
1271
|[[Image:Review_239195285364_3774_Fig_11e_patata_m4.png|260px|]]
1272
|-
1273
| (a)
1274
| (b)
1275
| (c)
1276
|- style="text-align: center; font-size: 75%;"
1277
| colspan="3" | '''Figure 5.6:''' Three of the meshes used in the first example.
1278
|}
1279
      
1280
1281
We set
1282
1283
* the source term <math display="inline">q=20.2</math>;
1284
* the heat diffusion constant <math display="inline">k=80.2</math>;
1285
* Newmann boundary condition <math display="inline">h=100</math> along the inner elliptical boundary;
1286
* Dirichlet boundary condition <math display="inline">u=10</math> along the external boundary.
1287
1288
The numerical results for the maximum temperature and flux magnitude values are exemplified in Table [[#5.2 Second example|5.2]].
1289
1290
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1291
|-
1292
| 
1293
{| style="text-align: left; margin:auto;width: 100%;" 
1294
|-
1295
| style="text-align: center;" | <math> \begin{array}{|c|c|c|c|c|c|c|} &  & & &  \\ {Mesh} & \# { nodes} &\# {elements} & {Max. Temp. Value} & {Max. Flux Magnitude}  \\ &&&{DEC \quad FEML}&{DEC \quad FEML} \\
1296
\mbox{ Figure 5.2(a)} &162 &268  & {129.07 \quad 128.92} & {467.22  \quad 467.25}\\
1297
\mbox{ Figure 5.2(b)} &678 &1223 & {129.71 \quad 129.70} & {547.87  \quad 548.42} \\
1298
\mbox{ Figure 5.2(c)} &9489 &18284  & {130.00  \quad 130.00} & {591.58  \quad 591.63} \\
1299
&27453 &53532  & {130.01 \quad 130.01} & {597.37  \quad 597.38}  \\
1300
&651406 &1295960 &{130.02  \quad 130.02} & {602.19  \quad 602.19}  \\ \end{array} </math>
1301
|}
1302
|- style="text-align: center; font-size: 75%;"
1303
| colspan="3" | '''Table 5.2:''' Numerical simulation results.
1304
|}
1305
1306
1307
The temperature and flux-magnitude distribution fields are shown in Figure [[#5.2 Second example|5.7]].
1308
1309
{| class="floating_imageSCP" style="text-align: center; margin: 1em auto; width: 100%;max-width: 100%; "
1310
|-
1311
|[[Image:Review_239195285364_1649_Fig_14_patata2_contour_temp.png|260px]]
1312
|[[Image:Review_239195285364_8633_Fig_14_patata2_contour_flux.png|260px]]
1313
|-
1314
| (a)
1315
| (b)
1316
|- style="text-align: center; font-size: 75%;"
1317
| colspan="3" | '''Figure 5.7:''' Temperature and Flux magnitude distribution fields.
1318
|}
1319
1320
1321
Figure [[#5.2 Second example|5.8]] shows the graphs of the temperature and the flux magnitude along a horizontal line crossing the elliptical boundary for the first two meshes.
1322
1323
{| class="floating_imageSCP" style="text-align: center; margin: 1em auto; width: 100%;max-width: 100%;"
1324
|-
1325
| (a)
1326
|[[Image:Review_239195285364_4170_Fig_15b_patata2_m2.png|260px]]
1327
|[[Image:Review_239195285364_1673_Fig_15b_patata2_m2_flux.png|260px]]
1328
|-
1329
| (b)
1330
|[[Image:Review_239195285364_8948_Fig_15c_patata2_m3.png|260px]]
1331
|[[Image:Review_239195285364_8624_Fig_15c_patata2_m3_flux.png|260px]]
1332
|- style="text-align: center; font-size: 75%;"
1333
| colspan="3" style="padding-top:10px;"| '''Figure 5.8:''' Temperature and Flux magnitude graphs along a horizontal line crossing the elliptical boundary for the first two meshes.
1334
|}
1335
1336
1337
As can be seen from Table [[#5.2 Second example|5.2]] and Figure [[#5.2 Second example|5.8]],  the performance of DEC is very similar to that of FEML in this example.  As expected, when the mesh is refined, the two methods converge to the same values.
1338
1339
==6 Conclusions==
1340
1341
DEC is a relatively recent discretization scheme for PDE which takes into account the geometric and analytic  features of the operators and the domains involved. The main contributions of this paper are the following:
1342
1343
<ol>
1344
1345
<li>We have presented 2D DEC in a simplified manner,  avoiding references to the theory of differential forms and motivating geometrically the new operators.  </li>
1346
<li>We have carried out a numerical comparison between DEC and FEML by solving the 2D Poisson equation  on two cirved domains. The numerical experiments show the solutions obtained with DEC on coarse meshes are as good or better as those of FEML.  On the other hand, the experiements also show numerical convergence. </li>
1347
<li>The computational cost of DEC is similar to that of FEML. </li>
1348
1349
</ol>
1350
1351
==Acknowledgements==
1352
1353
The second named author was partially supported by a grant of CONACYT, and would like to thank the International Centre for Numerical Methods in Engineering (CIMNE) and the University of Swansea for their hospitality. We gratefully acknowledge the support of NVIDIA Corporation with the donation of the Titan X Pascal GPU used for this research.
1354
1355
==Bibliography==
1356
1357
<div class="auto" style="width: auto; margin-left: auto; margin-right: auto;font-size: 85%;">
1358
1359
<div id="cite-1"></div>
1360
[[#citeF-1|[1]]]  D. Arnold, R. Falk, and R. Winther: “Finite element exterior calculus: From Hodge theory to numerical stability”, Bulletin of the American Mathematical Society 47 (2010), no. 2, 281-354.  
1361
1362
<div id="cite-2"></div>
1363
[[#citeF-2|[2]]]  S. Botello, M.A. Moreles, E. Oñate: “Modulo de Aplicaciones del Método de los Elementos Finitos para resolver la ecuación de Poisson: MEFIPOISS.”   Aula CIMNE-CIMAT, Septiembre 2010, ISBN 978-84-96736-95-5.  
1364
1365
<div id="cite-3"></div>
1366
[[#citeF-3|[3]]]  E. Cartan: ``Sur certaines expressions différentielles et le probleme de Pfaff". Annales Scientifiques de l'École Normale Supérieure. Série 3. Paris: Gauthier-Villars. 16: 239-332 (1899)   
1367
1368
<div id="cite-4"></div>
1369
[[#citeF-4|[4]]]   K. Crane, F. de Goes, M. Desbrun, P. Schröder: “Digital geometry processing with discrete exterior calculus.” ACM SIGGRAPH 2013 Courses. ACM, 2013.  
1370
1371
<div id="cite-5"></div>
1372
[[#citeF-5|[5]]]  I. Dassios, A. P. Jivkov, A. Abu-Muharib, P. James:  “A mathematical model for plasticity and damage: A discrete calculus formulation.” Journal of Computational and Applied Mathematics 312 (2017): 27-38.  
1373
1374
<div id="cite-6"></div>
1375
[[#citeF-6|[6]]] M. Desbrun, E. Kanso,  Y. Tong: “Discrete differential forms for computational modeling.” In SIGGRAPH 06: ACM SIGGRAPH 2006 Courses, pages 39&#8211;54, New York, NY, USA, 2006. ACM.  
1376
1377
<div id="cite-7"></div>
1378
[[#citeF-7|[7]]]  M. Griebel, C. Rieger, A. Schier: “Upwind Schemes for Scalar Advection-Dominated Problems in the Discrete Exterior Calculus.” Transport Processes at Fluidic Interfaces. Birkhäuser, Cham, 2017. 145-175.  
1379
1380
<div id="cite-8"></div>
1381
[[#citeF-8|[8]]]  A. N. Hirani: “Discrete exterior calculus”. Diss. California Institute of Technology, 2003.  
1382
1383
<div id="cite-9"></div>
1384
[[#citeF-9|[9]]]  A. N. Hirani, K. B. Nakshatrala, J. H. Chaudhry:  “Numerical method for Darcy flow derived using Discrete Exterior Calculus.” International Journal for Computational Methods in Engineering Science and Mechanics 16.3 (2015): 151-169.  
1385
1386
<div id="cite-10"></div>
1387
[[#citeF-10|[10]]]  A. N. Hirani, K. Kalyanaraman, E. B. VanderZee:  “Delaunay Hodge star”, Comput. Aided Des. 45 (2013) 540-544].  
1388
1389
<div id="cite-11"></div>
1390
[[#citeF-11|[11]]]   M. S. Mohamed, A. N. Hirani, R. Samtaney:  “Discrete exterior calculus discretization of incompressible Navier-Stokes equations over surface simplicial meshes.” Journal of Computational Physics 312 (2016): 175-191.  
1391
1392
<div id="cite-12"></div>
1393
[[#citeF-12|[12]]]  E. Oñate:  “4 - 2D Solids. Linear Triangular and Rectangular Elements,” in Structural Analysis with the Finite Element Method.  Linear Statics, Volume 1: Basis and Solids, CIMNE-Springer, Barcelona, 2009. Pages 117-157, ISBN 978-1-4020-8733-2.    
1394
1395
<div id="cite-13"></div>
1396
[[#citeF-13|[13]]]  O. C. Zienkiewicz, R. L. Taylor and J. Z. Zhu: “3 - Generalization of the finite element concepts. Galerkin-weighted residual and variational approaches,” In The Finite Element Method Set (Sixth Edition), Butterworth-Heinemann, Oxford, 2005, Pages 54-102, ISBN 9780750664318.
1397
</div>
1398

Return to Herrera et al 2018b.

Back to Top

Document information

Published on 08/02/19
Accepted on 04/10/18
Submitted on 06/06/18

Volume 35, Issue 1, 2019
DOI: 10.23967/j.rimni.2018.11.003
Licence: CC BY-NC-SA license

Document Score

0

Times cited: 1
Views 519
Recommendations 0

Share this document

claim authorship

Are you one of the authors of this document?