Meshes

To understand the usage of esys.finley one needs to have an understanding of how the finite element meshes are defined. Figure 7.1 shows an example of the subdivision of an ellipse into so called elements . In this case, triangles have been used but other forms of subdivisions can be constructed, e.g. into quadrilaterals or, in the three dimensional case, into tetrahedrons and hexahedrons. The idea of the finite element method is to approximate the solution by a function which is a polynomial of a certain order and is continuous across it boundary to neighbor elements. In the example of Figure 7.1 a linear polynomial is used on each triangle. As one can see, the triangulation is quite a poor approximation of the ellipse. It can be improved by introducing a midpoint on each element edge then positioning those nodes located on an edge expected to describe the boundary, onto the boundary. In this case the triangle gets a curved edge which requires a parametrization of the triangle using a quadratic polynomial. For this case, the solution is also approximated by a piecewise quadratic polynomial (which explains the name isoparametrical elements), see Reference [28,5] for more details.

The union of all elements defines the domain of the PDE. Each element is defined by the nodes used to describe its shape. In Figure 7.1 the element, which has type Tri3, with element reference number $ 19$ is defined by the nodes with reference numbers $ 9$, $ 11$ and 0 . Notice that the order is counterclockwise. The coefficients of the PDE are evaluated at integration nodes with each individual element. For quadrilateral elements a Gauss quadrature scheme is used. In the case of triangular elements a modified form is applied. The boundary of the domain is also subdivided into elements. In Figure 7.1 line elements with two nodes are used. The elements are also defined by their describing nodes, e.g. the face element reference number $ 20$ which has type Line2 is defined by the nodes with the reference numbers $ 11$ and 0. Again the order is crucial, if moving from the first to second node the domain has to lie on the left hand side (in the case of a two dimension surface element the domain has to lie on the left hand side when moving counterclockwise). If the gradient on the surface of the domain is to be calculated rich face elements face to be used. Rich elements on a face are identical to interior elements but with a modified order of nodes such that the 'first' face of the element aligns with the surface of the domain. In Figure 7.1 elements of the type Tri3Face are used. The face element reference number $ 20$ as a rich face element is defined by the nodes with reference numbers $ 11$, 0 and $ 9$. Notice that the face element $ 20$ is identical to the interior element $ 19$ except that, in this case, the order of the node is different to align the first edge of the triangle (which is the edge starting with the first node) with the boundary of the domain.

Be aware that face elements and elements in the interior of the domain must match, i.e. a face element must be the face of an interior element or, in case of a rich face element, it must be identical to an interior element. If no face elements are specified esys.finley implicitly assumes homogeneous natural boundary conditions , i.e. d=0 and y=0, on the entire boundary of the domain. For inhomogeneous natural boundary conditions , the boundary must be described by face elements.

If discontinuities of the PDE solution are considered contact elements are introduced to describe the contact region $ \Gamma^{contact}$ even if $ d^{contact}$ and $ y^{contact}$ are zero. Figure 7.2 shows a simple example of a mesh of rectangular elements around a contact region $ \Gamma^{contact}$ . The contact region is described by the elements $ 4$, $ 3$ and $ 6$. Their element type is Line2_Contact. The nodes $ 9$, $ 12$, $ 6$, $ 5$ define contact element $ 4$, where the coordinates of nodes $ 12$ and $ 5$ and nodes $ 4$ and $ 6$ are identical with the idea that nodes $ 12$ and $ 9$ are located above and nodes $ 5$ and $ 6$ below the contact region. Again, the order of the nodes within an element is crucial. There is also the option of using rich elements if the gradient is to be calculated on the contact region. Similarly to the rich face elements these are constructed from two interior elements by reordering the nodes such that the 'first' face of the element above and the 'first' face of the element below the contact regions line up. The rich version of element $ 4$ is of type Rec4Face_Contact and is defined by the nodes $ 9$, $ 12$, $ 16$, $ 18$, $ 6$, $ 5$, 0 and $ 2$.

Table 7.1 shows the interior element types and the corresponding element types to be used on the face and contacts. Figure 7.3, Figure 7.4 and Figure 7.5 show the ordering of the nodes within an element.


Table 7.1: Finley elements and corresponding elements to be used on domain faces and contacts. The rich types have to be used if the gradient of function is to be calculated on faces and contacts, respectively.
\begin{table}\begin{tablev}{l\vert llll}{textrm}{interior}{face}{rich face}{cont...
...\it Hex20Face_Contact}\index{finley!Hex20Face_Contact}}
\end{tablev}
\end{table}


