1 |
% $Id$ |
2 |
|
3 |
|
4 |
\chapter{ The module \finley} |
5 |
\label{CHAPTER ON FINLEY} |
6 |
|
7 |
\begin{figure} |
8 |
\centerline{\includegraphics[width=\figwidth]{FinleyMesh}} |
9 |
\caption{Subdivision of an Ellipse into triangles order 1 (\finleyelement{Tri3})} |
10 |
\label{FINLEY FIG 0} |
11 |
\end{figure} |
12 |
|
13 |
\begin{figure} |
14 |
\centerline{\includegraphics[width=\figwidth]{FinleyContact}} |
15 |
\caption{Mesh around a contact region (\finleyelement{Rec4})} |
16 |
\label{FINLEY FIG 01} |
17 |
\end{figure} |
18 |
|
19 |
\declaremodule{extension}{finley} \modulesynopsis{Solving linear, steady partial differential equations using |
20 |
finite elements} |
21 |
|
22 |
{\it finley} is a library of C functions solving linear, steady partial differential equations |
23 |
\index{partial differential equations} (PDEs) or systems of PDEs using isoparametrical finite |
24 |
elements \index{FEM!isoparametrical}. |
25 |
It supports unstructured, 1D, 2D and 3D meshes. The module \finley provides an access to the |
26 |
library through the \LinearPDE class of \escript supporting its full functionality. {\it finley} |
27 |
is parallelized using the OpenMP \index{OpenMP} paradigm. |
28 |
|
29 |
\subsection{Meshes} |
30 |
To understand the usage of \finley one needs to have an understanding of how the finite element meshes |
31 |
\index{FEM!mesh} are defined. \fig{FINLEY FIG 0} shows an example of the |
32 |
subdivision of an ellipse into so called elements \index{FEM!elements} \index{element}. |
33 |
In this case, triangles have been used but other forms of subdivisions |
34 |
can be constructed, e.g. into quadrilaterals or, in the three dimensional case, into tetrahedrons |
35 |
and hexahedrons. The idea of the finite element method is to approximate the solution by a function |
36 |
which is a polynomial of a certain order and is continuous across it boundary to neighbour elements. |
37 |
In the example of \fig{FINLEY FIG 0} a linear polynomial is used on each triangle. As one can see, the triangulation |
38 |
is quite a poor approximation of the ellipse. It can be improved by introducing a midpoint on each element edge then |
39 |
positioning those nodes located on an edge expected to describe the boundary, onto the boundary. |
40 |
In this case the triangle gets a curved edge which requires a parametrization of the triangle using a |
41 |
quadratic polynomial. For this case, the solution is also approximated by a piecewise quadratic polynomial |
42 |
(which explains the name isoparametrical elements), see \Ref{Zienc,NumHand} for more details. |
43 |
|
44 |
The union of all elements defines the domain of the PDE. |
45 |
Each element is defined by the nodes used to describe its shape. In \fig{FINLEY FIG 0} the element, |
46 |
which has type \finleyelement{Tri3}, |
47 |
with element reference number $19$ \index{element!reference number} is defined by the nodes |
48 |
with reference numbers $9$, $11$ and $0$ \index{node!reference number}. Notice that the order is counterclockwise. |
49 |
The coefficients of the PDE are evaluated at integration nodes with each individual element. |
50 |
For quadrilateral elements a Gauss quadrature scheme is used. In the case of triangular elements a |
51 |
modified form is applied. The boundary of the domain is also subdivided into elements. \index{element!face} In \fig{FINLEY FIG 0} |
52 |
line elements with two nodes are used. The elements are also defined by their describing nodes, e.g. |
53 |
the face element reference number $20$ which has type \finleyelement{Line2} is defined by the nodes |
54 |
with the reference numbers $11$ and $0$. Again the order is crucial, if moving from the first |
55 |
to second node the domain has to lie on the left hand side (in the case of a two dimension surface element |
56 |
the domain has to lie on the left hand side when moving counterclockwise). If the gradient on the |
57 |
surface of the domain is to be calculated rich face elements face to be used. Rich elements on a face |
58 |
are identical to interior elements but with a modified order of nodes such that the 'first' face of the element aligns |
59 |
with the surface of the domian. In \fig{FINLEY FIG 0} |
60 |
elements of the type \finleyelement{Tri3Face} are used. |
61 |
The face element reference number $20$ as a rich face element is defined by the nodes |
62 |
with reference numbers $11$, $0$ and $9$. Notice that the face element $20$ is identical to the |
63 |
interior element $19$ except that, in this case, the order of the node is different to align the first |
64 |
edge of the triangle (which is the edge starting with the first node) with the boundary of the domain. |
65 |
|
66 |
Be aware that face elements and elements in the interior of the domain must match, i.e. a face element must be the face |
67 |
of an interior element or, in case of a rich face element, it must be identical to an interior element. |
68 |
If no face elements are specified |
69 |
\finley implicitly assumes homogeneous natural boundary conditions \index{natural boundary conditions!homogeneous}, |
70 |
i.e. \var{d}=$0$ and \var{y}=$0$, on the entire boundary of the domain. For |
71 |
inhomogeneous natural boundary conditions \index{natural boundary conditions!inhomogeneous}, |
72 |
the boundary must be described by face elements. |
73 |
|
74 |
If discontinuities of the PDE solution are considered contact elements |
75 |
\index{element!contact}\index{contact conditions} are introduced to describe the contact region $\Gamma^{contact}$ |
76 |
even if $d^{contact}$ and $y^{contact}$ are zero. \fig{FINLEY FIG 01} shows a simple example of a mesh |
77 |
of rectangular elements around a contact region $\Gamma^{contact}$ \index{element!contact}. |
78 |
The contact region is described by the |
79 |
elements $4$, $3$ and $6$. Their element type is \finleyelement{Line2_Contact}. |
80 |
The nodes $9$, $12$, $6$, $5$ define contact element $4$, where the coordinates of nodes $12$ and $5$ and |
81 |
nodes $4$ and $6$ are identical with the idea that nodes $12$ and $9$ are located above and |
82 |
nodes $5$ and $6$ below the contact region. |
83 |
Again, the order of the nodes within an element is crucial. There is also the option of using rich elements |
84 |
if the gradient is to be calculated on the contact region. Similarly to the rich face elements |
85 |
these are constructed from two interior elements by reordering the nodes such that |
86 |
the 'first' face of the element above and the 'first' face of the element below the |
87 |
contact regions line up. The rich version of element |
88 |
$4$ is of type \finleyelement{Rec4Face_Contact} and is defined by the nodes $9$, $12$, $16$, $18$, $6$, $5$, $0$ and |
89 |
$2$. |
90 |
|
91 |
\tab{FINLEY TAB 1} shows the interior element types and the corresponding element types to be used |
92 |
on the face and contacts. \fig{FINLEY.FIG:1}, \fig{FINLEY.FIG:2} and \fig{FINLEY.FIG:4} show the ordering of |
93 |
the nodes within an element. |
94 |
|
95 |
\begin{table} |
96 |
\begin{tablev}{l|llll}{textrm}{interior}{face}{rich face}{contact}{rich contact} |
97 |
\linev{\finleyelement{Line2}}{\finleyelement{Point1}}{\finleyelement{Line2Face}}{\finleyelement{Point1_Contact}}{\finleyelement{Line2Face_Contact}} |
98 |
\linev{\finleyelement{Line3}}{\finleyelement{Point1}}{\finleyelement{Line3Face}}{\finleyelement{Point1_Contact}}{\finleyelement{Line3Face_Contact}} |
99 |
\linev{\finleyelement{Tri3}}{\finleyelement{Line2}}{\finleyelement{Tri3Face}}{\finleyelement{Line2_Contact}}{\finleyelement{Tri3Face_Contact}} |
100 |
\linev{\finleyelement{Tri6}}{\finleyelement{Line3}}{\finleyelement{Tri6Face}}{\finleyelement{Line3_Contact}}{\finleyelement{Tri6Face_Contact}} |
101 |
\linev{\finleyelement{Rec4}}{\finleyelement{Line2}}{\finleyelement{Rec4Face}}{\finleyelement{Line2_Contact}}{\finleyelement{Rec4Face_Contact}} |
102 |
\linev{\finleyelement{Rec8}}{\finleyelement{Line3}}{\finleyelement{Rec8Face}}{\finleyelement{Line3_Contact}}{\finleyelement{Rec8Face_Contact}} |
103 |
\linev{\finleyelement{Rec9}}{\finleyelement{Line3}}{\finleyelement{Rec9Face}}{\finleyelement{Line3_Contact}}{\finleyelement{Rec9Face_Contact}} |
104 |
\linev{\finleyelement{Tet4}}{\finleyelement{Tri6}}{\finleyelement{Tet4Face}}{\finleyelement{Tri6_Contact}}{\finleyelement{Tet4Face_Contact}} |
105 |
\linev{\finleyelement{Tet10}}{\finleyelement{Tri9}}{\finleyelement{Tet10Face}}{\finleyelement{Tri9_Contact}}{\finleyelement{Tet10Face_Contact}} |
106 |
\linev{\finleyelement{Hex8}}{\finleyelement{Rec4}}{\finleyelement{Hex8Face}}{\finleyelement{Rec4_Contact}}{\finleyelement{Hex8Face_Contact}} |
107 |
\linev{\finleyelement{Hex20}}{\finleyelement{Rec8}}{\finleyelement{Hex20Face}}{\finleyelement{Rec8_Contact}}{\finleyelement{Hex20Face_Contact}} |
108 |
\end{tablev} |
109 |
\caption{Finley elements and corresponding elements to be used on domain faces and contacts. |
110 |
The rich types have to be used if the gradient of function is to be calculated on faces and contacts, resepctively.} |
111 |
\label{FINLEY TAB 1} |
112 |
\end{table} |
113 |
|
114 |
The native \finley file format is defined as follows. |
115 |
Each node \var{i} has \var{dim} spatial coordinates \var{Node[i]}, a reference number |
116 |
\var{Node_ref[i]}, a degree of freedom \var{Node_DOF[i]} and tag \var{Node_tag[i]}. |
117 |
In most cases \var{Node_DOF[i]}=\var{Node_ref[i]} however, for periodic boundary conditions, |
118 |
\var{Node_DOF[i]} is chosen differently, see example below. The tag can be used to mark nodes sharing |
119 |
the same properties. Element \var{i} is defined by the \var{Element_numNodes} nodes \var{Element_Nodes[i]} |
120 |
which is a list of node reference numbers. The order is crucial. |
121 |
It has a reference number \var{Element_ref[i]} and a tag \var{Element_tag[i]}. The tag |
122 |
can be used to mark elements sharing the same properties. For instance elements above |
123 |
a contact region are marked with $2$ and elements below a contact region are marked with $1$. |
124 |
\var{Element_Type} and \var{Element_Num} give the element type and the number of elements in the mesh. |
125 |
Analogue notations are used for face and contact elements. The following Python script |
126 |
prints the mesh definition in the \finley file format: |
127 |
\begin{python} |
128 |
print "%s\n"%mesh_name |
129 |
# node coordinates: |
130 |
print "%dD-nodes %d\n"%(dim,numNodes) |
131 |
for i in range(numNodes): |
132 |
print "%d %d %d"%(Node_ref[i],Node_DOF[i],Node_tag[i]) |
133 |
for j in range(dim): print " %e"%Node[i][j] |
134 |
print "\n" |
135 |
# interior elements |
136 |
print "%s %d\n"%(Element_Type,Element_Num) |
137 |
for i in range(Element_Num): |
138 |
print "%d %d"%(Element_ref[i],Element_tag[i]) |
139 |
for j in range(Element_numNodes): print " %d"%Element_Nodes[i][j] |
140 |
print "\n" |
141 |
# face elements |
142 |
print "%s %d\n"%(FaceElement_Type,FaceElement_Num) |
143 |
for i in range(FaceElement_Num): |
144 |
print "%d %d"%(FaceElement_ref[i],FaceElement_tag[i]) |
145 |
for j in range(FaceElement_numNodes): print " %d"%FaceElement_Nodes[i][j] |
146 |
print "\n" |
147 |
# contact elements |
148 |
print "%s %d\n"%(ContactElement_Type,ContactElement_Num) |
149 |
for i in range(ContactElement_Num): |
150 |
print "%d %d"%(ContactElement_ref[i],ContactElement_tag[i]) |
151 |
for j in range(ContactElement_numNodes): print " %d"%ContactElement_Nodes[i][j] |
152 |
print "\n" |
153 |
# point sources (not supported yet) |
154 |
write("Point1 0",face_element_typ,numFaceElements) |
155 |
\end{python} |
156 |
|
157 |
The following example of a mesh file defines the mesh shown in \fig{FINLEY FIG 01}: |
158 |
\begin{verbatim} |
159 |
Example 1 |
160 |
2D Nodes 16 |
161 |
0 0 0 0. 0. |
162 |
2 2 0 0.33 0. |
163 |
3 3 0 0.66 0. |
164 |
7 4 0 1. 0. |
165 |
5 5 0 0. 0.5 |
166 |
6 6 0 0.33 0.5 |
167 |
8 8 0 0.66 0.5 |
168 |
10 10 0 1.0 0.5 |
169 |
12 12 0 0. 0.5 |
170 |
9 9 0 0.33 0.5 |
171 |
13 13 0 0.66 0.5 |
172 |
15 15 0 1.0 0.5 |
173 |
16 16 0 0. 1.0 |
174 |
18 18 0 0.33 1.0 |
175 |
19 19 0 0.66 1.0 |
176 |
20 20 0 1.0 1.0 |
177 |
Rec4 6 |
178 |
0 1 0 2 6 5 |
179 |
1 1 2 3 8 6 |
180 |
2 1 3 7 10 8 |
181 |
5 2 12 9 18 16 |
182 |
7 2 13 19 18 9 |
183 |
10 2 20 19 13 15 |
184 |
Line2 0 |
185 |
Line2_Contact 3 |
186 |
4 0 9 12 6 5 |
187 |
3 0 13 9 8 6 |
188 |
6 0 15 13 10 8 |
189 |
Point1 0 |
190 |
\end{verbatim} |
191 |
Notice that the order in which the nodes and elements are given is arbitrary. |
192 |
In the case that rich contact elements are used the contact element section gets |
193 |
the form |
194 |
\begin{verbatim} |
195 |
Rec4Face_Contact 3 |
196 |
4 0 9 12 16 18 6 5 0 2 |
197 |
3 0 13 9 18 19 8 6 2 3 |
198 |
6 0 15 13 19 20 10 8 3 7 |
199 |
\end{verbatim} |
200 |
Periodic boundary condition \index{boundary conditions!periodic} can be introduced by altering \var{Node_DOF}. |
201 |
It allows identification of nodes even if they have different physical locations. For instance, to |
202 |
enforce periodic boundary conditions at the face $x_0=0$ and $x_0=1$ one identifies |
203 |
the degrees of freedom for nodes $0$, $5$, $12$ and $16$ with the degrees of freedom for |
204 |
$7$, $10$, $15$ and $20$, respectively. The node section of the \finley mesh gets now the form: |
205 |
\begin{verbatim} |
206 |
2D Nodes 16 |
207 |
0 0 0 0. 0. |
208 |
2 2 0 0.33 0. |
209 |
3 3 0 0.66 0. |
210 |
7 0 0 1. 0. |
211 |
5 5 0 0. 0.5 |
212 |
6 6 0 0.33 0.5 |
213 |
8 8 0 0.66 0.5 |
214 |
10 5 0 1.0 0.5 |
215 |
12 12 0 0. 0.5 |
216 |
9 9 0 0.33 0.5 |
217 |
13 13 0 0.66 0.5 |
218 |
15 12 0 1.0 0.5 |
219 |
16 16 0 0. 1.0 |
220 |
18 18 0 0.33 1.0 |
221 |
19 19 0 0.66 1.0 |
222 |
20 16 0 1.0 1.0 |
223 |
\end{verbatim} |
224 |
|
225 |
|
226 |
\include{finleyelements} |
227 |
|
228 |
\subsection{Linear Solvers in \LinearPDE} |
229 |
Currently \finley supports the linear solvers \PCG, \GMRES, \PRESTWENTY and \BiCGStab. |
230 |
For \GMRES the options \var{trancation} and \var{restart} of the \method{getSolution} can be |
231 |
used to control the trunction and restart during iteration. Default values are |
232 |
\var{truncation}=5 and \var{restart}=20. |
233 |
The default solver is \BiCGStab but if the symmetry flag is set \PCG is the default solver. |
234 |
\finley supports the solver options \var{iter_max} which specifies the maximum number of iterations steps, |
235 |
\var{verbose}=\True or \False and \var{preconditioner}=\constant{JACOBI} or \constant {ILU0}. |
236 |
In some installations \finley supports the \Direct solver and the |
237 |
solver options \var{reordering}=\constant{util.NO_REORDERING}, |
238 |
\constant{util.MINIMUM_FILL_IN} or \constant{util.NESTED_DISSECTION} (default is \constant{util.NO_REORDERING}), |
239 |
\var{drop_tolerance} specifying the threshold for values to be dropped in the |
240 |
incomplete elimation process (default is 0.01) and \var{drop_storage} specifying the maximum increase |
241 |
in storage allowed in the |
242 |
incomplete elimation process (default is 1.20). |
243 |
|
244 |
\subsection{Functions} |
245 |
\begin{funcdesc}{Mesh}{fileName,integrationOrder=-1} |
246 |
creates a \Domain object form the FEM mesh defined in |
247 |
file \var{fileName}. The file must be given the \finley file format. |
248 |
If \var{integrationOrder} is positive, a numerical integration scheme |
249 |
chosen which is accurate on each element up to a polynomial of |
250 |
degree \var{integrationOrder} \index{integration order}. Otherwise |
251 |
an appropriate integration order is chosen independently. |
252 |
\end{funcdesc} |
253 |
|
254 |
\begin{funcdesc}{Interval}{n0,order=1,l0=1.,integrationOrder=-1, \\ |
255 |
periodic0=\False,useElementsOnFace=\False} |
256 |
Generates a \Domain object representing a interval $[0,l0]$. The interval is filled with |
257 |
\var{n0} elements. |
258 |
For \var{order}=1 and \var{order}=2 |
259 |
\finleyelement{Line2} and |
260 |
\finleyelement{Line3} are used, respectively. |
261 |
In the case of \var{useElementsOnFace}=\False, |
262 |
\finleyelement{Point1} are used to describe the boundary points. |
263 |
In the case of \var{useElementsOnFace}=\True (this option should be used if gradients |
264 |
are calculated on domain faces), |
265 |
\finleyelement{Line2} and |
266 |
\finleyelement{Line3} are used on both ends of the interval. |
267 |
If \var{integrationOrder} is positive, a numerical integration scheme |
268 |
chosen which is accurate on each element up to a polynomial of |
269 |
degree \var{integrationOrder} \index{integration order}. Otherwise |
270 |
an appropriate integration order is chosen independently. If |
271 |
\var{periodic0}=\True, periodic boundary conditions \index{periodic boundary conditions} |
272 |
along the $x_0$-directions are enforced. That means when for any solution of a PDE solved by \finley |
273 |
the value at $x_0=0$ will be identical to the values at $x_0=\var{l0}$. |
274 |
\end{funcdesc} |
275 |
|
276 |
\begin{funcdesc}{Rectangle}{n0,n1,order=1,l0=1.,l1=1., integrationOrder=-1, \\ |
277 |
periodic0=\False,periodic1=\False,useElementsOnFace=\False} |
278 |
Generates a \Domain object representing a two dimensional rectangle between |
279 |
$(0,0)$ and $(l0,l1)$ with orthogonal edges. The rectangle is filled with |
280 |
\var{n0} elements along the $x_0$-axis and |
281 |
\var{n1} elements along the $x_1$-axis. |
282 |
For \var{order}=1 and \var{order}=2 |
283 |
\finleyelement{Rec4} and |
284 |
\finleyelement{Rec8} are used, respectively. |
285 |
In the case of \var{useElementsOnFace}=\False, |
286 |
\finleyelement{Line2} and |
287 |
\finleyelement{Line3} are used to subdivide the edges of the rectangle, respectively. |
288 |
In the case of \var{useElementsOnFace}=\True (this option should be used if gradients |
289 |
are calculated on domain faces), |
290 |
\finleyelement{Rec4Face} and |
291 |
\finleyelement{Rec8Face} are used on the edges, respectively. |
292 |
If \var{integrationOrder} is positive, a numerical integration scheme |
293 |
chosen which is accurate on each element up to a polynomial of |
294 |
degree \var{integrationOrder} \index{integration order}. Otherwise |
295 |
an appropriate integration order is chosen independently. If |
296 |
\var{periodic0}=\True, periodic boundary conditions \index{periodic boundary conditions} |
297 |
along the $x_0$-directions are enforced. That means when for any solution of a PDE solved by \finley |
298 |
the value on the line $x_0=0$ will be identical to the values on $x_0=\var{l0}$. |
299 |
Correspondingly, |
300 |
\var{periodic1}=\False sets periodic boundary conditions |
301 |
in $x_1$-direction. |
302 |
\end{funcdesc} |
303 |
|
304 |
\begin{funcdesc}{Brick}{n0,n1,n2,order=1,l0=1.,l1=1.,l2=1., integrationOrder=-1, \\ |
305 |
periodic0=\False,periodic1=\False,periodic2=\False,useElementsOnFace=\False} |
306 |
Generates a \Domain object representing a three dimensional brick between |
307 |
$(0,0,0)$ and $(l0,l1,l2)$ with orthogonal faces. The brick is filled with |
308 |
\var{n0} elements along the $x_0$-axis, |
309 |
\var{n1} elements along the $x_1$-axis and |
310 |
\var{n2} elements along the $x_2$-axis. |
311 |
For \var{order}=1 and \var{order}=2 |
312 |
\finleyelement{Hex8} and |
313 |
\finleyelement{Hex20} are used, respectively. |
314 |
In the case of \var{useElementsOnFace}=\False, |
315 |
\finleyelement{Rec4} and |
316 |
\finleyelement{Rec8} are used to subdivide the faces of the brick, respectively. |
317 |
In the case of \var{useElementsOnFace}=\True (this option should be used if gradients |
318 |
are calculated on domain faces), |
319 |
\finleyelement{Hex8Face} and |
320 |
\finleyelement{Hex20Face} are used on the brick faces, respectively. |
321 |
If \var{integrationOrder} is positive, a numerical integration scheme |
322 |
chosen which is accurate on each element up to a polynomial of |
323 |
degree \var{integrationOrder} \index{integration order}. Otherwise |
324 |
an appropriate integration order is chosen independently. If |
325 |
\var{periodic0}=\True, periodic boundary conditions \index{periodic boundary conditions} |
326 |
along the $x_0$-directions are enforced. That means when for any solution of a PDE solved by \finley |
327 |
the value on the plane $x_0=0$ will be identical to the values on $x_0=\var{l0}$. Correspondingly, |
328 |
\var{periodic1}=\False and \var{periodic2}=\False sets periodic boundary conditions |
329 |
in $x_1$-direction and $x_2$-direction, respectively. |
330 |
\end{funcdesc} |
331 |
|
332 |
\begin{funcdesc}{GlueFaces}{meshList,safetyFactor=0.2,tolerance=1.e-13} |
333 |
Generates a new \Domain object from the list \var{mehList} of \finley meshes. |
334 |
Nodes in face elements whose difference of coordinates is less then \var{tolerance} times the |
335 |
diameter of the domain are merged. The corresponding face elements are removed from the mesh. |
336 |
|
337 |
TODO: explain \var{safetyFactor} and show an example. |
338 |
\end{funcdesc} |
339 |
|
340 |
\begin{funcdesc}{JoinFaces}{meshList,safetyFactor=0.2,tolerance=1.e-13} |
341 |
Generates a new \Domain object from the list \var{mehList} of \finley meshes. |
342 |
Face elements whose nodes coordinates have difference is less then \var{tolerance} times the |
343 |
diameter of the domain are combined to form a contact element \index{element!contact} |
344 |
The corresponding face elements are removed from the mesh. |
345 |
|
346 |
TODO: explain \var{safetyFactor} and show an example. |
347 |
\end{funcdesc} |