<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN" "http://www.w3.org/TR/html4/loose.dtd">
<html>
<!-- Created by GNU Texinfo 6.5, http://www.gnu.org/software/texinfo/ -->
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8">
<title>Real Life Example (GNU Octave (version 5.1.0))</title>
<meta name="description" content="Real Life Example (GNU Octave (version 5.1.0))">
<meta name="keywords" content="Real Life Example (GNU Octave (version 5.1.0))">
<meta name="resource-type" content="document">
<meta name="distribution" content="global">
<meta name="Generator" content="makeinfo">
<link href="index.html#Top" rel="start" title="Top">
<link href="Concept-Index.html#Concept-Index" rel="index" title="Concept Index">
<link href="index.html#SEC_Contents" rel="contents" title="Table of Contents">
<link href="Sparse-Matrices.html#Sparse-Matrices" rel="up" title="Sparse Matrices">
<link href="Numerical-Integration.html#Numerical-Integration" rel="next" title="Numerical Integration">
<link href="Iterative-Techniques.html#Iterative-Techniques" rel="prev" title="Iterative Techniques">
<style type="text/css">
<!--
a.summary-letter {text-decoration: none}
blockquote.indentedblock {margin-right: 0em}
blockquote.smallindentedblock {margin-right: 0em; font-size: smaller}
blockquote.smallquotation {font-size: smaller}
div.display {margin-left: 3.2em}
div.example {margin-left: 3.2em}
div.lisp {margin-left: 3.2em}
div.smalldisplay {margin-left: 3.2em}
div.smallexample {margin-left: 3.2em}
div.smalllisp {margin-left: 3.2em}
kbd {font-style: oblique}
pre.display {font-family: inherit}
pre.format {font-family: inherit}
pre.menu-comment {font-family: serif}
pre.menu-preformatted {font-family: serif}
pre.smalldisplay {font-family: inherit; font-size: smaller}
pre.smallexample {font-size: smaller}
pre.smallformat {font-family: inherit; font-size: smaller}
pre.smalllisp {font-size: smaller}
span.nolinebreak {white-space: nowrap}
span.roman {font-family: initial; font-weight: normal}
span.sansserif {font-family: sans-serif; font-weight: normal}
ul.no-bullet {list-style: none}
-->
</style>
<link rel="stylesheet" type="text/css" href="octave.css">
</head>
<body lang="en">
<a name="Real-Life-Example"></a>
<div class="header">
<p>
Previous: <a href="Iterative-Techniques.html#Iterative-Techniques" accesskey="p" rel="prev">Iterative Techniques</a>, Up: <a href="Sparse-Matrices.html#Sparse-Matrices" accesskey="u" rel="up">Sparse Matrices</a> [<a href="index.html#SEC_Contents" title="Table of contents" rel="contents">Contents</a>][<a href="Concept-Index.html#Concept-Index" title="Index" rel="index">Index</a>]</p>
</div>
<hr>
<a name="Real-Life-Example-using-Sparse-Matrices"></a>
<h3 class="section">22.4 Real Life Example using Sparse Matrices</h3>
<p>A common application for sparse matrices is in the solution of Finite
Element Models. Finite element models allow numerical solution of
partial differential equations that do not have closed form solutions,
typically because of the complex shape of the domain.
</p>
<p>In order to motivate this application, we consider the boundary value
Laplace equation. This system can model scalar potential fields, such
as heat or electrical potential. Given a medium
Omega with boundary dOmega. At all points on the dOmega
the boundary conditions are known, and we wish to calculate the potential in
Omega.
Boundary conditions may specify the potential (Dirichlet
boundary condition), its normal derivative across the boundary
(Neumann boundary condition), or a weighted sum of the potential and
its derivative (Cauchy boundary condition).
</p>
<p>In a thermal model, we want to calculate the temperature in
Omega
and know the boundary temperature (Dirichlet condition)
or heat flux (from which we can calculate the Neumann condition
by dividing by the thermal conductivity at the boundary). Similarly,
in an electrical model, we want to calculate the voltage in
Omega
and know the boundary voltage (Dirichlet) or current
(Neumann condition after diving by the electrical conductivity).
In an electrical model, it is common for much of the boundary
to be electrically isolated; this is a Neumann boundary condition
with the current equal to zero.
</p>
<p>The simplest finite element models will divide
Omega
into simplexes (triangles in 2D, pyramids in 3D).
We take as a 3-D example a cylindrical liquid filled tank with a small
non-conductive ball from the EIDORS project<a name="DOCF11" href="#FOOT11"><sup>11</sup></a>. This is model is designed to reflect
an application of electrical impedance tomography, where current patterns
are applied to such a tank in order to image the internal conductivity
distribution. In order to describe the FEM geometry, we have a matrix of
vertices <code>nodes</code> and simplices <code>elems</code>.
</p>
<p>The following example creates a simple rectangular 2-D electrically
conductive medium with 10 V and 20 V imposed on opposite sides
(Dirichlet boundary conditions). All other edges are electrically
isolated.
</p>
<div class="example">
<pre class="example"> node_y = [1;1.2;1.5;1.8;2]*ones(1,11);
node_x = ones(5,1)*[1,1.05,1.1,1.2, ...
1.3,1.5,1.7,1.8,1.9,1.95,2];
nodes = [node_x(:), node_y(:)];
[h,w] = size (node_x);
elems = [];
for idx = 1:w-1
widx = (idx-1)*h;
elems = [elems; ...
widx+[(1:h-1);(2:h);h+(1:h-1)]'; ...
widx+[(2:h);h+(2:h);h+(1:h-1)]' ];
endfor
E = size (elems,1); # No. of simplices
N = size (nodes,1); # No. of vertices
D = size (elems,2); # dimensions+1
</pre></div>
<p>This creates a N-by-2 matrix <code>nodes</code> and a E-by-3 matrix
<code>elems</code> with values, which define finite element triangles:
</p>
<div class="example">
<pre class="example"> nodes(1:7,:)'
1.00 1.00 1.00 1.00 1.00 1.05 1.05 …
1.00 1.20 1.50 1.80 2.00 1.00 1.20 …
elems(1:7,:)'
1 2 3 4 2 3 4 …
2 3 4 5 7 8 9 …
6 7 8 9 6 7 8 …
</pre></div>
<p>Using a first order FEM, we approximate the electrical conductivity
distribution in
Omega
as constant on each simplex (represented by the vector <code>conductivity</code>).
Based on the finite element geometry, we first calculate a system (or
stiffness) matrix for each simplex (represented as 3-by-3 elements on the
diagonal of the element-wise system matrix <code>SE</code>). Based on <code>SE</code>
and a N-by-DE connectivity matrix <code>C</code>, representing the connections
between simplices and vertices, the global connectivity matrix <code>S</code> is
calculated.
</p>
<div class="example">
<pre class="example"> ## Element conductivity
conductivity = [1*ones(1,16), ...
2*ones(1,48), 1*ones(1,16)];
## Connectivity matrix
C = sparse ((1:D*E), reshape (elems', ...
D*E, 1), 1, D*E, N);
## Calculate system matrix
Siidx = floor ([0:D*E-1]'/D) * D * ...
ones(1,D) + ones(D*E,1)*(1:D) ;
Sjidx = [1:D*E]'*ones (1,D);
Sdata = zeros (D*E,D);
dfact = factorial (D-1);
for j = 1:E
a = inv ([ones(D,1), ...
nodes(elems(j,:), :)]);
const = conductivity(j) * 2 / ...
dfact / abs (det (a));
Sdata(D*(j-1)+(1:D),:) = const * ...
a(2:D,:)' * a(2:D,:);
endfor
## Element-wise system matrix
SE = sparse(Siidx,Sjidx,Sdata);
## Global system matrix
S = C'* SE *C;
</pre></div>
<p>The system matrix acts like the conductivity
<code>S</code>
in Ohm’s law
<code>S * V = I</code>.
Based on the Dirichlet and Neumann boundary conditions, we are able
to solve for the voltages at each vertex <code>V</code>.
</p>
<div class="example">
<pre class="example"> ## Dirichlet boundary conditions
D_nodes = [1:5, 51:55];
D_value = [10*ones(1,5), 20*ones(1,5)];
V = zeros (N,1);
V(D_nodes) = D_value;
idx = 1:N; # vertices without Dirichlet
# boundary condns
idx(D_nodes) = [];
## Neumann boundary conditions. Note that
## N_value must be normalized by the
## boundary length and element conductivity
N_nodes = [];
N_value = [];
Q = zeros (N,1);
Q(N_nodes) = N_value;
V(idx) = S(idx,idx) \ ( Q(idx) - ...
S(idx,D_nodes) * V(D_nodes));
</pre></div>
<p>Finally, in order to display the solution, we show each solved voltage
value in the z-axis for each simplex vertex.
See <a href="#fig_003afemmodel">Figure 22.6</a>.
</p>
<div class="example">
<pre class="example"> elemx = elems(:,[1,2,3,1])';
xelems = reshape (nodes(elemx, 1), 4, E);
yelems = reshape (nodes(elemx, 2), 4, E);
velems = reshape (V(elemx), 4, E);
plot3 (xelems,yelems,velems,"k");
print "grid.eps";
</pre></div>
<div class="float"><a name="fig_003afemmodel"></a>
<div align="center"><img src="grid.png" alt="grid">
</div>
<div class="float-caption"><p><strong>Figure 22.6: </strong>Example finite element model the showing triangular elements.
The height of each vertex corresponds to the solution value.</p></div></div>
<div class="footnote">
<hr>
<h4 class="footnotes-heading">Footnotes</h4>
<h3><a name="FOOT11" href="#DOCF11">(11)</a></h3>
<p>EIDORS - Electrical
Impedance Tomography and Diffuse optical Tomography Reconstruction Software
<a href="http://eidors3d.sourceforge.net">http://eidors3d.sourceforge.net</a></p>
</div>
<hr>
<div class="header">
<p>
Previous: <a href="Iterative-Techniques.html#Iterative-Techniques" accesskey="p" rel="prev">Iterative Techniques</a>, Up: <a href="Sparse-Matrices.html#Sparse-Matrices" accesskey="u" rel="up">Sparse Matrices</a> [<a href="index.html#SEC_Contents" title="Table of contents" rel="contents">Contents</a>][<a href="Concept-Index.html#Concept-Index" title="Index" rel="index">Index</a>]</p>
</div>
</body>
</html>
distrib
>
Mageia
>
7
>
i586
>
media
>
core-updates
>
by-pkgid
>
641ebb3060c35990cc021d8f7aaf9aca
>
files
>
444