The native esys.finley file format is defined as follows. Each node i has dim spatial coordinates Node[i], a reference number Node_ref[i], a degree of freedom Node_DOF[i] and tag Node_tag[i]. In most cases Node_DOF[i]=Node_ref[i] however, for periodic boundary conditions, Node_DOF[i] is chosen differently, see example below. The tag can be used to mark nodes sharing the same properties. Element i is defined by the Element_numNodes nodes Element_Nodes[i] which is a list of node reference numbers. The order is crucial. It has a reference number Element_ref[i] and a tag Element_tag[i]. The tag can be used to mark elements sharing the same properties. For instance elements above a contact region are marked with $ 2$ and elements below a contact region are marked with $ 1$. Element_Type and Element_Num give the element type and the number of elements in the mesh. Analogue notations are used for face and contact elements. The following Python script prints the mesh definition in the esys.finley file format:
\begin{python}
print ''%s\n''%mesh_name
...

The following example of a mesh file defines the mesh shown in Figure 7.2:

Example 1
2D Nodes 16
0   0 0 0.   0.
2   2 0 0.33 0.
3   3 0 0.66 0.
7   4 0 1.   0.
5   5 0 0.   0.5
6   6 0 0.33 0.5
8   8 0 0.66 0.5
10 10 0 1.0  0.5
12 12 0 0.   0.5
9   9 0 0.33 0.5
13 13 0 0.66 0.5
15 15 0 1.0  0.5
16 16 0 0.   1.0
18 18 0 0.33 1.0
19 19 0 0.66 1.0
20 20 0 1.0  1.0
Rec4 6
 0 1  0  2  6  5
 1 1  2  3  8  6
 2 1  3  7 10  8
 5 2 12  9 18 16
 7 2 13 19 18  9
10 2 20 19 13 15
Line2 0
Line2_Contact 3
 4 0  9 12  6 5
 3 0 13  9  8 6
 6 0 15 13 10 8
Point1 0
Notice that the order in which the nodes and elements are given is arbitrary. In the case that rich contact elements are used the contact element section gets the form
Rec4Face_Contact 3
 4 0  9 12 16 18  6  5  0  2
 3 0 13  9 18 19  8  6  2  3
 6 0 15 13 19 20 10  8  3  7
Periodic boundary condition can be introduced by altering Node_DOF. It allows identification of nodes even if they have different physical locations. For instance, to enforce periodic boundary conditions at the face $ x_0=0$ and $ x_0=1$ one identifies the degrees of freedom for nodes 0, $ 5$, $ 12$ and $ 16$ with the degrees of freedom for $ 7$, $ 10$, $ 15$ and $ 20$, respectively. The node section of the esys.finley mesh gets now the form:
2D Nodes 16
0   0 0 0.   0.
2   2 0 0.33 0.
3   3 0 0.66 0.
7   0 0 1.   0.
5   5 0 0.   0.5
6   6 0 0.33 0.5
8   8 0 0.66 0.5
10  5 0 1.0  0.5
12 12 0 0.   0.5
9   9 0 0.33 0.5
13 13 0 0.66 0.5
15 12 0 1.0  0.5
16 16 0 0.   1.0
18 18 0 0.33 1.0
19 19 0 0.66 1.0
20 16 0 1.0  1.0

(0,0) (0,0)2 (1,0)(1,0)28 (30,0)2

(0,0) (0,0)2 (1,0)(1,0)28 (30,0)2

(0,30) (0,0)2 (0,1)(0,1)28 (0,30)2

(0,30) (0,0)2 (0,1)(0,1)28 (0,30)2

(0,30) (0,0)(-1,1)30 (0,0)2 (-30,30)2

(0,15) (0,0)(-4,3)20 (0,0)2 (-20,15)2

(0,-15) (0,0)(-4,-3)20 (0,0)2 (-20,-15)2

(0,-15) (-0.7,-0.7)(-4,-3)18.7 (0,0)2 (-20,-15)2

(0,15) (0,0)(-2,3)10 (0,0)2 (-10,15)2

(0,15) (0,0)(-2,3)9.4 (0,0)2 (-10,15)2

(0,0) (0,0)2 (10,0)2 (1,0)(1,0)28 (20,0)2 (30,0)2

(0,0) (0,0)2 (10,0)2 (1,0)(10,0)3(1,0)8 (20,0)2 (30,0)2

(0,30) (0,0)2 (0,10)2 (0,1)(0,1)28 (0,20)2 (0,30)2

(0,30) (0,0)2 (0,10)2 (0,1)(0,10)3(0,1)8 (0,20)2 (0,30)2

(0,30) (0,0)(-1,1)30 (0,0)2 (-10,10)2 (-20,20)2 (-30,30)2

