/[escript]/trunk/doc/user/finley.tex
ViewVC logotype

Annotation of /trunk/doc/user/finley.tex

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26