- Email: [email protected]
Multigrid microsoftware for engineers on a Macintosh computer
Advances in Engineering Software 14 (1992) 219-233
Multigrid microsoftware for engineers on a Macintosh computer C. Le Carlier de Veslud, G. Maurice & R. Kouitat Laboratoire de Science et G~nie des Surfaces, Ecole des Mines, Parc de Saurupt, 54042 Nancy Cedex, France
We present the 2-D version of a multigrid finite element software to solve engineering problems on microcomputers. Some 3-D developments are also presented. This software uses several automatic mesh generators starting from the nodes on the boundaries. It is possible to use a logarithmic scale to determine the nodes in the neighbouring rapid variation of the unknown. Our meshing technique allows to automatically modify the shapes of the irregular elements using the centre of gravity of his neighbours. A 3-D mesh generator has been also developed. It is provided with an automatic process to refine a network in the desired region and several levels of zoom. It is also self adaptable to the problem because it can give maximum information in a precise region. It works with linear and quadratic triangles and rectangles in the 2D version and with linear parallelepipeds in the 3D version. The shapes are visualised in space with hidden surfaces removed. The user data input routine with its thorough use of menus and the mouse makes this program very user friendly. At the moment, this program can only be used for a limited set of physical situations, since it only treats problems with plane elasticity and with axisymmetric deformations and constraints. The program is being extended to include problems in nonlinear elasticity as well as fluid flows with yield stress. Key words: multigrid method, finite elements, meshing.
1 INTRODUCTION Among the methods for solving partial differential equations, the finite elements method stands out with its performance and adaptability, which tends to place it as the favoured general method for solving linear systems of considerable size, often ill conditioned. The classical iterative method for solving these systems often leads complex calculations, and since the convergence rate is close to unity, a considerable computing time is required before a solution is obtained. In the 60s appeared the multigrid methods 9 and in the 80s they were developed in this context.~ It consists of a process coupled with an iterative numerical method which considerably accelerates the convergence. Their experimentation on linear, as well as nonlinear classical problems, has rapidly lead to results which may be optimal. 2 These methods are based on a sequence of numerical processes which must be fitted judiciously in order to give rise to a good package. Then the solution to each problem is a particular case. Advances in Engineering Software 0965-9978/92/$05.00 © 1992 Elsevier Science Publishers Ltd.
219
The finite difference method, simple and symmetric in space, has become the first field of application for most multigrid packages. In the case of the finite element method, though it has a greater adaptability, the space symmetry is last. 3-5'8 So the finite element method is not well adapted to the multigrid method and is used nowadays with simple or particular elements. Its implementation on microcomputer leads to the solution of difficult problems of memory and conviviality. Therefore, it is important to solve the question of data structuring and that of the creation of algorithms which allow the implementation of the method for various elements, taking in consideration a small memory and a low speed of execution. In this paper, we describe the finite element multigrid software that we have developed on a Macintosh microcomputer. We also try to put in light our original contribution. The study is divided into five parts. In the first we recall the basics of multigrid methods. In the second part we describe the main part of the program, the data structure, the different methods we have used and the main problems we have encountered and solved in the combination of the multigrid method with the finite elements method. In the following two sections we present the
220
C. L e Carlier de Veslud, G. Maurice, R. Kouitat
heart of our work: the two mesh generators called ISOPAR and D E L A U N A Y . In the last part we present the solution of some applications with this software.
2 A VIEW ON MULTIGRID METHODS
The use of the finite differences or the finite elements methods for obtaining solutions to certain problems, sometimes gives rise to big linear systems of type: 4'7 Lu =
f
where L is the discrete operator of the problem, u is the unknown vector, f is a known vector. The multigrid method is an iterative process for finding the solutions of the discrete problems. It uses a sequence of grids more and more coarse. This method requires the exact solution of the linear system on the coarsest grid, obtained by a direct method (as opposed to an iterative method which yields only approximate solutions). 1° The multigrid method is efficient for big linear systems. Depending on the size of the system, the computation of the numerical solution of the discrete problem L. = f, can be obtained using two kinds of methods. Among the most common direct methods, one can mention the Gauss methods, the Choleski and the conjuguate gradient methods. However, many disadvantages appear, especially for large systems: a high sensibility to compute time errors leading to erroneous solutions; a large occupation of memory; a very high computing time and number of operations. Then generally iterative methods are preferred. These methods give an approximation to the solution which is hoped to be better as the iterative processes. They use a fixed point model with the aid of the operator splitting: L = M - N. Therefore, one is lead to matrice iterations of the type: Mun+l
=
Southwell tried to develop methods for solving discrete problems using a coarse discretization." In the early 60s, Fedorenko combined these two aspects in a method called 'multigrid' for the solution of Poisson's equation, r2 Then Bakhavalor undertook the study of problems with variable coefficients, j3 Next, Astrachanov proposed a multigrid formulation for finite element problems and he showed the low cost of the methodJ 4 From the 70s, Brandt devised general theories which include earlier works and he proposed efficient solution methods for certain problems, essentially those of fluid mechanics. ~5 From the end of the 70s considerable work has been done on multigrids by Hackbush, ~': Stfiben & Trotenberg, J6 as well as by Maitre in France. ~7 Principle o f the method
In order to describe the main properties of the method, let us use the same one-dimensional example as hackbush. 2 d2u dx: + f
=
x ~ ]0, 1[, u(o)
0
=
u(1)
=
0
The discretization of the problem in N steps with size h and the use of boundary conditions lead to a system of the form: Lh llh :
f t,
First stage."
To fix our ideas, let us solve the linear system by the relaxed Jacobi method (~o = 1/2). This method leads to a slow decrease of the global numerical error which is a discrete signal involving high and low frequencies. Both types of frequencies have different behaviour during iterations. It is shown mathematically that high frequencies present a convergence rate of about 0. Therefore, their effect is eliminated after some iterations. During this stage, where smoothing takes place, the error has not decreased significantly but has taken a regular aspect.
Nun + f
where u, is the solution at the nth iteration of the method and f is the second member of the initial equation. M is chosen so that the successive linear system are easily solved. The most known are those of Jacobi and GaussSeidel. Although a judicious choice of M and N gives rise to mathematical convergence, numerically the convergence is quite slow and divergence may occur. As far back as 1940, the Fourier analysis has contributed to the discovery of one of the essential properties of iterative methods: the rapid elimination of high frequencies of Un during iterations, that is the property known as smoothing, jj Elsewhere, Schr6der with his reduction method and
Second stage:
Let fih be an approximation ofuh obtained after v smoothing iterations, and, eh =fih - uh. One obtains Lheh = Lhfih-- Lhuh = Lhfih -- fh = dh. Since fih and fh are known, so is d h and then eh is the solution of the system Lheh = dh which is similar to the base system. Let consider a much coarser discretization with steps of size 2h (Fig. 1). 0 []
"~1"--2h--"'~4 h ~ o
[]
o
[] Fig. 1.
o
1 []
o
[]
Multigrid
m i c r o s o f t w a r e f o r e n g i n e e r s on a M a c i n t o s h
starting grid
starting grid
coarse grid
intermediate grid
O smoothing '5~Restriction, 7( Prolongation,[ ] exactresolution
coarse grid
Fig. 2.
Fig.
On this grid, the error is a solution of L2he2h = d2h where L2h is the same form a s L h and dEh is the restriction of dh on the coarse grid. Since the number of unknowns is small compared with those of the former grid one can compute exactly e2h = L~Jd2h and then by a p interpolation recover a new approximation of e h. So fl h : I1h ~ pe2h constitutes a new approximation of Uh which is much more precise than fih since the initial error has been corrected. This second stage is called the correction on a coarse grid. This cycle based on the use of two grids is the simplest multigrid method. It involves a smoothing followed by correction on a coarse grid (restriction on the coarse grid, exact calculation of the e r r o r LEh : L2h~d2h , extension of e2h on the sharp grid, calculation of fib). A useful way to represent this cycle is sketched in Fig. 2.
221
computer
4.
approximation fil = Pfit-~ is obtained. The multigrid process can be applied at any step for more efficiency. The full multigrid process has three main advantages: (1) l is not fixed a p r i o r i . This makes the mesh refinement a dynamic process. (2) A better accuracy than that of the discretization error can be obtained by using two distinct operators Lh and lh for the based iterations and for the calculus of residues (a method with defect correction). (3) Very good results are obtained for nonlinear problems for which the iterative process requires a good initial estimate of the solution. Therefore, the number of operations tends to be optimal. 3 THE MAIN PART OF THE PROGRAM: DIVQUAD
Multigrid method Frequently, in large size finite element systems the convergence acceleration obtained from the bigrid cycle is insufficient and the size of the linear systems on the coarse grid is too big to be solved by a direct method. Therefore, one combines smoothing and restriction until a sufficiently coarse grid is obtained, which allows an exact calculation of the error. From this point one goes back to the original sharp grid by multiple extensions. Figure 3 represents a cycle with 3 grids. Two or more corrections can be done at each level in order t. As for every iterative method the efficacy of the multigrid version the method constructs a good initial approximation in the following way. Consider a multigrid process with ! + 1 grids numbered from 0 to 1 where 0 stands for the most coarse grid; on which the exact calculation is performed; and 1 for the sharpest grid. At the zero level u0 = Lo ~f0 can be obtained exactly. A coarse approximation fi~ = pu0 on grid 1 is obtained by interpolation of u0 on it. The improvement of this solution by the use of a multigrid process leads to the vector fi~. Then u2 = PUt is calculated and what has been done for u~ is repeated and so on. At level ! a good initial
starting grid intermediategrid coarse grid
Fig.
3.
The multigrid methods as well as their implementation on computers, which requires a more particular input data structure than that of finite elements, still stand as research work. The difficulty is increased by the wish for generality of the program. DIVQUAD; written in F O R T R A N 77 language and implemented in an Apple Macintosh computer partly solves this problem. For the moment it treats plane elasticity problems (plane stresses or deformations and axisymmetric). Input data structure The multigrid method is widely recursive. The full multigrid version offers the advantage to vary the size which is unknown a p r i o r i . Therefore, it gives a solution to the problem of memory saving. From a strictly geometrical point of the view, we note that the method starts on the most coarse grid and refines step after step. The exact geometric definition of the domain is therefore necessary from the start especially concerning curved boundaries. Unknown and second member The unknown vector U and the second member B are dimensioned as large size 'tables'. After, they are cut into parts relative to each grid as one does along the progression of grids as shown in Fig. 5 below. With regard to the unknown vector U, one goes from grid i + 1 and uses the space of grid i for the solution of the correction equation in the descending phase (restriction). In the ascending phase the space of the grid i is used
222
C. Le Carlier de Veslud, G. Maurice, R. Kouitat gridi +
gridi + I ++
networkof the grid 1 +
networkof the grid 2
VectorU NTI elements NT2 elements
Fig. 5.
Fig. 7.
for correction by interpolation on grid i + 1. For the vector B, the space of the grid i + 1 is intended for the second member while that of grid i is for the default and
ensure the existence of a knot of grids wherever the force applies). This file is fixed and is common to all grids. There is an evolution of the situation during computation in the case of the two other types of forces. Effectively, the notion of loaded elements or loaded edges varies with the grids since these entities are subdivided. Then we have two files which vary as one goes along the display of grids. At the most coarse level they contain only the bare minimum for the description of the loading. Later the loading on the second grid is added, then that of the third and so on. As before, a set of addresses allows the reading of the loading which corresponds to the grid as shown in Fig. 8. It is to be noted that the most coarse grid must be sufficiently accurate in order to represent the loading of the problem under study.
SO o n .
When a new grid i + 1 is constructed geometrically the vector Ui+t is deduced from U~ by interpolation and Bi+ ~ is constructed according to 'opposing forces'. Therefore, whenever a new grid is constructed, the main parameters (numbers of points, number of elements, etc.) are memorised (saved). Therefore, the number of unknowns on all grids in use is known. An elementary procedure gives at each level the starting and ending addresses of the grid space. Mesh generation
In this case as well there is the problem of how to differentiate the grids. With regard to points the solution is simple: when a new grid is constructed the only points that are saved are those which appear for the first time making the table of points into the structure of the sketch on Fig. 6 below. The first N t points are in grid 1 and the first Nz points in grid 2. The points having numbers between NI and N2 have appeared recently and are in grid 2 but not in grid 1. Such an artifice cannot be applied for mesh generations which are saved in consecutive files as can be seen in Fig. 7. As previously mentioned, an archivage procedure allows the recovery of starting and ending addresses of points and elements at any level. Applied forces
In the case of a plane problem applied forces are of three types: punctual forces, lineal forces and surface forces. With regards to the first forces there are no particular difficulties and they are saved in a file having the size of the most coarse grid (however care must be taken to
grid n°2j N2 points )
Boundary conditions
Two types of boundary conditions must be distinguished: (1) Geometrical boundary conditions that are tied to a point and which allow the placing of points on the boundary and during displays of grids to create new points in conformity with boundaries (curved boundaries for example); (2) Numerical boundary conditions which express the imposed displacements at some points (lockings for example). However, these boundaries also arise during the display of a grid. For example if a new point is created on the boundary between two locked points, it is also locked unless there is an explicit different order. Geometrical boundary conditions are defined through an edge file describing the external contour. The edges defined by the mesh generators and are divided into two types: straight edges and curved edges. Each point is File
Loadedelement
<;>
t
gridn° 1 Loading
4 Loadedelement s
Y
Loading gridn° 2 !
Loading grid n° 3
grin n° l ( N1 points ) Fig. 6.
Fig. 8.
Multigrid microsoftware for engineers on a Macintosh computer IFR NO2X 1"
NO2Y UX UY 1" t
I geometric boundary code
223
I eventaalyimposeddisplacements
Fig. 9.
assigned a boundary code related if needed to the edges. Internal points have a null boundary code. Numerical boundary conditions are stored in memory by two locking codes (one on the x and the other on the y) which are null if there is no imposed displacement and have the value 1 otherwise. This information is gathered in a file having the same size as that of the points and varying according to the grids. Each point involves the locking codes in x and y according to Fig. 9.
Stocking of the stiffness matrix
We have used a M O R S E type stocking. The topologically no null elements of the matrix are compacted line by line then put end to end in a one-dimensional table. A table o f pointers gives the location of the beginning and the end of each line. A second table gives the null location of the compacted term for each line. It is a particularly efficient method for large sparse matrices since it allows a rapid pass through the file which is well adapted to the solution method (conjugate gradient for example). Elsewhere the two tables of 'pointers' give precious indications on the mesh generation (very useful for example in the reduction of the band width). The materials and the stresses
The materials are defined again each time the program is used. It is possible to use many materials the features of which are stocked in sequential files. The stresses are put into files in such a way that their screen drawing is possible later. Grid divisions
The basic idea of the software is to closely relate the solution and the mesh generation. For this reason we use the 'full multigrid' process described previously. The solution of the problem on the coarse grid permits the calculation of stresses at Gauss points. These are projected on the edges by a least square method. Since now, the variations of stresses on each element compared with the mean variations on the whole domain give us precious information on the mesh generation quality. The elements that are considered to be to stressed are then divided into four or nine and the others are adjusted such that a link is made with the total mesh generation. This process leads to a new data base according to the new elements and features so created (components,
/\ undivided startingmeshing
divided meshing
Fig. 10.
boundaries, limits, etc.). With regard to linear elements, the created points must eventually be projected on curved boundaries (necessity of boundary codes). The quadratic elements already involve the notion of curvature under their isoparametric form so that the projection is not necessary (Fig. 11). 8 After having verified that each new point does not already exist, the new elements are then always constructed in the same order so as to make location easier. At this level, files of linear and surface forces are made complete in order to take into account new edges and new elements created in the new grid. This new grid is stocked according to the data structure described previously. A similar method is used during the ascending of grids in the multigrid process. The corrected error on the coarse grid is interpolated on the sharp grid by using proper shape functions of each element. Then, one obtains a correction scheme which can be used for any element. Numerical smoothers and solvers
The smoothing is actually done by a relaxed GaussSeidel method with 5 or 10 iterations. A smoothing by a relaxed Jacobi method is also possible but it is less efficient. However, we must note that this last method allows smoothing without the necessity of stocking the stiffness matrix which is computed again after each iteration. The solution at the most coarse level, is obtained by using one of the most efficient numerical methods which is a preconditioned conjugate gradient method with elementsto dividein 4
\\ initial meshing
\\
final meshing afterdivis ~n in 4 and connectionto the rest of the meshing Fig. 11.
224
C. Le Carlier de Veslud, G. Maurice, R. Kouitat Re~ling of the coarse [ meshing - Resolution
coarse grkl []
of the grid [ in~rpolation ] coarsegrid . . . .
refinement
Smoothing [ smoothgrkl -- - ~ cornputm~o!t~ defect resign
coarse
grid . . . . . . . . . . . .
Resolutionof the correction srr~o~ grid '-----~----~,i I equation-Coxrec~n [ co~se grid
l
The precisionof the ][ computin~is enough?
~
Yes
l( -~y~b~s [] Exit ~ol~i~ 0
Graphicoutputoftheresults[
Smoot~ ~xnl:~n
Fig. 12.
incomplete Choleski factorization (ICCG). Nevertheless, this method is not always stable for elements of degree 2. Therefore, the solution is done by a conjugate gradient method with successive relaxation (SSORCG), a method which is inefficient but more stable. Drawing the isostresses
The net of isostresses is drawn from the values of principal stresses computed at the Gauss points of each element. Then, they are projected and smoothed at nodes by a least square method. The isostresses are therefore drawn on the deformed mesh generation from their values at the nodes. From these values, we have built an auto-adaptive mesher giving a new meshing refined where the stresses show a rapid variation. This auto adaptive mesher is used to obtain the new meshing shown in Fig. 15 from the values of the isostresses shown in Fig. 14. Flow chart of the program
Fig. 13.
four levels of enlargement. Elsewhere, a large number of safety checks have been set in order to detect most of the anomalies (no forces, no material, etc.) as well as a sufficiently flexible structure which allows the user to retrace his steps whenever an error occurs. Examples
To illustrate the program we have determined the displacements and the isostresses of a hook with its left part is sealed in a wall. The Young's modulus is 70 000 N/mm 2 and the Poisson ratio is 0.3. We impose the force (1414 N, 1414 N) on the left node of the upper part. On the Fig. 13 below, we show the initial meshing and the deformated meshing in bold line. On the Fig. 14 we show the isostresses for this problem. In the Fig. 15 we show the new meshing obtained with the self-adaptive mesher and the new displacements.
~ - ~
Iso-su-esses with
"l"mscacriteria In Fig. 12, we present the main organization of a full multigrid cycle in association with a finite element method for the calculation of the structure. Presentation of the software
The program has been particularly centered on a conversational working mode which allows the simplest and the most visually possible way of use. Thereby, it is managed by menus and the data inputs are done on screen in a visual way (capture of various entities). The mesh generation and its deformed version are displayed graphically on the screen. They are related to a zoom which allows
Fig. 14.
Multigrid microsoftware for engineers on a Macintosh computer
225
!j
H~raber or"~des: 153
H~r~berol eler~nts: 12"/
4 blocks and 20 nodesof reference
Fig. 17.
Fig. 15.
4 THE MESH GENERATOR
ISOPAR
ISOPAR is a semi-automatic mesh generator. Its principal aim is to avoid the tiresome work of the data input. This data must be implemented in order to be consistent with the data structure of the central program. This mesh generator must be able to use multivarious finite elements and must be as conversational as possible. Principle of the meshing generation
This semi-automatic mesh generator uses a method developed from a partition of the domain in elementary parts or blocks having curved boundaries. From this decomposition, a mesh generation is automatically done, in each part, using different types of elements: triangles with 3 or 6 nodes and quadrilaterals with 4, 8 or 9 nodes. The notion of block is simple to adapt and flexible enough to various situations. For example, a circle can be seen as quadrilateral with curved boundaries, as it is shown in Fig. 16 for a quarter of a ring. Figure 17 shows a square plate with a circular hole. It is divided into 4 elementary blocks. After the domain has been decomposed into blocks, the software divides each block into quadrilaterals according to the principle of isoparametric elements. The user may ask two mesh generation densities; one with
isoparametric quadrilaterals with 4 nodes if linear elements are needed in part of the domain, and the other with isoparametric quadrilaterals with 8 nodes. Let us briefly recall the principle of mesh generation with isoparametric elements. Each real block of reference with 8 nodes is the image of a fixed block of reference, centred on [ - 1, 1]. The quadilateral mesh generator, for example, is easily obtained in the reference block. It is then transposed in the real block as shown in Fig. 18. The lines drawn in the reference element by the mesh generation lines are curves with x = constant or y = constant, In both elements the nets of drawn curves are orthogonal (hence the name 'isoparametric'). This leads to a better generation which is more balanced and of better quality. The mesh generation of each block is done automatically independent of the other blocks. Then the multiple nodes that are created at the interface are automatically suppressed in order to conform with the mesh generation in the finite element point of view. Now we have a mesh generator for the domain in quadrilaterals which can be divided into triangles. Generally speaking, to obtain the desired element, one undertakes a post-treatment from the available element (Fig. 19) To obtain 3 nodes triangle, elements of type 1 are cut in two parts according to the shortest diagonal; In the case of quadrilaterals with 4 nodes, nothing is changed; For 6 nodes triangles, first a node is added at the center
b
Thei:xa;trde2t~licOnitftheblo~,k 1t
~ ~,
O Elementof reference
Fig. 16.
Realelement Fig. 18.
226
C. Le Carlier de Veslud, G. Maurice, R. Kouitat
3noaes 1quadrilateral (Elementtype1) so~e s~ I(Elementtype2)
b~)~ 2nU'i~gles Fig. 19.
of the element which is then cut according to the shortest diagonal; For 8 nodes quadrilaterals, nothing is changed; For 9 nodes quadrilaterals, a node is added at the center of the element. A geometrical boundary code is allocated to every new created point. This boundary code is null for internal points, while for external points it sends to an edge number. Generally speaking, boundary edges are defined during the entering of blocks and are localized by a code number.
translation, dilatation). After the mesh generation, many utilities allow to cut out the resulting contour according to user convenience by a direct destruction on the screen of points and elements. This allows the creation of holes and slots in the domain. Elsewhere the mesh generation can be corrected locally by a manual displacement and an automatic paste of certain nodes, which allows to weld elements between them. Finally, it is possible to project internal or external edges on curved or straight boundaries defined by the user. In the mesh generation, the coordinates and the contour are updated automatically by the software. In order to render the software as interactive as possible, there is a central process of visualisation and graphic capture which works when needed. This process is managed by on-the-screen menus and is guided by a number of safetys which limit aberrations. Finally, the mesh generation can be visualised in different ways: side view, perspective view, shrink view (dilatation of element around its center of mass) or a combination of all these modes. After everything has been done the mesh generation is saved in a file. General chart flow t
Readingof the file of blocks or creation of this file-Entry or correction of nodes density
Input data structure The necessary points for the definition of blocks are stocked in the central memory either from a previously saved file or from a new definition of points and blocks which will be saved. New blocks can be created either by a direct 'catch' or by using geometrical transformations (symmetry, rotation). For the geometrical boundaries there is a file where all the necessary information is saved. This file is updated any time a new block is created. During the mesh generation blocks are locally saved and then destroyed at the end of this operation. From this step, the work is done only on the nodes of the mesh generation since the set of geometrical information is available. The nodes of the mesh generation since the set of geometrical information is available. The nodes and the elements are stocked in the central memory. At the output time they are then saved in a file which is compatible with the central program DIVQUAD.
File of the blocks
t I
Management of the screen (Zoom etc...)
"*
Meshing block after block Elimination of common nodes in interfaces
-*l
Possiblycutting in triangles
1 I Modification of the ~-~ meshing (Destruction, displacement, pasting)
t Final meshing -- Output
I
f
I I Saving of the ~ meshingfile
Some remarks
Some examples
The catch of blocks is made easier on the one hand by the automatic calculus of elementary shapes (such as rectangles or ovals) and on the other hand by the duplication of shapes from old blocks with the aid of some simple geometrical transformations (rotation, symmetry,
The fabric shown in the Fig. 20 has been made from a rectangle and a circle. The nodes which are in the neighbourhood of the boundaries have been welded by a 'cut and paste'. Those on the two horizontal boundaries have been projected on circles. Finally, some elements and
Multigrid microsoftware for engineers on a Macintosh computer
Fig. 20. nodes at the center of the mesh generation have been destroyed. This part represented in the Fig. 21 has been achieved from a circle. The central elements have been destroyed to create gaps. A side of each gap has been projected on a circular boundary.
5 THE MESH GENERATOR DELAUNAY The initial aim of this mesh generator was to rapidly generate a net of triangular elements from a set of given points with given values in order to draw the family of isovalues. A deep study has allowed the creation of a true automatic and interactive mesh generator which uses triangles with 3 and 6 nodes and which is compatible with the central program DIVQUAD. Principle of the mesh generator The first step in the mesh generation consists in the definition of internal and external contours of the domain. This operation is done from the points marked on the screen. The external contour is a closed univoque Nodes: 121
N~x~/~/~/Element
s: 184
227
curve defining the interior of the domain. The internal contours are any lines or curves defining the changes of mediums, or holes, etc. Therefore the contours are a set of curved or straight edges indispensable for the geometrical definition of the problem. With this data the program automatically generates intermediate points on the edges with respect to their sizes. The repartition of these points can be changed as needed by the user. Each edge can be divided linearly or logarithmically with respect to one of the extremities (this allows a progressive adaption of the mesh generation to the problem). The second step consists in the mesh generation of the set of these points by using a Delaunay's method. ~s Roughly speaking this method consists of testing the affiliation of each point to the previously formed triangles or to circles circumscribed about these triangles. According to the response, some triangles are destroyed and others are created in order to tie the point to the set of existing triangles. The process is initialized by creating a triangle with the first three points and by beginning with the fourth point. This method gives a regular mesh generation which is coherent in the finite element point of view, but valid only for the convex envelop of points. This may give rise to some problems when dealing with non convex contours. This question will be studied later in this paper. The third and last step consists in the test of homogeneity by the comparison of the sizes of elements between them. The elements that are considered to be too big are given a supplementary point in their center of mass. Then the new points are included in the mesh generation by a Delaunay's method. At this level the mesh generation is smoothed in order to increase its regularity. This smoothing is repeated until no point is added and the final mesh generation is obtained. What remains to be done is to force external edges which could have disappeared in the mesh generation and to suppress external elements on external contours in order to obtain the desired mesh generation of the domain. At this level many utilities are available. First of all, it is possible to refine the mesh generation of a part of the domain by dividing each element it contains into four parts. The remaining elements are modified such that the mesh generation is coherent. Following this, it is possible to replace the mesh generation with 3 node triangles by a one with 4 nodes quadrilaterals. This is done by the regrouping and division of the base mesh generation. Finally, if desired the user may transform linear elements (3 nodes and 4 nodes) into quadratic isoparametric elements (6 nodes and 8 nodes). It is to be noted that all these transformations take into account the definition of the boundaries, in particular that of curved boundaries. Input data structure
Fig. 21.
The edges which are necessary for the definition of the
C. Le Carlier de Veslud, G. Maurice, R. Kouitat
228
contours are stored in files beginning with the external ones according to the following process: Starting and ending points; Radius of curvature (A = 0 if the edge is straight); The coordinates of the center of the circle; The number of intermediate points on the edge (initially set to 0). In the central memory there is a table which manages the beginning and the end of each contour in the file. It is to be noted that when an external contour is defined an automatic process tests its validity. The set of points is stored in the central memory and is governed by three main parameters: (a) The number of points that strictly define the contours. This number is fixed; (b) The total number of points on the contour. This number depends on the density of chosen intermediate nodes and may vary; (c) The total number of points which vary during the mesh generation;
befo me after meshing
Fig. 23.
such that the elements are either completely in or completely out of the contours. The idea is to substitute locally for the default mesh generation an equivalent one containing the missed edge. Following are the fundamental steps: (a) Determination of all the elements cutting the edge (AB) and the corresponding points: i~
i2
These three parameters and the file of edges allow the definition of geometrical boundary codes which are indispensable for the multigrid methods. The mesh generation is kept in the central memory and is then saved in a file once the program has terminated. Convexity problems As seen previously, the Delaunay's method works only on a convex envelop of points. That is if the contour is not convex there will be some elements out of the contours. In the case where the elements are completely out of the contours they are easily detected and destroyed as shown in Fig. 22 below. However, there are cases where the elements are neither completely out nor completely in the mesh generation so that they can be destroyed by mistake. Analytically this means that an edge of the contour is not more in the mesh generation as shown in Fig. 23. The addition of points ~k has suppressed the edge. The edge (AB) is recovered by calling an edge forcing process
i3
~2
(b) Suppression of those triangles. Therefore one has two series of points: the one of points 'under' (A B) and that of points on the top of (AB). The treatment being symmetric, one will consider only those points on the top. i2
A
X
~
(c) Determination of the nearest point (1'0)to (AB) and creation of the element AioB. ii
.
i3
elements to destroy
(d) It is found that work is done recursively on two edges; Aio and ioB for which the order of critical points is low such that the convergence of the method is guaranteed. It is continued likewise until there is no intermediate point. before destruction
after destruction
Fig. 22.
In working so 'on the top' and 'under' (A B), the edge is recovered and finally the resulting elements are either completely in or completely out in relation to (AB).
Multigrid microsoftware for engineers on a Macintosh computer There is no reason for the elements to be suppressed, but for the latter it is necessary to smooth in order to improve the quality of the mesh generation. In conclusion this method when repeatedly applied on all edges of the contour, gives the exact mesh generation of the contour. From this point, one can for example, transform 3 node elements into 6 node elements.
229
of points: 89
General unfolding
of nodes: 141
The chart flow of the program is:
I Savingof I - ~ ~ [ - d theda~a I ~ / "
Envy of , ~ poin~ definitmnoftheinner outerboxt.~e~ies
Au~ma~ definitionof the
divisions and modfficalior~
Meshing from the boundary points
1
[
i
Addkngpoints of meshing, ,smmothing,%~s~sof
Menu
points chosen by the user is shown, and in Fig. 31 the calculated mesh generation.
homogeneity pf the meshing Scz~en
Fig. 24.
Is the meshix~ homogeneus
Zoom
6 3-D M E S H G E N E R A T O R
Yes
Forcing of the exterioredges • elimi~ion o! the exterior dements and srr~othi~
dk
[ Savi~ the. Dle s]~Jl~
_
Our study has been extended in 3-D with a software which includes two semi-automatic mesh generators.
I Post trai~ementin ~alXglaS6 x~des
I
Examples
In the case of Fig. 24 the mesh genereation of the domain initially had few elements. Then in expectation of a high variation of the solution in some parts of the domain, a part of the mesh generation was divided twice. Figure 25 shows the progressive refinement of the mesh generation up to the critical sector. In Figs 26 to 29 the mesh generation of a diesel engine piston's head is shown. It is designed to be studied in axisymmetric representation. In Fig. 26 the contour of the fabric is represented as given by the user. The results of mesh generation using two different densities of points are given in Figs 27 and 28. Figure 29 is a zoom of Fig. 28. Finally we present an example of logarithmic repartition of points in a triangle where two edges have been imposed to be close enough. On Fig. 30 the repartition of
Fig. 25.
230
C. Le Carlier de Veslud, G. Maurice, R. Kouitat
Diesel engine piston head: external contour
J
T
/
1 -!_
Number of nodes: 82 Number of elements: 91
1
Fig. 26.
Fig. 29.
Number of nodes: 134
The first one uses the same conformal transformation as the mesh generator ISOPAR. The basic principle is the same as in 2-D: the domain is decomposed in elementary 'cubes' with curved boundaries (20-nodes isoparametric elements). This definition is large enough to take into account many cases (so a sphere can be considered as a 'cube with curved boundaries'). According to this data and as previously mentioned the mesh generation is undertaken in each block with respect to the given node density and independently to the other blocks. Then
number od elements: 162
Fig. 27. Number of points: 71 Number of nodes: 116
[] i
I I b
,
Fig. 28.
"2--
O
~
Fig. 30.
O
O
1D
Multigrid microsoftware for engineers on a Macintosh computer
231
Fig. 31.
Fig. 33.
nodes in common interfaces are eliminated. In its actual version the mesh generator works with 8-node (straight edges) or 20-node (curved edges) isoparametric cubic elements. A post treatment allows a decomposition in 4-node tetrahedrons. The second one uses cylindrical topology. From the mesh generation of a plane surface the mesh generator constructs layers which are deduced from the reference surface by a defined geometrical transformation (translation, rotation, dilatation, compound of precedings). Degenerated cases which can appear (for example a point
on axes of rotation) are detected and treated if possible. The software works with 8-node quadrilaterals, 6-node prisms and 4-node tetrahedrons. The results of the two mesh generators can be visualised in 3-D in perspective from a point of view fixed by the user and which can be changed as needed (Fig. 32). The software allows the visualisation of the mesh generation and the external contour in normal or in shrink view which can be zoomed 5 times. Finally, it is possible to have a solid representation of the whole by elimination of hidden parts using a 'painter' method. After the calculation has been done the whole can be zoomed as needed.
Meshing vith 20 nodes elements Number of nodes: 208 Number of el~ments: 27
This cylinder is particular: the first end is a square and the second is a circle. The meshing has been made with cubic isoparametric elements with 20 nodes.
Fig. 32.
Fig. 34.
232
C. Le Carlier de Veslud, G. Maurice, R. Kouitat
Fig. 37. Fig. 35. In Fig. 33 we show an example of mesh generation using cylindrical topology. The mesh generation of the base figure, which is a ring quarter, consists of four 4-node quadrilateral elements. Following, the whole is transformed by a rotation of 180 ° around Oy axis. One can remark that the mesh generation of the whole domain consists of 8-node quadrilaterals, except in the neighbourhood of the axis of rotation where there are 6-node prisms. In Fig. 34 the mesh generation in cylindrical topology is shown with hidden parts of a helicoidal gearing. Figures 36 to 38 show three quarters of an engine piston's head. The first figure shows the mesh generation
skin in an exploded view. The two following figures take into account the hidden parts.
CONCLUSION We have presented a software for the solution of engineering equations. It consists of about 27000 program lines in its actual version. The software includes an important graphic part which is very convivial for the user. It is developed on a microcomputer and has a very personal character, since it is provided with automatic process or is auto-adaptable for the problem to be solved. This software is based on the most efficient numerical methods in association with the multigrid process. Therefore, it can treat problems with large structures. It is concerned in a modular way so that it can be adapted to various situations. In its present version, the only applications that have been developed are in elasticity. Studies on nonlinear elasticity, homogenization and flow of fluids with yield stress are being carried out.
ACKNOWLEDGEMENT Thanks is given to the Ministry of Research and Technology of the French Government for their financial support for this v:ork under contract no 86.T.0467. REFERENCES
Fig. 36.
1. Hackbusch, W. Multigrid Methods and Application Springer-Verlag, 1985. 2. Hackbusch, W. & Trottenberg, U. Multigrid methods, Proceedings Kfln Porz 1981, Springer-Verlag, 1986.
Multigrid microsoftware for engineers on a Macintosh computer
3. Briggs, W.L. A multigrid tutorial, SIAM, 1987. 4. Golub, G. & Meurant, G. R6solution num~rique de grands syst6mes, Eyrolles, 1983. 5. Modulef. Une biblioth6que modulaire d'616ments finis, INRIA, 1988. 6. Thouzot, G. & Dhatt, G. Une pr6sentation de la m6thode des 616ments finis-Maloine, 1984. 7. George, A. & Liu, J.W. Computer solution of large sparse definite systems, Prentice Hall, Series in Computational Mathematics, Cleve Moler, advisor, 1981. 8. Zienkiewicz, O.C. & Taylor, R.L. The Finite Elements Method, Vol. 1 and 2, Fourth edition. 9. Fedorenko, R.P. A relaxation method for dolving elliptic differences equations, Computational Math. and Math. Phys., 1962, 1(5), 1092-1094. 10. Hackbusch, W. Introduction to multigrid method for the numerical solution of boundary value problems. Computational Methods for Turbulent, Transonic and Viscous Flows. Proceedings of a lecture series, Von Karman Institute, 1981. 11. Southwell, R.V. Relaxation Methods in Theoretical Physics, Clarendon Press, Oxford, 12, 1956. 12. Schroder, J., Trottenberg & Reutersberg, H. R6duktionverfahren fur Diff6renzengleichungen bei Randvertaufgaben, Num~r. Math., 1976, 26, 429-459.
233
13. Bachvalov, N.S. On the convergence of a relaxation method with natural contraints of elliptic operator, Comp. Math. Math. Phys., 1966, 6(5), 101-135. 14. Astachancev, G.P. An iterative method for solving elliptic net problems, Comp. Math. Math. Phys., 1971, 11(2), 171182. 15. Brandt, A. & Dinar, N. Multigrid Solutions to Elliptic Flow Problems for P.D.E. ed. S. Parter, Academic Press, New York, 1979. 16. Stuben, K. & Trottenberg, U. On the construction of fast solvers for elliptic equations, Computational Fluids Dynamics, Lectures Notes 1982, 04-Von Karman Institute for Fluids Dynamics, Rhode Saint Genese, 1982. 17. Maitre, J.F. & Musy, F. The contraction number of a class of two level methods, an exact evaluation of some finite elements subspace and model problem, Lecture Notes in Mathematics, 1982, 960, 535-544. 18. George, P.L. & Hermeline, F. Maillage de Delaunay d'un poly~dre convexe en dimension d. Extension ~, un polyh6dre quelconque. Rapport de recherche, no. 969. lnstitut National de Recherche en lnformatique et en Automatique, 1989.
Multigrid microsoftware for engineers on a Macintosh computer C. Le Carlier de Veslud, G. Maurice & R. Kouitat Laboratoire de Science et G~nie des Surfaces, Ecole des Mines, Parc de Saurupt, 54042 Nancy Cedex, France
We present the 2-D version of a multigrid finite element software to solve engineering problems on microcomputers. Some 3-D developments are also presented. This software uses several automatic mesh generators starting from the nodes on the boundaries. It is possible to use a logarithmic scale to determine the nodes in the neighbouring rapid variation of the unknown. Our meshing technique allows to automatically modify the shapes of the irregular elements using the centre of gravity of his neighbours. A 3-D mesh generator has been also developed. It is provided with an automatic process to refine a network in the desired region and several levels of zoom. It is also self adaptable to the problem because it can give maximum information in a precise region. It works with linear and quadratic triangles and rectangles in the 2D version and with linear parallelepipeds in the 3D version. The shapes are visualised in space with hidden surfaces removed. The user data input routine with its thorough use of menus and the mouse makes this program very user friendly. At the moment, this program can only be used for a limited set of physical situations, since it only treats problems with plane elasticity and with axisymmetric deformations and constraints. The program is being extended to include problems in nonlinear elasticity as well as fluid flows with yield stress. Key words: multigrid method, finite elements, meshing.
1 INTRODUCTION Among the methods for solving partial differential equations, the finite elements method stands out with its performance and adaptability, which tends to place it as the favoured general method for solving linear systems of considerable size, often ill conditioned. The classical iterative method for solving these systems often leads complex calculations, and since the convergence rate is close to unity, a considerable computing time is required before a solution is obtained. In the 60s appeared the multigrid methods 9 and in the 80s they were developed in this context.~ It consists of a process coupled with an iterative numerical method which considerably accelerates the convergence. Their experimentation on linear, as well as nonlinear classical problems, has rapidly lead to results which may be optimal. 2 These methods are based on a sequence of numerical processes which must be fitted judiciously in order to give rise to a good package. Then the solution to each problem is a particular case. Advances in Engineering Software 0965-9978/92/$05.00 © 1992 Elsevier Science Publishers Ltd.
219
The finite difference method, simple and symmetric in space, has become the first field of application for most multigrid packages. In the case of the finite element method, though it has a greater adaptability, the space symmetry is last. 3-5'8 So the finite element method is not well adapted to the multigrid method and is used nowadays with simple or particular elements. Its implementation on microcomputer leads to the solution of difficult problems of memory and conviviality. Therefore, it is important to solve the question of data structuring and that of the creation of algorithms which allow the implementation of the method for various elements, taking in consideration a small memory and a low speed of execution. In this paper, we describe the finite element multigrid software that we have developed on a Macintosh microcomputer. We also try to put in light our original contribution. The study is divided into five parts. In the first we recall the basics of multigrid methods. In the second part we describe the main part of the program, the data structure, the different methods we have used and the main problems we have encountered and solved in the combination of the multigrid method with the finite elements method. In the following two sections we present the
220
C. L e Carlier de Veslud, G. Maurice, R. Kouitat
heart of our work: the two mesh generators called ISOPAR and D E L A U N A Y . In the last part we present the solution of some applications with this software.
2 A VIEW ON MULTIGRID METHODS
The use of the finite differences or the finite elements methods for obtaining solutions to certain problems, sometimes gives rise to big linear systems of type: 4'7 Lu =
f
where L is the discrete operator of the problem, u is the unknown vector, f is a known vector. The multigrid method is an iterative process for finding the solutions of the discrete problems. It uses a sequence of grids more and more coarse. This method requires the exact solution of the linear system on the coarsest grid, obtained by a direct method (as opposed to an iterative method which yields only approximate solutions). 1° The multigrid method is efficient for big linear systems. Depending on the size of the system, the computation of the numerical solution of the discrete problem L. = f, can be obtained using two kinds of methods. Among the most common direct methods, one can mention the Gauss methods, the Choleski and the conjuguate gradient methods. However, many disadvantages appear, especially for large systems: a high sensibility to compute time errors leading to erroneous solutions; a large occupation of memory; a very high computing time and number of operations. Then generally iterative methods are preferred. These methods give an approximation to the solution which is hoped to be better as the iterative processes. They use a fixed point model with the aid of the operator splitting: L = M - N. Therefore, one is lead to matrice iterations of the type: Mun+l
=
Southwell tried to develop methods for solving discrete problems using a coarse discretization." In the early 60s, Fedorenko combined these two aspects in a method called 'multigrid' for the solution of Poisson's equation, r2 Then Bakhavalor undertook the study of problems with variable coefficients, j3 Next, Astrachanov proposed a multigrid formulation for finite element problems and he showed the low cost of the methodJ 4 From the 70s, Brandt devised general theories which include earlier works and he proposed efficient solution methods for certain problems, essentially those of fluid mechanics. ~5 From the end of the 70s considerable work has been done on multigrids by Hackbush, ~': Stfiben & Trotenberg, J6 as well as by Maitre in France. ~7 Principle o f the method
In order to describe the main properties of the method, let us use the same one-dimensional example as hackbush. 2 d2u dx: + f
=
x ~ ]0, 1[, u(o)
0
=
u(1)
=
0
The discretization of the problem in N steps with size h and the use of boundary conditions lead to a system of the form: Lh llh :
f t,
First stage."
To fix our ideas, let us solve the linear system by the relaxed Jacobi method (~o = 1/2). This method leads to a slow decrease of the global numerical error which is a discrete signal involving high and low frequencies. Both types of frequencies have different behaviour during iterations. It is shown mathematically that high frequencies present a convergence rate of about 0. Therefore, their effect is eliminated after some iterations. During this stage, where smoothing takes place, the error has not decreased significantly but has taken a regular aspect.
Nun + f
where u, is the solution at the nth iteration of the method and f is the second member of the initial equation. M is chosen so that the successive linear system are easily solved. The most known are those of Jacobi and GaussSeidel. Although a judicious choice of M and N gives rise to mathematical convergence, numerically the convergence is quite slow and divergence may occur. As far back as 1940, the Fourier analysis has contributed to the discovery of one of the essential properties of iterative methods: the rapid elimination of high frequencies of Un during iterations, that is the property known as smoothing, jj Elsewhere, Schr6der with his reduction method and
Second stage:
Let fih be an approximation ofuh obtained after v smoothing iterations, and, eh =fih - uh. One obtains Lheh = Lhfih-- Lhuh = Lhfih -- fh = dh. Since fih and fh are known, so is d h and then eh is the solution of the system Lheh = dh which is similar to the base system. Let consider a much coarser discretization with steps of size 2h (Fig. 1). 0 []
"~1"--2h--"'~4 h ~ o
[]
o
[] Fig. 1.
o
1 []
o
[]
Multigrid
m i c r o s o f t w a r e f o r e n g i n e e r s on a M a c i n t o s h
starting grid
starting grid
coarse grid
intermediate grid
O smoothing '5~Restriction, 7( Prolongation,[ ] exactresolution
coarse grid
Fig. 2.
Fig.
On this grid, the error is a solution of L2he2h = d2h where L2h is the same form a s L h and dEh is the restriction of dh on the coarse grid. Since the number of unknowns is small compared with those of the former grid one can compute exactly e2h = L~Jd2h and then by a p interpolation recover a new approximation of e h. So fl h : I1h ~ pe2h constitutes a new approximation of Uh which is much more precise than fih since the initial error has been corrected. This second stage is called the correction on a coarse grid. This cycle based on the use of two grids is the simplest multigrid method. It involves a smoothing followed by correction on a coarse grid (restriction on the coarse grid, exact calculation of the e r r o r LEh : L2h~d2h , extension of e2h on the sharp grid, calculation of fib). A useful way to represent this cycle is sketched in Fig. 2.
221
computer
4.
approximation fil = Pfit-~ is obtained. The multigrid process can be applied at any step for more efficiency. The full multigrid process has three main advantages: (1) l is not fixed a p r i o r i . This makes the mesh refinement a dynamic process. (2) A better accuracy than that of the discretization error can be obtained by using two distinct operators Lh and lh for the based iterations and for the calculus of residues (a method with defect correction). (3) Very good results are obtained for nonlinear problems for which the iterative process requires a good initial estimate of the solution. Therefore, the number of operations tends to be optimal. 3 THE MAIN PART OF THE PROGRAM: DIVQUAD
Multigrid method Frequently, in large size finite element systems the convergence acceleration obtained from the bigrid cycle is insufficient and the size of the linear systems on the coarse grid is too big to be solved by a direct method. Therefore, one combines smoothing and restriction until a sufficiently coarse grid is obtained, which allows an exact calculation of the error. From this point one goes back to the original sharp grid by multiple extensions. Figure 3 represents a cycle with 3 grids. Two or more corrections can be done at each level in order t. As for every iterative method the efficacy of the multigrid version the method constructs a good initial approximation in the following way. Consider a multigrid process with ! + 1 grids numbered from 0 to 1 where 0 stands for the most coarse grid; on which the exact calculation is performed; and 1 for the sharpest grid. At the zero level u0 = Lo ~f0 can be obtained exactly. A coarse approximation fi~ = pu0 on grid 1 is obtained by interpolation of u0 on it. The improvement of this solution by the use of a multigrid process leads to the vector fi~. Then u2 = PUt is calculated and what has been done for u~ is repeated and so on. At level ! a good initial
starting grid intermediategrid coarse grid
Fig.
3.
The multigrid methods as well as their implementation on computers, which requires a more particular input data structure than that of finite elements, still stand as research work. The difficulty is increased by the wish for generality of the program. DIVQUAD; written in F O R T R A N 77 language and implemented in an Apple Macintosh computer partly solves this problem. For the moment it treats plane elasticity problems (plane stresses or deformations and axisymmetric). Input data structure The multigrid method is widely recursive. The full multigrid version offers the advantage to vary the size which is unknown a p r i o r i . Therefore, it gives a solution to the problem of memory saving. From a strictly geometrical point of the view, we note that the method starts on the most coarse grid and refines step after step. The exact geometric definition of the domain is therefore necessary from the start especially concerning curved boundaries. Unknown and second member The unknown vector U and the second member B are dimensioned as large size 'tables'. After, they are cut into parts relative to each grid as one does along the progression of grids as shown in Fig. 5 below. With regard to the unknown vector U, one goes from grid i + 1 and uses the space of grid i for the solution of the correction equation in the descending phase (restriction). In the ascending phase the space of the grid i is used
222
C. Le Carlier de Veslud, G. Maurice, R. Kouitat gridi +
gridi + I ++
networkof the grid 1 +
networkof the grid 2
VectorU NTI elements NT2 elements
Fig. 5.
Fig. 7.
for correction by interpolation on grid i + 1. For the vector B, the space of the grid i + 1 is intended for the second member while that of grid i is for the default and
ensure the existence of a knot of grids wherever the force applies). This file is fixed and is common to all grids. There is an evolution of the situation during computation in the case of the two other types of forces. Effectively, the notion of loaded elements or loaded edges varies with the grids since these entities are subdivided. Then we have two files which vary as one goes along the display of grids. At the most coarse level they contain only the bare minimum for the description of the loading. Later the loading on the second grid is added, then that of the third and so on. As before, a set of addresses allows the reading of the loading which corresponds to the grid as shown in Fig. 8. It is to be noted that the most coarse grid must be sufficiently accurate in order to represent the loading of the problem under study.
SO o n .
When a new grid i + 1 is constructed geometrically the vector Ui+t is deduced from U~ by interpolation and Bi+ ~ is constructed according to 'opposing forces'. Therefore, whenever a new grid is constructed, the main parameters (numbers of points, number of elements, etc.) are memorised (saved). Therefore, the number of unknowns on all grids in use is known. An elementary procedure gives at each level the starting and ending addresses of the grid space. Mesh generation
In this case as well there is the problem of how to differentiate the grids. With regard to points the solution is simple: when a new grid is constructed the only points that are saved are those which appear for the first time making the table of points into the structure of the sketch on Fig. 6 below. The first N t points are in grid 1 and the first Nz points in grid 2. The points having numbers between NI and N2 have appeared recently and are in grid 2 but not in grid 1. Such an artifice cannot be applied for mesh generations which are saved in consecutive files as can be seen in Fig. 7. As previously mentioned, an archivage procedure allows the recovery of starting and ending addresses of points and elements at any level. Applied forces
In the case of a plane problem applied forces are of three types: punctual forces, lineal forces and surface forces. With regards to the first forces there are no particular difficulties and they are saved in a file having the size of the most coarse grid (however care must be taken to
grid n°2j N2 points )
Boundary conditions
Two types of boundary conditions must be distinguished: (1) Geometrical boundary conditions that are tied to a point and which allow the placing of points on the boundary and during displays of grids to create new points in conformity with boundaries (curved boundaries for example); (2) Numerical boundary conditions which express the imposed displacements at some points (lockings for example). However, these boundaries also arise during the display of a grid. For example if a new point is created on the boundary between two locked points, it is also locked unless there is an explicit different order. Geometrical boundary conditions are defined through an edge file describing the external contour. The edges defined by the mesh generators and are divided into two types: straight edges and curved edges. Each point is File
Loadedelement
<;>
t
gridn° 1 Loading
4 Loadedelement s
Y
Loading gridn° 2 !
Loading grid n° 3
grin n° l ( N1 points ) Fig. 6.
Fig. 8.
Multigrid microsoftware for engineers on a Macintosh computer IFR NO2X 1"
NO2Y UX UY 1" t
I geometric boundary code
223
I eventaalyimposeddisplacements
Fig. 9.
assigned a boundary code related if needed to the edges. Internal points have a null boundary code. Numerical boundary conditions are stored in memory by two locking codes (one on the x and the other on the y) which are null if there is no imposed displacement and have the value 1 otherwise. This information is gathered in a file having the same size as that of the points and varying according to the grids. Each point involves the locking codes in x and y according to Fig. 9.
Stocking of the stiffness matrix
We have used a M O R S E type stocking. The topologically no null elements of the matrix are compacted line by line then put end to end in a one-dimensional table. A table o f pointers gives the location of the beginning and the end of each line. A second table gives the null location of the compacted term for each line. It is a particularly efficient method for large sparse matrices since it allows a rapid pass through the file which is well adapted to the solution method (conjugate gradient for example). Elsewhere the two tables of 'pointers' give precious indications on the mesh generation (very useful for example in the reduction of the band width). The materials and the stresses
The materials are defined again each time the program is used. It is possible to use many materials the features of which are stocked in sequential files. The stresses are put into files in such a way that their screen drawing is possible later. Grid divisions
The basic idea of the software is to closely relate the solution and the mesh generation. For this reason we use the 'full multigrid' process described previously. The solution of the problem on the coarse grid permits the calculation of stresses at Gauss points. These are projected on the edges by a least square method. Since now, the variations of stresses on each element compared with the mean variations on the whole domain give us precious information on the mesh generation quality. The elements that are considered to be to stressed are then divided into four or nine and the others are adjusted such that a link is made with the total mesh generation. This process leads to a new data base according to the new elements and features so created (components,
/\ undivided startingmeshing
divided meshing
Fig. 10.
boundaries, limits, etc.). With regard to linear elements, the created points must eventually be projected on curved boundaries (necessity of boundary codes). The quadratic elements already involve the notion of curvature under their isoparametric form so that the projection is not necessary (Fig. 11). 8 After having verified that each new point does not already exist, the new elements are then always constructed in the same order so as to make location easier. At this level, files of linear and surface forces are made complete in order to take into account new edges and new elements created in the new grid. This new grid is stocked according to the data structure described previously. A similar method is used during the ascending of grids in the multigrid process. The corrected error on the coarse grid is interpolated on the sharp grid by using proper shape functions of each element. Then, one obtains a correction scheme which can be used for any element. Numerical smoothers and solvers
The smoothing is actually done by a relaxed GaussSeidel method with 5 or 10 iterations. A smoothing by a relaxed Jacobi method is also possible but it is less efficient. However, we must note that this last method allows smoothing without the necessity of stocking the stiffness matrix which is computed again after each iteration. The solution at the most coarse level, is obtained by using one of the most efficient numerical methods which is a preconditioned conjugate gradient method with elementsto dividein 4
\\ initial meshing
\\
final meshing afterdivis ~n in 4 and connectionto the rest of the meshing Fig. 11.
224
C. Le Carlier de Veslud, G. Maurice, R. Kouitat Re~ling of the coarse [ meshing - Resolution
coarse grkl []
of the grid [ in~rpolation ] coarsegrid . . . .
refinement
Smoothing [ smoothgrkl -- - ~ cornputm~o!t~ defect resign
coarse
grid . . . . . . . . . . . .
Resolutionof the correction srr~o~ grid '-----~----~,i I equation-Coxrec~n [ co~se grid
l
The precisionof the ][ computin~is enough?
~
Yes
l( -~y~b~s [] Exit ~ol~i~ 0
Graphicoutputoftheresults[
Smoot~ ~xnl:~n
Fig. 12.
incomplete Choleski factorization (ICCG). Nevertheless, this method is not always stable for elements of degree 2. Therefore, the solution is done by a conjugate gradient method with successive relaxation (SSORCG), a method which is inefficient but more stable. Drawing the isostresses
The net of isostresses is drawn from the values of principal stresses computed at the Gauss points of each element. Then, they are projected and smoothed at nodes by a least square method. The isostresses are therefore drawn on the deformed mesh generation from their values at the nodes. From these values, we have built an auto-adaptive mesher giving a new meshing refined where the stresses show a rapid variation. This auto adaptive mesher is used to obtain the new meshing shown in Fig. 15 from the values of the isostresses shown in Fig. 14. Flow chart of the program
Fig. 13.
four levels of enlargement. Elsewhere, a large number of safety checks have been set in order to detect most of the anomalies (no forces, no material, etc.) as well as a sufficiently flexible structure which allows the user to retrace his steps whenever an error occurs. Examples
To illustrate the program we have determined the displacements and the isostresses of a hook with its left part is sealed in a wall. The Young's modulus is 70 000 N/mm 2 and the Poisson ratio is 0.3. We impose the force (1414 N, 1414 N) on the left node of the upper part. On the Fig. 13 below, we show the initial meshing and the deformated meshing in bold line. On the Fig. 14 we show the isostresses for this problem. In the Fig. 15 we show the new meshing obtained with the self-adaptive mesher and the new displacements.
~ - ~
Iso-su-esses with
"l"mscacriteria In Fig. 12, we present the main organization of a full multigrid cycle in association with a finite element method for the calculation of the structure. Presentation of the software
The program has been particularly centered on a conversational working mode which allows the simplest and the most visually possible way of use. Thereby, it is managed by menus and the data inputs are done on screen in a visual way (capture of various entities). The mesh generation and its deformed version are displayed graphically on the screen. They are related to a zoom which allows
Fig. 14.
Multigrid microsoftware for engineers on a Macintosh computer
225
!j
H~raber or"~des: 153
H~r~berol eler~nts: 12"/
4 blocks and 20 nodesof reference
Fig. 17.
Fig. 15.
4 THE MESH GENERATOR
ISOPAR
ISOPAR is a semi-automatic mesh generator. Its principal aim is to avoid the tiresome work of the data input. This data must be implemented in order to be consistent with the data structure of the central program. This mesh generator must be able to use multivarious finite elements and must be as conversational as possible. Principle of the meshing generation
This semi-automatic mesh generator uses a method developed from a partition of the domain in elementary parts or blocks having curved boundaries. From this decomposition, a mesh generation is automatically done, in each part, using different types of elements: triangles with 3 or 6 nodes and quadrilaterals with 4, 8 or 9 nodes. The notion of block is simple to adapt and flexible enough to various situations. For example, a circle can be seen as quadrilateral with curved boundaries, as it is shown in Fig. 16 for a quarter of a ring. Figure 17 shows a square plate with a circular hole. It is divided into 4 elementary blocks. After the domain has been decomposed into blocks, the software divides each block into quadrilaterals according to the principle of isoparametric elements. The user may ask two mesh generation densities; one with
isoparametric quadrilaterals with 4 nodes if linear elements are needed in part of the domain, and the other with isoparametric quadrilaterals with 8 nodes. Let us briefly recall the principle of mesh generation with isoparametric elements. Each real block of reference with 8 nodes is the image of a fixed block of reference, centred on [ - 1, 1]. The quadilateral mesh generator, for example, is easily obtained in the reference block. It is then transposed in the real block as shown in Fig. 18. The lines drawn in the reference element by the mesh generation lines are curves with x = constant or y = constant, In both elements the nets of drawn curves are orthogonal (hence the name 'isoparametric'). This leads to a better generation which is more balanced and of better quality. The mesh generation of each block is done automatically independent of the other blocks. Then the multiple nodes that are created at the interface are automatically suppressed in order to conform with the mesh generation in the finite element point of view. Now we have a mesh generator for the domain in quadrilaterals which can be divided into triangles. Generally speaking, to obtain the desired element, one undertakes a post-treatment from the available element (Fig. 19) To obtain 3 nodes triangle, elements of type 1 are cut in two parts according to the shortest diagonal; In the case of quadrilaterals with 4 nodes, nothing is changed; For 6 nodes triangles, first a node is added at the center
b
Thei:xa;trde2t~licOnitftheblo~,k 1t
~ ~,
O Elementof reference
Fig. 16.
Realelement Fig. 18.
226
C. Le Carlier de Veslud, G. Maurice, R. Kouitat
3noaes 1quadrilateral (Elementtype1) so~e s~ I(Elementtype2)
b~)~ 2nU'i~gles Fig. 19.
of the element which is then cut according to the shortest diagonal; For 8 nodes quadrilaterals, nothing is changed; For 9 nodes quadrilaterals, a node is added at the center of the element. A geometrical boundary code is allocated to every new created point. This boundary code is null for internal points, while for external points it sends to an edge number. Generally speaking, boundary edges are defined during the entering of blocks and are localized by a code number.
translation, dilatation). After the mesh generation, many utilities allow to cut out the resulting contour according to user convenience by a direct destruction on the screen of points and elements. This allows the creation of holes and slots in the domain. Elsewhere the mesh generation can be corrected locally by a manual displacement and an automatic paste of certain nodes, which allows to weld elements between them. Finally, it is possible to project internal or external edges on curved or straight boundaries defined by the user. In the mesh generation, the coordinates and the contour are updated automatically by the software. In order to render the software as interactive as possible, there is a central process of visualisation and graphic capture which works when needed. This process is managed by on-the-screen menus and is guided by a number of safetys which limit aberrations. Finally, the mesh generation can be visualised in different ways: side view, perspective view, shrink view (dilatation of element around its center of mass) or a combination of all these modes. After everything has been done the mesh generation is saved in a file. General chart flow t
Readingof the file of blocks or creation of this file-Entry or correction of nodes density
Input data structure The necessary points for the definition of blocks are stocked in the central memory either from a previously saved file or from a new definition of points and blocks which will be saved. New blocks can be created either by a direct 'catch' or by using geometrical transformations (symmetry, rotation). For the geometrical boundaries there is a file where all the necessary information is saved. This file is updated any time a new block is created. During the mesh generation blocks are locally saved and then destroyed at the end of this operation. From this step, the work is done only on the nodes of the mesh generation since the set of geometrical information is available. The nodes of the mesh generation since the set of geometrical information is available. The nodes and the elements are stocked in the central memory. At the output time they are then saved in a file which is compatible with the central program DIVQUAD.
File of the blocks
t I
Management of the screen (Zoom etc...)
"*
Meshing block after block Elimination of common nodes in interfaces
-*l
Possiblycutting in triangles
1 I Modification of the ~-~ meshing (Destruction, displacement, pasting)
t Final meshing -- Output
I
f
I I Saving of the ~ meshingfile
Some remarks
Some examples
The catch of blocks is made easier on the one hand by the automatic calculus of elementary shapes (such as rectangles or ovals) and on the other hand by the duplication of shapes from old blocks with the aid of some simple geometrical transformations (rotation, symmetry,
The fabric shown in the Fig. 20 has been made from a rectangle and a circle. The nodes which are in the neighbourhood of the boundaries have been welded by a 'cut and paste'. Those on the two horizontal boundaries have been projected on circles. Finally, some elements and
Multigrid microsoftware for engineers on a Macintosh computer
Fig. 20. nodes at the center of the mesh generation have been destroyed. This part represented in the Fig. 21 has been achieved from a circle. The central elements have been destroyed to create gaps. A side of each gap has been projected on a circular boundary.
5 THE MESH GENERATOR DELAUNAY The initial aim of this mesh generator was to rapidly generate a net of triangular elements from a set of given points with given values in order to draw the family of isovalues. A deep study has allowed the creation of a true automatic and interactive mesh generator which uses triangles with 3 and 6 nodes and which is compatible with the central program DIVQUAD. Principle of the mesh generator The first step in the mesh generation consists in the definition of internal and external contours of the domain. This operation is done from the points marked on the screen. The external contour is a closed univoque Nodes: 121
N~x~/~/~/Element
s: 184
227
curve defining the interior of the domain. The internal contours are any lines or curves defining the changes of mediums, or holes, etc. Therefore the contours are a set of curved or straight edges indispensable for the geometrical definition of the problem. With this data the program automatically generates intermediate points on the edges with respect to their sizes. The repartition of these points can be changed as needed by the user. Each edge can be divided linearly or logarithmically with respect to one of the extremities (this allows a progressive adaption of the mesh generation to the problem). The second step consists in the mesh generation of the set of these points by using a Delaunay's method. ~s Roughly speaking this method consists of testing the affiliation of each point to the previously formed triangles or to circles circumscribed about these triangles. According to the response, some triangles are destroyed and others are created in order to tie the point to the set of existing triangles. The process is initialized by creating a triangle with the first three points and by beginning with the fourth point. This method gives a regular mesh generation which is coherent in the finite element point of view, but valid only for the convex envelop of points. This may give rise to some problems when dealing with non convex contours. This question will be studied later in this paper. The third and last step consists in the test of homogeneity by the comparison of the sizes of elements between them. The elements that are considered to be too big are given a supplementary point in their center of mass. Then the new points are included in the mesh generation by a Delaunay's method. At this level the mesh generation is smoothed in order to increase its regularity. This smoothing is repeated until no point is added and the final mesh generation is obtained. What remains to be done is to force external edges which could have disappeared in the mesh generation and to suppress external elements on external contours in order to obtain the desired mesh generation of the domain. At this level many utilities are available. First of all, it is possible to refine the mesh generation of a part of the domain by dividing each element it contains into four parts. The remaining elements are modified such that the mesh generation is coherent. Following this, it is possible to replace the mesh generation with 3 node triangles by a one with 4 nodes quadrilaterals. This is done by the regrouping and division of the base mesh generation. Finally, if desired the user may transform linear elements (3 nodes and 4 nodes) into quadratic isoparametric elements (6 nodes and 8 nodes). It is to be noted that all these transformations take into account the definition of the boundaries, in particular that of curved boundaries. Input data structure
Fig. 21.
The edges which are necessary for the definition of the
C. Le Carlier de Veslud, G. Maurice, R. Kouitat
228
contours are stored in files beginning with the external ones according to the following process: Starting and ending points; Radius of curvature (A = 0 if the edge is straight); The coordinates of the center of the circle; The number of intermediate points on the edge (initially set to 0). In the central memory there is a table which manages the beginning and the end of each contour in the file. It is to be noted that when an external contour is defined an automatic process tests its validity. The set of points is stored in the central memory and is governed by three main parameters: (a) The number of points that strictly define the contours. This number is fixed; (b) The total number of points on the contour. This number depends on the density of chosen intermediate nodes and may vary; (c) The total number of points which vary during the mesh generation;
befo me after meshing
Fig. 23.
such that the elements are either completely in or completely out of the contours. The idea is to substitute locally for the default mesh generation an equivalent one containing the missed edge. Following are the fundamental steps: (a) Determination of all the elements cutting the edge (AB) and the corresponding points: i~
i2
These three parameters and the file of edges allow the definition of geometrical boundary codes which are indispensable for the multigrid methods. The mesh generation is kept in the central memory and is then saved in a file once the program has terminated. Convexity problems As seen previously, the Delaunay's method works only on a convex envelop of points. That is if the contour is not convex there will be some elements out of the contours. In the case where the elements are completely out of the contours they are easily detected and destroyed as shown in Fig. 22 below. However, there are cases where the elements are neither completely out nor completely in the mesh generation so that they can be destroyed by mistake. Analytically this means that an edge of the contour is not more in the mesh generation as shown in Fig. 23. The addition of points ~k has suppressed the edge. The edge (AB) is recovered by calling an edge forcing process
i3
~2
(b) Suppression of those triangles. Therefore one has two series of points: the one of points 'under' (A B) and that of points on the top of (AB). The treatment being symmetric, one will consider only those points on the top. i2
A
X
~
(c) Determination of the nearest point (1'0)to (AB) and creation of the element AioB. ii
.
i3
elements to destroy
(d) It is found that work is done recursively on two edges; Aio and ioB for which the order of critical points is low such that the convergence of the method is guaranteed. It is continued likewise until there is no intermediate point. before destruction
after destruction
Fig. 22.
In working so 'on the top' and 'under' (A B), the edge is recovered and finally the resulting elements are either completely in or completely out in relation to (AB).
Multigrid microsoftware for engineers on a Macintosh computer There is no reason for the elements to be suppressed, but for the latter it is necessary to smooth in order to improve the quality of the mesh generation. In conclusion this method when repeatedly applied on all edges of the contour, gives the exact mesh generation of the contour. From this point, one can for example, transform 3 node elements into 6 node elements.
229
of points: 89
General unfolding
of nodes: 141
The chart flow of the program is:
I Savingof I - ~ ~ [ - d theda~a I ~ / "
Envy of , ~ poin~ definitmnoftheinner outerboxt.~e~ies
Au~ma~ definitionof the
divisions and modfficalior~
Meshing from the boundary points
1
[
i
Addkngpoints of meshing, ,smmothing,%~s~sof
Menu
points chosen by the user is shown, and in Fig. 31 the calculated mesh generation.
homogeneity pf the meshing Scz~en
Fig. 24.
Is the meshix~ homogeneus
Zoom
6 3-D M E S H G E N E R A T O R
Yes
Forcing of the exterioredges • elimi~ion o! the exterior dements and srr~othi~
dk
[ Savi~ the. Dle s]~Jl~
_
Our study has been extended in 3-D with a software which includes two semi-automatic mesh generators.
I Post trai~ementin ~alXglaS6 x~des
I
Examples
In the case of Fig. 24 the mesh genereation of the domain initially had few elements. Then in expectation of a high variation of the solution in some parts of the domain, a part of the mesh generation was divided twice. Figure 25 shows the progressive refinement of the mesh generation up to the critical sector. In Figs 26 to 29 the mesh generation of a diesel engine piston's head is shown. It is designed to be studied in axisymmetric representation. In Fig. 26 the contour of the fabric is represented as given by the user. The results of mesh generation using two different densities of points are given in Figs 27 and 28. Figure 29 is a zoom of Fig. 28. Finally we present an example of logarithmic repartition of points in a triangle where two edges have been imposed to be close enough. On Fig. 30 the repartition of
Fig. 25.
230
C. Le Carlier de Veslud, G. Maurice, R. Kouitat
Diesel engine piston head: external contour
J
T
/
1 -!_
Number of nodes: 82 Number of elements: 91
1
Fig. 26.
Fig. 29.
Number of nodes: 134
The first one uses the same conformal transformation as the mesh generator ISOPAR. The basic principle is the same as in 2-D: the domain is decomposed in elementary 'cubes' with curved boundaries (20-nodes isoparametric elements). This definition is large enough to take into account many cases (so a sphere can be considered as a 'cube with curved boundaries'). According to this data and as previously mentioned the mesh generation is undertaken in each block with respect to the given node density and independently to the other blocks. Then
number od elements: 162
Fig. 27. Number of points: 71 Number of nodes: 116
[] i
I I b
,
Fig. 28.
"2--
O
~
Fig. 30.
O
O
1D
Multigrid microsoftware for engineers on a Macintosh computer
231
Fig. 31.
Fig. 33.
nodes in common interfaces are eliminated. In its actual version the mesh generator works with 8-node (straight edges) or 20-node (curved edges) isoparametric cubic elements. A post treatment allows a decomposition in 4-node tetrahedrons. The second one uses cylindrical topology. From the mesh generation of a plane surface the mesh generator constructs layers which are deduced from the reference surface by a defined geometrical transformation (translation, rotation, dilatation, compound of precedings). Degenerated cases which can appear (for example a point
on axes of rotation) are detected and treated if possible. The software works with 8-node quadrilaterals, 6-node prisms and 4-node tetrahedrons. The results of the two mesh generators can be visualised in 3-D in perspective from a point of view fixed by the user and which can be changed as needed (Fig. 32). The software allows the visualisation of the mesh generation and the external contour in normal or in shrink view which can be zoomed 5 times. Finally, it is possible to have a solid representation of the whole by elimination of hidden parts using a 'painter' method. After the calculation has been done the whole can be zoomed as needed.
Meshing vith 20 nodes elements Number of nodes: 208 Number of el~ments: 27
This cylinder is particular: the first end is a square and the second is a circle. The meshing has been made with cubic isoparametric elements with 20 nodes.
Fig. 32.
Fig. 34.
232
C. Le Carlier de Veslud, G. Maurice, R. Kouitat
Fig. 37. Fig. 35. In Fig. 33 we show an example of mesh generation using cylindrical topology. The mesh generation of the base figure, which is a ring quarter, consists of four 4-node quadrilateral elements. Following, the whole is transformed by a rotation of 180 ° around Oy axis. One can remark that the mesh generation of the whole domain consists of 8-node quadrilaterals, except in the neighbourhood of the axis of rotation where there are 6-node prisms. In Fig. 34 the mesh generation in cylindrical topology is shown with hidden parts of a helicoidal gearing. Figures 36 to 38 show three quarters of an engine piston's head. The first figure shows the mesh generation
skin in an exploded view. The two following figures take into account the hidden parts.
CONCLUSION We have presented a software for the solution of engineering equations. It consists of about 27000 program lines in its actual version. The software includes an important graphic part which is very convivial for the user. It is developed on a microcomputer and has a very personal character, since it is provided with automatic process or is auto-adaptable for the problem to be solved. This software is based on the most efficient numerical methods in association with the multigrid process. Therefore, it can treat problems with large structures. It is concerned in a modular way so that it can be adapted to various situations. In its present version, the only applications that have been developed are in elasticity. Studies on nonlinear elasticity, homogenization and flow of fluids with yield stress are being carried out.
ACKNOWLEDGEMENT Thanks is given to the Ministry of Research and Technology of the French Government for their financial support for this v:ork under contract no 86.T.0467. REFERENCES
Fig. 36.
1. Hackbusch, W. Multigrid Methods and Application Springer-Verlag, 1985. 2. Hackbusch, W. & Trottenberg, U. Multigrid methods, Proceedings Kfln Porz 1981, Springer-Verlag, 1986.
Multigrid microsoftware for engineers on a Macintosh computer
3. Briggs, W.L. A multigrid tutorial, SIAM, 1987. 4. Golub, G. & Meurant, G. R6solution num~rique de grands syst6mes, Eyrolles, 1983. 5. Modulef. Une biblioth6que modulaire d'616ments finis, INRIA, 1988. 6. Thouzot, G. & Dhatt, G. Une pr6sentation de la m6thode des 616ments finis-Maloine, 1984. 7. George, A. & Liu, J.W. Computer solution of large sparse definite systems, Prentice Hall, Series in Computational Mathematics, Cleve Moler, advisor, 1981. 8. Zienkiewicz, O.C. & Taylor, R.L. The Finite Elements Method, Vol. 1 and 2, Fourth edition. 9. Fedorenko, R.P. A relaxation method for dolving elliptic differences equations, Computational Math. and Math. Phys., 1962, 1(5), 1092-1094. 10. Hackbusch, W. Introduction to multigrid method for the numerical solution of boundary value problems. Computational Methods for Turbulent, Transonic and Viscous Flows. Proceedings of a lecture series, Von Karman Institute, 1981. 11. Southwell, R.V. Relaxation Methods in Theoretical Physics, Clarendon Press, Oxford, 12, 1956. 12. Schroder, J., Trottenberg & Reutersberg, H. R6duktionverfahren fur Diff6renzengleichungen bei Randvertaufgaben, Num~r. Math., 1976, 26, 429-459.
233
13. Bachvalov, N.S. On the convergence of a relaxation method with natural contraints of elliptic operator, Comp. Math. Math. Phys., 1966, 6(5), 101-135. 14. Astachancev, G.P. An iterative method for solving elliptic net problems, Comp. Math. Math. Phys., 1971, 11(2), 171182. 15. Brandt, A. & Dinar, N. Multigrid Solutions to Elliptic Flow Problems for P.D.E. ed. S. Parter, Academic Press, New York, 1979. 16. Stuben, K. & Trottenberg, U. On the construction of fast solvers for elliptic equations, Computational Fluids Dynamics, Lectures Notes 1982, 04-Von Karman Institute for Fluids Dynamics, Rhode Saint Genese, 1982. 17. Maitre, J.F. & Musy, F. The contraction number of a class of two level methods, an exact evaluation of some finite elements subspace and model problem, Lecture Notes in Mathematics, 1982, 960, 535-544. 18. George, P.L. & Hermeline, F. Maillage de Delaunay d'un poly~dre convexe en dimension d. Extension ~, un polyh6dre quelconque. Rapport de recherche, no. 969. lnstitut National de Recherche en lnformatique et en Automatique, 1989.