(0,15) (0,0)(-4,3)20 (0,0)2 (-6.66,5)2 (-13.33,10)2 (-20,15)2

(0,-15) (0,0)(-4,-3)20 (0,0)2 (-6.66,-5)2 (-13.33,-10)2 (-20,-15)2

(0,-15) (-0.7,-0.7)(-6.66,-5)3(-4,-3)5.1 (0,0)2 (-6.66,-5)2 (-13.33,-10)2 (-20,-15)2

(0,15) (0,0)(-2,3)10 (0,0)2 (-3.33,5)2 (-6.66,10)2 (-10,15)2

(0,15) (-0.6,0.8)(-3.33,5)3(-2,3)2.35 (0,0)2 (-3.33,5)2 (-6.66,10)2 (-10,15)2

(0,0) (0,0)2 (15,0)2 (1,0)(1,0)28 (30,0)2

(0,0) (0,0)2 (15,0)2 (1,0)(15,0)2(1,0)13 (30,0)2

(0,30) (0,0)2 (0,15)2 (0,1)(0,1)28 (0,30)2

(0,30) (0,0)2 (0,15)2 (0,1)(0,15)2(0,1)13 (0,30)2

(0,30) (0,0)(-1,1)30 (0,0)2 (-15,15)2 (-30,30)2

(0,15) (0,0)(-4,3)20 (0,0)2 (-10,7.5)2 (-20,15)2

(0,-15) (0,0)(-4,-3)20 (0,0)2 (-10,-7.5)2 (-20,-15)2

(0,-15) (-0.7,-0.7)(-10,-7.5)2(-4,-3)8.4 (0,0)2 (-10,-7.5)2 (-20,-15)2

(0,15) (0,0)(-2,3)10 (0,0)2 (-5,7.5)2 (-10,15)2

(0,15) (-0.6,0.8)(-5,7.5)2(-2,3)3.9 (0,0)2 (-5,7.5)2 (-10,15)2

Figure 7.3: Elements of order 1
\begin{figure}\par
\begin{picture}(170,210)
\thicklines
\par
\put(20,205){\cir...
...{{\it 4}}
\put(143,64){{\it 7}}
\put(143,34){{\it 3}}
\end{picture}
\end{figure}

Figure 7.4: Elements of order 2
\begin{figure}
\setlength{\unitlength}{1mm}
\par
\begin{picture}(170,210)
\thi...
...
\put(94.5,56.5){{\it 20}}
\put(132.5,56.5){{\it 18}}
\end{picture}
\end{figure}

Figure 7.5: Additional shape functions
\begin{figure}\begin{picture}(170,210)
\thicklines
\par
\put(50,095){\usebox{\...
...{{\it 6}}
\put(83,094){{\it 2}}
\put(64,105){{\it 9}}
\end{picture}
\end{figure}


Table: Solvers available for esys.finley and the PASO package and the relevant options in SolverOptions. MKL supports SolverOptions.MINIMUMFILLIN and SolverOptions.NESTEDDISSECTION reordering. Currently the UMFPACK interface does not support any reordering.
setSolverMethod DIRECT PCG GMRES TFQMR MINRES PRES20 BICGSTAB LUMPING
setReordering $ \checkmark$
setRestart $ \checkmark$ $ 20$
setTruncation $ \checkmark$ $ 5$
setIterMax $ \checkmark$ $ \checkmark$ $ \checkmark$ $ \checkmark$ $ \checkmark$ $ \checkmark$
setTolerance $ \checkmark$ $ \checkmark$ $ \checkmark$ $ \checkmark$ $ \checkmark$ $ \checkmark$
setAbsoluteTolerance $ \checkmark$ $ \checkmark$ $ \checkmark$ $ \checkmark$ $ \checkmark$ $ \checkmark$
setReordering $ \checkmark$



Table: Preconditioners available for esys.finley and the PASO package and the relevant options in SolverOptions.
setPreconditioner NO_PRECONDITIONER AMG JACOBI GAUSS_SEIDEL REC_ILU RILU ILU0 DIRECT
status: later later $ \checkmark$ $ \checkmark$ $ \checkmark$ later $ \checkmark$ later
setCoarsening $ \checkmark$
setLevelMax $ \checkmark$
setCoarseningThreshold $ \checkmark$
setMinCoarseMatrixSize $ \checkmark$
setNumSweeps $ \checkmark$ $ \checkmark$
setNumPreSweeps $ \checkmark$
setNumPostSweeps $ \checkmark$
setInnerTolerance
setDropTolerance
setDropStorage
setRelaxationFactor $ \checkmark$
adaptInnerTolerance
setInnerIterMax




Subsections
esys@esscc.uq.edu.au