-
Notifications
You must be signed in to change notification settings - Fork 3
Expand file tree
/
Copy path5_functions.html
More file actions
391 lines (373 loc) · 37.4 KB
/
Copy path5_functions.html
File metadata and controls
391 lines (373 loc) · 37.4 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
<!DOCTYPE html>
<html lang="en" data-content_root="./">
<head>
<meta charset="utf-8" />
<meta name="viewport" content="width=device-width, initial-scale=1.0" /><meta name="viewport" content="width=device-width, initial-scale=1" />
<title>5. Functions in finite element spaces — Finite element course 2026.0 documentation</title>
<link rel="stylesheet" type="text/css" href="_static/pygments.css?v=03e43079" />
<link rel="stylesheet" type="text/css" href="_static/fenics.css?v=7793b76c" />
<link rel="stylesheet" type="text/css" href="_static/proof.css" />
<link rel="stylesheet" type="text/css" href="_static/sphinx-design.min.css?v=95c83b7e" />
<script src="_static/documentation_options.js?v=9c5a5fc5"></script>
<script src="_static/doctools.js?v=9a2dae69"></script>
<script src="_static/sphinx_highlight.js?v=dc90522c"></script>
<script src="_static/proof.js"></script>
<script src="_static/design-tabs.js?v=36754332"></script>
<script async="async" src="https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML"></script>
<link rel="index" title="Index" href="genindex.html" />
<link rel="search" title="Search" href="search.html" />
<link rel="next" title="6. Assembling and solving finite element problems" href="6_finite_element_problems.html" />
<link rel="prev" title="4. Function spaces: associating data with meshes" href="4_function_spaces.html" />
<!--[if lte IE 6]>
<link rel="stylesheet" href="_static/ie6.css" type="text/css" media="screen" charset="utf-8" />
<![endif]-->
<link rel="stylesheet" href="_static/featured.css">
<link rel="shortcut icon" href="_static/icon.ico" />
</head><body>
<div class="wrapper">
<a href="index.html"><img src="_static/banner.svg" width="900px" alt="FInAT Project Banner" /></a>
<div id="access">
<div class="menu">
<ul>
<li class="page_item"><a href="https://github.com/finite-element/finite-element-course" title="GitHub">GitHub</a></li>
</ul>
</div><!-- .menu -->
</div><!-- #access -->
</div><!-- #wrapper -->
<div class="document">
<div class="documentwrapper">
<div class="bodywrapper">
<div class="body" role="main">
<section id="functions-in-finite-element-spaces">
<span id="secfunctions"></span><h1><span class="section-number">5. </span>Functions in finite element spaces<a class="headerlink" href="#functions-in-finite-element-spaces" title="Link to this heading">¶</a></h1>
<p>Recall that the general form of a function in a finite element space is:</p>
<div class="math notranslate nohighlight" id="equation-function">
<span class="eqno">(5.1)<a class="headerlink" href="#equation-function" title="Link to this equation">¶</a></span>\[f(x) = \sum_i f_i \phi_i(x)\]</div>
<p>Where the <span class="math notranslate nohighlight">\(\phi_i(x)\)</span> are now the global basis functions achieved by
stitching together the local basis functions defined by the
<a class="reference internal" href="2_finite_elements.html#secfinitelement"><span class="std std-ref">finite element</span></a>.</p>
<section id="a-python-implementation-of-functions-in-finite-element-spaces">
<h2><span class="section-number">5.1. </span>A python implementation of functions in finite element spaces<a class="headerlink" href="#a-python-implementation-of-functions-in-finite-element-spaces" title="Link to this heading">¶</a></h2>
<details class="sd-sphinx-override sd-dropdown sd-card sd-mb-3">
<summary class="sd-summary-title sd-card-header">
<span class="sd-summary-text">A video recording of the following material is available here.</span><span class="sd-summary-state-marker sd-summary-chevron-right"><svg version="1.1" width="1.5em" height="1.5em" class="sd-octicon sd-octicon-chevron-right" viewBox="0 0 24 24" aria-hidden="true"><path d="M8.72 18.78a.75.75 0 0 1 0-1.06L14.44 12 8.72 6.28a.751.751 0 0 1 .018-1.042.751.751 0 0 1 1.042-.018l6.25 6.25a.75.75 0 0 1 0 1.06l-6.25 6.25a.75.75 0 0 1-1.06 0Z"></path></svg></span></summary><div class="sd-summary-content sd-card-body docutils">
<div class="vimeo docutils container">
<iframe src="https://player.vimeo.com/video/495447092"
frameborder="0" allow="autoplay; fullscreen"
allowfullscreen></iframe></div>
<p class="sd-card-text">Imperial students can also <a class="reference external" href="https://imperial.cloud.panopto.eu/Panopto/Pages/Viewer.aspx?id=6e7570e9-cc2c-49d2-908b-aca000b63449">watch this video on Panopto</a></p>
</div>
</details><p>The <a class="reference internal" href="fe_utils.html#fe_utils.function_spaces.Function" title="fe_utils.function_spaces.Function"><code class="xref py py-class docutils literal notranslate"><span class="pre">Function</span></code></a> class provides a
simple implementation of function storage. The input is a
<a class="reference internal" href="fe_utils.html#fe_utils.function_spaces.FunctionSpace" title="fe_utils.function_spaces.FunctionSpace"><code class="xref py py-class docutils literal notranslate"><span class="pre">FunctionSpace</span></code></a> which defines the
mesh and finite element to be employed, to which the
<a class="reference internal" href="fe_utils.html#fe_utils.function_spaces.Function" title="fe_utils.function_spaces.Function"><code class="xref py py-class docutils literal notranslate"><span class="pre">Function</span></code></a> adds an array of degree of
freedom values, one for each node in the
<a class="reference internal" href="fe_utils.html#fe_utils.function_spaces.FunctionSpace" title="fe_utils.function_spaces.FunctionSpace"><code class="xref py py-class docutils literal notranslate"><span class="pre">FunctionSpace</span></code></a>.</p>
</section>
<section id="interpolating-values-into-finite-element-spaces">
<h2><span class="section-number">5.2. </span>Interpolating values into finite element spaces<a class="headerlink" href="#interpolating-values-into-finite-element-spaces" title="Link to this heading">¶</a></h2>
<details class="sd-sphinx-override sd-dropdown sd-card sd-mb-3">
<summary class="sd-summary-title sd-card-header">
<span class="sd-summary-text">A video recording of the following material is available here.</span><span class="sd-summary-state-marker sd-summary-chevron-right"><svg version="1.1" width="1.5em" height="1.5em" class="sd-octicon sd-octicon-chevron-right" viewBox="0 0 24 24" aria-hidden="true"><path d="M8.72 18.78a.75.75 0 0 1 0-1.06L14.44 12 8.72 6.28a.751.751 0 0 1 .018-1.042.751.751 0 0 1 1.042-.018l6.25 6.25a.75.75 0 0 1 0 1.06l-6.25 6.25a.75.75 0 0 1-1.06 0Z"></path></svg></span></summary><div class="sd-summary-content sd-card-body docutils">
<div class="vimeo docutils container">
<iframe src="https://player.vimeo.com/video/495447659"
frameborder="0" allow="autoplay; fullscreen"
allowfullscreen></iframe></div>
<p class="sd-card-text">Imperial students can also <a class="reference external" href="https://imperial.cloud.panopto.eu/Panopto/Pages/Viewer.aspx?id=6d13d669-dc8c-43d1-8aae-aca000b634ed">watch this video on Panopto</a></p>
</div>
</details><p>Suppose we have a function <span class="math notranslate nohighlight">\(g(x): \mathbb{R}^n \rightarrow \mathbb{R}\)</span>
which we wish to approximate as a function <span class="math notranslate nohighlight">\(f(x)\)</span> in some finite
element space <span class="math notranslate nohighlight">\(V\)</span>. In other words, we want to find the <span class="math notranslate nohighlight">\(f_i\)</span> such that:</p>
<div class="math notranslate nohighlight" id="equation-5-functions-0">
<span class="eqno">(5.2)<a class="headerlink" href="#equation-5-functions-0" title="Link to this equation">¶</a></span>\[\sum_i f_i \phi_i(x) \approx g(x)\]</div>
<p>The simplest way to do this is to <em>interpolate</em> <span class="math notranslate nohighlight">\(g(x)\)</span> onto <span class="math notranslate nohighlight">\(V\)</span>. In
other words, we evaluate:</p>
<div class="math notranslate nohighlight" id="equation-5-functions-1">
<span class="eqno">(5.3)<a class="headerlink" href="#equation-5-functions-1" title="Link to this equation">¶</a></span>\[f_i = n_i(g(x))\]</div>
<p>where <span class="math notranslate nohighlight">\(n_i\)</span> is the node associated with <span class="math notranslate nohighlight">\(\phi_i\)</span> as considered over the
entire mesh (<em>globally</em>). Since we are only concerned with point
evaluation nodes, this is equivalent to:</p>
<div class="math notranslate nohighlight" id="equation-5-functions-2">
<span class="eqno">(5.4)<a class="headerlink" href="#equation-5-functions-2" title="Link to this equation">¶</a></span>\[f_i = g(x_i)\]</div>
<p>where <span class="math notranslate nohighlight">\(x_i\)</span> is the coordinate vector of the point defining the node
<span class="math notranslate nohighlight">\(n_i\)</span>. This looks straightforward, however the <span class="math notranslate nohighlight">\(x_i\)</span> are the <em>global</em>
node points, and so far we have only defined the node points in
<em>local</em> coordinates on the reference element.</p>
<section id="changing-coordinates-between-reference-and-physical-space">
<span id="coordinates"></span><h3><span class="section-number">5.2.1. </span>Changing coordinates between reference and physical space<a class="headerlink" href="#changing-coordinates-between-reference-and-physical-space" title="Link to this heading">¶</a></h3>
<p>We’ll refer to coordinates on the global mesh as being in <em>physical
space</em> while those on the reference element are in <em>local
space</em>. We’ll use case to distinguish local and global objects, so
local coordinates will be written as <span class="math notranslate nohighlight">\(X\)</span> and global coordinates as
<span class="math notranslate nohighlight">\(x\)</span>. The key observation is that within each cell, the global
coordinates are the linear interpolation of the global coordinate
values at the cell vertices. In other words, if <span class="math notranslate nohighlight">\(\{\Psi_j\}\)</span> is the
local basis for the <strong>linear</strong> lagrange elements on the reference cell and
<span class="math notranslate nohighlight">\(\hat{x}_j\)</span> are the corresponding global vertex locations on a cell <span class="math notranslate nohighlight">\(c\)</span>
then:</p>
<div class="math notranslate nohighlight" id="equation-change">
<span class="eqno">(5.5)<a class="headerlink" href="#equation-change" title="Link to this equation">¶</a></span>\[x = \sum_j \hat{x}_j \Psi_j(X) \quad \forall x \in c.\]</div>
<p>Remember that we know the location of the nodes in local coordinates,
and we have the
<a class="reference internal" href="fe_utils.html#fe_utils.finite_elements.FiniteElement.tabulate" title="fe_utils.finite_elements.FiniteElement.tabulate"><code class="xref py py-meth docutils literal notranslate"><span class="pre">tabulate()</span></code></a> method to
evaluate all the basis functions of an element at a known set of
points. So if we write:</p>
<div class="math notranslate nohighlight" id="equation-foo0">
<span class="eqno">(5.6)<a class="headerlink" href="#equation-foo0" title="Link to this equation">¶</a></span>\[A_{i,j} = \Psi_j(X_i)\]</div>
<p>where <span class="math notranslate nohighlight">\(\{X_i\}\)</span> are the local node points of our finite element, then
the corresponding global node points <span class="math notranslate nohighlight">\(\{x_i\}\)</span> are given by</p>
<div class="math notranslate nohighlight" id="equation-foo1">
<span class="eqno">(5.7)<a class="headerlink" href="#equation-foo1" title="Link to this equation">¶</a></span>\[x = A\cdot \hat{x}\]</div>
<p>where <span class="math notranslate nohighlight">\(\hat{x}\)</span> is the <span class="math notranslate nohighlight">\((\dim+1, \dim)\)</span> array whose rows are the current
element vertex coordinates. Specifically, <span class="math notranslate nohighlight">\(x\)</span> is the <span class="math notranslate nohighlight">\((\textrm{nodes}, \dim)\)</span> array whose
rows are the global coordinates of the nodes in the current
element. We can then apply <span class="math notranslate nohighlight">\(g()\)</span> to each row of <span class="math notranslate nohighlight">\(x\)</span> in turn and record
the result as the <a class="reference internal" href="fe_utils.html#fe_utils.function_spaces.Function" title="fe_utils.function_spaces.Function"><code class="xref py py-class docutils literal notranslate"><span class="pre">Function</span></code></a> value
for that node.</p>
<div class="admonition hint">
<p class="admonition-title">Hint</p>
<p>The observant reader will notice that this algorithm is inefficient
because the function values at nodes on the boundaries of elements
are evaluated more than once. This can be avoided with a little
tedious bookkeeping but we will not concern ourselves with that
here.</p>
</div>
</section>
<section id="looking-up-cell-coordinates-and-values">
<h3><span class="section-number">5.2.2. </span>Looking up cell coordinates and values<a class="headerlink" href="#looking-up-cell-coordinates-and-values" title="Link to this heading">¶</a></h3>
<details class="sd-sphinx-override sd-dropdown sd-card sd-mb-3">
<summary class="sd-summary-title sd-card-header">
<span class="sd-summary-text">A video recording of the following material is available here.</span><span class="sd-summary-state-marker sd-summary-chevron-right"><svg version="1.1" width="1.5em" height="1.5em" class="sd-octicon sd-octicon-chevron-right" viewBox="0 0 24 24" aria-hidden="true"><path d="M8.72 18.78a.75.75 0 0 1 0-1.06L14.44 12 8.72 6.28a.751.751 0 0 1 .018-1.042.751.751 0 0 1 1.042-.018l6.25 6.25a.75.75 0 0 1 0 1.06l-6.25 6.25a.75.75 0 0 1-1.06 0Z"></path></svg></span></summary><div class="sd-summary-content sd-card-body docutils">
<div class="vimeo docutils container">
<iframe src="https://player.vimeo.com/video/495715631"
frameborder="0" allow="autoplay; fullscreen"
allowfullscreen></iframe></div>
<p class="sd-card-text">Imperial students can also <a class="reference external" href="https://imperial.cloud.panopto.eu/Panopto/Pages/Viewer.aspx?id=eb8b0e40-b476-4caa-a961-aca1009f6172">watch this video on Panopto</a></p>
</div>
</details><p>In the previous section we used the vertex coordinates of a cell to
find the node coordinates, and then we calculated
<a class="reference internal" href="fe_utils.html#fe_utils.function_spaces.Function" title="fe_utils.function_spaces.Function"><code class="xref py py-class docutils literal notranslate"><span class="pre">Function</span></code></a> values at those
points. The coordinates are stored in a single long list associated
with the <a class="reference internal" href="fe_utils.html#fe_utils.mesh.Mesh" title="fe_utils.mesh.Mesh"><code class="xref py py-class docutils literal notranslate"><span class="pre">Mesh</span></code></a>, and the
<a class="reference internal" href="fe_utils.html#fe_utils.function_spaces.Function" title="fe_utils.function_spaces.Function"><code class="xref py py-class docutils literal notranslate"><span class="pre">Function</span></code></a> contains a single long
list of values. We need to use <em>indirect addressing</em> to access these
values. This is best illustrated using some Python code.</p>
<p>Suppose <code class="docutils literal notranslate"><span class="pre">f</span></code> is a <a class="reference internal" href="fe_utils.html#fe_utils.function_spaces.Function" title="fe_utils.function_spaces.Function"><code class="xref py py-class docutils literal notranslate"><span class="pre">Function</span></code></a>.
For brevity, we write <code class="docutils literal notranslate"><span class="pre">fs</span> <span class="pre">=</span> <span class="pre">f.function_space</span></code>, the
<a class="reference internal" href="fe_utils.html#fe_utils.function_spaces.FunctionSpace" title="fe_utils.function_spaces.FunctionSpace"><code class="xref py py-class docutils literal notranslate"><span class="pre">FunctionSpace</span></code></a> associated with
<code class="docutils literal notranslate"><span class="pre">f</span></code>. Now, we first need a linear element and a corresponding
<a class="reference internal" href="fe_utils.html#fe_utils.function_spaces.FunctionSpace" title="fe_utils.function_spaces.FunctionSpace"><code class="xref py py-class docutils literal notranslate"><span class="pre">FunctionSpace</span></code></a>:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="n">cg1</span> <span class="o">=</span> <span class="n">fe_utils</span><span class="o">.</span><span class="n">LagrangeElement</span><span class="p">(</span><span class="n">fs</span><span class="o">.</span><span class="n">mesh</span><span class="o">.</span><span class="n">cell</span><span class="p">,</span> <span class="mi">1</span><span class="p">)</span>
<span class="n">cg1fs</span> <span class="o">=</span> <span class="n">fe_utils</span><span class="o">.</span><span class="n">FunctionSpace</span><span class="p">(</span><span class="n">fs</span><span class="o">.</span><span class="n">mesh</span><span class="p">,</span> <span class="n">cg1</span><span class="p">)</span>
</pre></div>
</div>
<p>Then the vertex indices of cell number <code class="docutils literal notranslate"><span class="pre">c</span></code> in the correct order for the linear Lagrange element are:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="n">cg1fs</span><span class="o">.</span><span class="n">cell_nodes</span><span class="p">[</span><span class="n">c</span><span class="p">,</span> <span class="p">:]</span>
</pre></div>
</div>
<p>and therefore the set of coordinate vectors for the vertices of
element <code class="docutils literal notranslate"><span class="pre">c</span></code> are:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="n">fs</span><span class="o">.</span><span class="n">mesh</span><span class="o">.</span><span class="n">vertex_coords</span><span class="p">[</span><span class="n">cg1fs</span><span class="o">.</span><span class="n">cell_nodes</span><span class="p">[</span><span class="n">c</span><span class="p">,</span> <span class="p">:],</span> <span class="p">:]</span>
</pre></div>
</div>
<p>That is, the <code class="docutils literal notranslate"><span class="pre">cg1fs.cell_nodes</span></code> array is used to look up the right
vertex coordinates. By a similar process we can access the values
associated with the nodes of element <code class="docutils literal notranslate"><span class="pre">c</span></code>:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="n">f</span><span class="o">.</span><span class="n">values</span><span class="p">[</span><span class="n">fs</span><span class="o">.</span><span class="n">cell_nodes</span><span class="p">[</span><span class="n">c</span><span class="p">,</span> <span class="p">:]]</span>
</pre></div>
</div>
</section>
<section id="a-python-implementation-of-interpolation">
<h3><span class="section-number">5.2.3. </span>A Python implementation of interpolation<a class="headerlink" href="#a-python-implementation-of-interpolation" title="Link to this heading">¶</a></h3>
<details class="sd-sphinx-override sd-dropdown sd-card sd-mb-3">
<summary class="sd-summary-title sd-card-header">
<span class="sd-summary-text">A video recording of the following material is available here.</span><span class="sd-summary-state-marker sd-summary-chevron-right"><svg version="1.1" width="1.5em" height="1.5em" class="sd-octicon sd-octicon-chevron-right" viewBox="0 0 24 24" aria-hidden="true"><path d="M8.72 18.78a.75.75 0 0 1 0-1.06L14.44 12 8.72 6.28a.751.751 0 0 1 .018-1.042.751.751 0 0 1 1.042-.018l6.25 6.25a.75.75 0 0 1 0 1.06l-6.25 6.25a.75.75 0 0 1-1.06 0Z"></path></svg></span></summary><div class="sd-summary-content sd-card-body docutils">
<div class="vimeo docutils container">
<iframe src="https://player.vimeo.com/video/495718032"
frameborder="0" allow="autoplay; fullscreen"
allowfullscreen></iframe></div>
<p class="sd-card-text">Imperial students can also <a class="reference external" href="https://imperial.cloud.panopto.eu/Panopto/Pages/Viewer.aspx?id=e74afdcc-7727-49cb-bf3d-aca100a2ec68">watch this video on Panopto</a></p>
</div>
</details><p>Putting together the change of coordinates with the right indirect
addressing, we can provide the
<a class="reference internal" href="fe_utils.html#fe_utils.function_spaces.Function" title="fe_utils.function_spaces.Function"><code class="xref py py-class docutils literal notranslate"><span class="pre">Function</span></code></a> class with a
<a class="reference internal" href="fe_utils.html#fe_utils.function_spaces.Function.interpolate" title="fe_utils.function_spaces.Function.interpolate"><code class="xref py py-meth docutils literal notranslate"><span class="pre">interpolate()</span></code></a> method which
interpolates a user-provided function onto the
<a class="reference internal" href="fe_utils.html#fe_utils.function_spaces.Function" title="fe_utils.function_spaces.Function"><code class="xref py py-class docutils literal notranslate"><span class="pre">Function</span></code></a>.</p>
<div class="proof proof-type-exercise" id="id2">
<div class="proof-title">
<span class="proof-type">Exercise 5.37</span>
</div><div class="proof-content">
<p>Read and understand the
<a class="reference internal" href="fe_utils.html#fe_utils.function_spaces.Function.interpolate" title="fe_utils.function_spaces.Function.interpolate"><code class="xref py py-meth docutils literal notranslate"><span class="pre">interpolate()</span></code></a> method. Use
<code class="docutils literal notranslate"><span class="pre">plot_sin_function</span></code> to investigate interpolating different
functions onto finite element spaces at differering resolutions and
polynomial degrees.</p>
</div></div><div class="admonition hint">
<p class="admonition-title">Hint</p>
<p>There is no implementation work associated with this exercise, but
the programming constructs used in
<a class="reference internal" href="fe_utils.html#fe_utils.function_spaces.Function.interpolate" title="fe_utils.function_spaces.Function.interpolate"><code class="xref py py-meth docutils literal notranslate"><span class="pre">interpolate()</span></code></a> will be
needed when you implement integration.</p>
</div>
</section>
</section>
<section id="integration">
<span id="id1"></span><h2><span class="section-number">5.3. </span>Integration<a class="headerlink" href="#integration" title="Link to this heading">¶</a></h2>
<details class="sd-sphinx-override sd-dropdown sd-card sd-mb-3">
<summary class="sd-summary-title sd-card-header">
<span class="sd-summary-text">A video recording of the following material is available here.</span><span class="sd-summary-state-marker sd-summary-chevron-right"><svg version="1.1" width="1.5em" height="1.5em" class="sd-octicon sd-octicon-chevron-right" viewBox="0 0 24 24" aria-hidden="true"><path d="M8.72 18.78a.75.75 0 0 1 0-1.06L14.44 12 8.72 6.28a.751.751 0 0 1 .018-1.042.751.751 0 0 1 1.042-.018l6.25 6.25a.75.75 0 0 1 0 1.06l-6.25 6.25a.75.75 0 0 1-1.06 0Z"></path></svg></span></summary><div class="sd-summary-content sd-card-body docutils">
<div class="vimeo docutils container">
<iframe src="https://player.vimeo.com/video/495447232"
frameborder="0" allow="autoplay; fullscreen"
allowfullscreen></iframe></div>
<p class="sd-card-text">Imperial students can also <a class="reference external" href="https://imperial.cloud.panopto.eu/Panopto/Pages/Viewer.aspx?id=4b2f4755-0682-4e36-ad3e-aca000b634b9">watch this video on Panopto</a></p>
</div>
</details><p>We now come to one of the fundamental operations in the finite element
method: integrating a <a class="reference internal" href="fe_utils.html#fe_utils.function_spaces.Function" title="fe_utils.function_spaces.Function"><code class="xref py py-class docutils literal notranslate"><span class="pre">Function</span></code></a> over
the domain. The full finite element method actually requires the
integration of expressions of unknown test and trial functions, but we
will start with the more straightforward case of integrating a single,
known, <a class="reference internal" href="fe_utils.html#fe_utils.function_spaces.Function" title="fe_utils.function_spaces.Function"><code class="xref py py-class docutils literal notranslate"><span class="pre">Function</span></code></a> over a domain
<span class="math notranslate nohighlight">\(\Omega\)</span>:</p>
<div class="math notranslate nohighlight" id="equation-5-functions-3">
<span class="eqno">(5.8)<a class="headerlink" href="#equation-5-functions-3" title="Link to this equation">¶</a></span>\[\int_\Omega f \mathrm{d} x \quad f \in V\]</div>
<p>where <span class="math notranslate nohighlight">\(\mathrm{d}x\)</span> should be understood as being the volume measure
with the correct dimension for the domain and <span class="math notranslate nohighlight">\(V\)</span> is some finite
element space over <span class="math notranslate nohighlight">\(\Omega\)</span>. We can express this integral as a sum of
integrals over individual cells:</p>
<div class="math notranslate nohighlight" id="equation-integral-sum">
<span class="eqno">(5.9)<a class="headerlink" href="#equation-integral-sum" title="Link to this equation">¶</a></span>\[\int_\Omega f \mathrm{d} x = \sum_{c\in\Omega} \int_c f \mathrm{d} x.\]</div>
<p>So we have in fact reduced the integration problem to the problem of
integrating <span class="math notranslate nohighlight">\(f\)</span> over each cell. In <a class="reference internal" href="1_quadrature.html"><span class="doc">a previous part</span></a>
of the module we implemented quadrature rules which enable us to
integrate over specified reference cells. If we can express the
integral over some arbitrary cell <span class="math notranslate nohighlight">\(c\)</span> as an integral over a reference
cell <span class="math notranslate nohighlight">\(c_0\)</span> then we are done. In fact this simply requires us to employ
the change of variables formula for integration:</p>
<div class="math notranslate nohighlight" id="equation-5-functions-4">
<span class="eqno">(5.10)<a class="headerlink" href="#equation-5-functions-4" title="Link to this equation">¶</a></span>\[\int_{c} f(x) \mathrm{d} x = \int_{c_0} f(X) |J|\mathrm{d} X\]</div>
<p>where <span class="math notranslate nohighlight">\(|J|\)</span> is the absolute value of the determinant of the Jacobian
matrix. <span class="math notranslate nohighlight">\(J\)</span> is given by:</p>
<div class="math notranslate nohighlight" id="equation-jacobian-def">
<span class="eqno">(5.11)<a class="headerlink" href="#equation-jacobian-def" title="Link to this equation">¶</a></span>\[J_{\alpha\beta} = \frac{\partial x_\alpha}{\partial X_\beta}.\]</div>
<div class="admonition hint">
<p class="admonition-title">Hint</p>
<p>We will generally adopt the convention of using Greek letters to
indicate indices in spatial dimensions, while we will use Roman
letters in the sequence <span class="math notranslate nohighlight">\(i,j,\ldots\)</span> for basis function indices. We
will continue to use <span class="math notranslate nohighlight">\(q\)</span> for the index over the quadrature points.</p>
</div>
<p>Evaluating <a class="reference internal" href="#equation-jacobian-def">(5.11)</a> depends on having an expression for <span class="math notranslate nohighlight">\(x\)</span> in
terms of <span class="math notranslate nohighlight">\(X\)</span>. Fortunately, <a class="reference internal" href="#equation-change">(5.5)</a> is exactly this expression,
and applying the usual rule for differentiating functions in finite
element spaces produces:</p>
<div class="math notranslate nohighlight" id="equation-jacobian">
<span class="eqno">(5.12)<a class="headerlink" href="#equation-jacobian" title="Link to this equation">¶</a></span>\[J_{\alpha\beta} = \sum_j (\tilde{x}_j)_\alpha \nabla_\beta\Psi_j(X)\]</div>
<p>where <span class="math notranslate nohighlight">\(\{\Psi_j\}\)</span> is once again the degree 1 Lagrange basis and
<span class="math notranslate nohighlight">\(\{\tilde{x}_j\}\)</span> are the coordinates of the corresponding vertices of
cell <span class="math notranslate nohighlight">\(c\)</span>. The presence of <span class="math notranslate nohighlight">\(X\)</span> in <a class="reference internal" href="#equation-jacobian">(5.12)</a> implies that the
Jacobian varies spatially across the reference cell. However since
<span class="math notranslate nohighlight">\(\{\Psi_j\}\)</span> is the degree 1 Lagrange basis, the gradients of the
basis functions are constant over the cell and so it does not matter
at which point in the cell the Jacobian is evaluated. For example we
might choose to evaluate the Jacobian at the cell origin <span class="math notranslate nohighlight">\(X=0\)</span>.</p>
<div class="admonition hint">
<p class="admonition-title">Hint</p>
<p>When using simplices with curved sides, and on all but the simplest
quadrilateral or hexahedral meshes, the change of coordinates
will not be affine. In that case, to preserve full accuracy it will be
necessary to compute the Jacobian at every quadrature
point. However, non-affine coordinate transforms are beyond the
scope of this course.</p>
</div>
<section id="expressing-the-function-in-the-finite-element-basis">
<h3><span class="section-number">5.3.1. </span>Expressing the function in the finite element basis<a class="headerlink" href="#expressing-the-function-in-the-finite-element-basis" title="Link to this heading">¶</a></h3>
<details class="sd-sphinx-override sd-dropdown sd-card sd-mb-3">
<summary class="sd-summary-title sd-card-header">
<span class="sd-summary-text">A video recording of the following material is available here.</span><span class="sd-summary-state-marker sd-summary-chevron-right"><svg version="1.1" width="1.5em" height="1.5em" class="sd-octicon sd-octicon-chevron-right" viewBox="0 0 24 24" aria-hidden="true"><path d="M8.72 18.78a.75.75 0 0 1 0-1.06L14.44 12 8.72 6.28a.751.751 0 0 1 .018-1.042.751.751 0 0 1 1.042-.018l6.25 6.25a.75.75 0 0 1 0 1.06l-6.25 6.25a.75.75 0 0 1-1.06 0Z"></path></svg></span></summary><div class="sd-summary-content sd-card-body docutils">
<div class="vimeo docutils container">
<iframe src="https://player.vimeo.com/video/495444345"
frameborder="0" allow="autoplay; fullscreen"
allowfullscreen></iframe></div>
<p class="sd-card-text">Imperial students can also <a class="reference external" href="https://imperial.cloud.panopto.eu/Panopto/Pages/Viewer.aspx?id=005f670f-28ac-47c3-9cf1-aca000b7133d">watch this video on Panopto</a></p>
</div>
</details><p>Let <span class="math notranslate nohighlight">\(\{\Phi_i(X)\}\)</span> be a <strong>local</strong> basis for <span class="math notranslate nohighlight">\(V\)</span> on the reference element
<span class="math notranslate nohighlight">\(c_0\)</span>. Then our integral becomes:</p>
<div class="math notranslate nohighlight" id="equation-5-functions-5">
<span class="eqno">(5.13)<a class="headerlink" href="#equation-5-functions-5" title="Link to this equation">¶</a></span>\[\int_c f(x)\mathrm{d}x = \int_{c_0} \sum_i F(M(c,i))\,\Phi_i(X)\, |J|\,\mathrm{d} X\]</div>
<p>where <span class="math notranslate nohighlight">\(F\)</span> is the vector of global coefficient values of <span class="math notranslate nohighlight">\(f\)</span>, and <span class="math notranslate nohighlight">\(M\)</span> is <a class="reference internal" href="4_function_spaces.html#cell-node"><span class="std std-ref">the cell node map</span></a>.</p>
</section>
<section id="numerical-quadrature">
<h3><span class="section-number">5.3.2. </span>Numerical quadrature<a class="headerlink" href="#numerical-quadrature" title="Link to this heading">¶</a></h3>
<p>The actual evaluation of the integral will employ the quadrature rules
we discussed in <a class="reference internal" href="1_quadrature.html"><span class="doc">a previous section</span></a>. Let
<span class="math notranslate nohighlight">\(\{X_q\}, \{w_q\}\)</span> be a quadrature rule of sufficient degree of precision that
the quadrature is exact. Then:</p>
<div class="math notranslate nohighlight" id="equation-integration">
<span class="eqno">(5.14)<a class="headerlink" href="#equation-integration" title="Link to this equation">¶</a></span>\[\int_c f(x)\mathrm{d}x = \sum_q \sum_i F(M(c,i))\,\Phi_i(X_q)\, |J|\,w_q\]</div>
</section>
<section id="implementing-integration">
<h3><span class="section-number">5.3.3. </span>Implementing integration<a class="headerlink" href="#implementing-integration" title="Link to this heading">¶</a></h3>
<div class="proof proof-type-exercise" id="id3">
<div class="proof-title">
<span class="proof-type">Exercise 5.38</span>
</div><div class="proof-content">
<p>Use <a class="reference internal" href="#equation-jacobian">(5.12)</a> to implement the
<a class="reference internal" href="fe_utils.html#fe_utils.mesh.Mesh.jacobian" title="fe_utils.mesh.Mesh.jacobian"><code class="xref py py-meth docutils literal notranslate"><span class="pre">jacobian()</span></code></a> method of
<a class="reference internal" href="fe_utils.html#fe_utils.mesh.Mesh" title="fe_utils.mesh.Mesh"><code class="xref py py-class docutils literal notranslate"><span class="pre">Mesh</span></code></a>. <code class="docutils literal notranslate"><span class="pre">test/test_09_jacobian.py</span></code> is
available for you to test your results.</p>
</div></div><div class="admonition hint">
<p class="admonition-title">Hint</p>
<p>The <span class="math notranslate nohighlight">\(\nabla_\beta\Psi_j(X)\)</span> factor in <a class="reference internal" href="#equation-jacobian">(5.12)</a> is the same
for every cell in the mesh. You could make your implementation more
efficient by precalculating this term in the <code class="xref py py-meth docutils literal notranslate"><span class="pre">__init__()</span></code>
method of <a class="reference internal" href="fe_utils.html#fe_utils.mesh.Mesh" title="fe_utils.mesh.Mesh"><code class="xref py py-class docutils literal notranslate"><span class="pre">Mesh</span></code></a>.</p>
</div>
<div class="proof proof-type-exercise" id="id4">
<div class="proof-title">
<span class="proof-type">Exercise 5.39</span>
</div><div class="proof-content">
<p>Use <a class="reference internal" href="#equation-integral-sum">(5.9)</a> and <a class="reference internal" href="#equation-integration">(5.14)</a> to implement
<a class="reference internal" href="fe_utils.html#fe_utils.function_spaces.Function.integrate" title="fe_utils.function_spaces.Function.integrate"><code class="xref py py-meth docutils literal notranslate"><span class="pre">integrate()</span></code></a>.
<code class="docutils literal notranslate"><span class="pre">test/test_10_integrate_function.py</span></code> may be used to test your
implementation.</p>
</div></div><div class="admonition hint">
<p class="admonition-title">Hint</p>
<p>Your method will need to:</p>
<ol class="arabic simple">
<li><p>Construct a suitable <a class="reference internal" href="fe_utils.html#fe_utils.quadrature.QuadratureRule" title="fe_utils.quadrature.QuadratureRule"><code class="xref py py-class docutils literal notranslate"><span class="pre">QuadratureRule</span></code></a>.</p></li>
<li><p><a class="reference internal" href="fe_utils.html#fe_utils.finite_elements.FiniteElement.tabulate" title="fe_utils.finite_elements.FiniteElement.tabulate"><code class="xref py py-meth docutils literal notranslate"><span class="pre">tabulate()</span></code></a> the
basis functions at each quadrature point.</p></li>
<li><p>Visit each cell in turn.</p></li>
<li><p>Construct the <a class="reference internal" href="fe_utils.html#fe_utils.mesh.Mesh.jacobian" title="fe_utils.mesh.Mesh.jacobian"><code class="xref py py-meth docutils literal notranslate"><span class="pre">jacobian()</span></code></a> for that cell
and take the absolute value of its determinant (<a class="reference external" href="https://numpy.org/doc/stable/reference/generated/numpy.absolute.html#numpy.absolute" title="(in NumPy v2.4)"><code class="xref py py-data docutils literal notranslate"><span class="pre">numpy.absolute</span></code></a>
and <a class="reference external" href="https://numpy.org/doc/stable/reference/generated/numpy.linalg.det.html#numpy.linalg.det" title="(in NumPy v2.4)"><code class="xref py py-func docutils literal notranslate"><span class="pre">numpy.linalg.det()</span></code></a> will be useful here).</p></li>
<li><p>Sum all of the arrays you have constructed over the correct
indices to a contribution to the integral (<a class="reference external" href="https://numpy.org/doc/stable/reference/generated/numpy.einsum.html#numpy.einsum" title="(in NumPy v2.4)"><code class="xref py py-func docutils literal notranslate"><span class="pre">numpy.einsum()</span></code></a>
may be useful for this).</p></li>
</ol>
</div>
<div class="admonition hint">
<p class="admonition-title">Hint</p>
<p>You might choose to read ahead before implementing
<a class="reference internal" href="fe_utils.html#fe_utils.function_spaces.Function.integrate" title="fe_utils.function_spaces.Function.integrate"><code class="xref py py-meth docutils literal notranslate"><span class="pre">integrate()</span></code></a>, since the
<a class="reference internal" href="fe_utils.html#fe_utils.utils.errornorm" title="fe_utils.utils.errornorm"><code class="xref py py-func docutils literal notranslate"><span class="pre">errornorm()</span></code></a> function is very similar and may provide a useful
template for your work.</p>
</div>
</section>
</section>
</section>
<div class="clearer"></div>
</div>
</div>
</div>
<div class="clearer"></div>
</div>
<div class="footer" role="contentinfo">
© Copyright 2014-2026, David A. Ham and Colin J. Cotter.
Created using <a href="https://www.sphinx-doc.org/">Sphinx</a> 7.4.7.
</div>
</body>
</html>