Wavelets in hybrid-mixed stress elements

Comput. Methods Appl. Mech. Engrg. 190 (2001) 3977±3998

www.elsevier.com/locate/cma

Wavelets in hybrid-mixed stress elements Luõs Manuel Santos Castro *, Jo~ ao Ant onio Teixeira de Freitas Departamento de Engenharia Civil, Instituto Superior T ecnico, Av Rovisco Pais, 1049-001 Lisboa, Portugal Received 15 July 1999; received in revised form 18 May 2000

Abstract The objective of this paper is to present the application of wavelets in the implementation of the stress model of the hybrid-mixed ®nite element formulation. Independent wavelet bases are used to approximate the stresses in the domain and the displacements both in the domain and on the boundary of elastostatic ®nite elements. Except for the kinematic boundary conditions, which are enforced directly, all the remaining equations of elastostatics are enforced in a weighted residual form so designed as to ensure that the discrete model embodies the relevant properties of the continuum ®elds they represent, namely static-kinematic duality and elastic reciprocity. Numerical examples are used to illustrate the characteristics of the numerical model being presented and to assess its accuracy and eciency. Ó 2001 Elsevier Science B.V. All rights reserved.

1. Introduction The non-conventional ®nite element formulations can be classi®ed into three groups, namely hybridmixed, hybrid and hybrid-Tre€tz, and two models can be distinguished in each formulation [1,2]. The alternative stress and displacement models are designed to produce, under certain conditions, statically and kinematically admissible solutions, respectively. What distinguishes the alternative formulations are the local conditions that the domain approximation bases are constrained to satisfy a priori. The approximation bases of the hybrid-Tre€tz formulation must locally solve all the fundamental domain conditions. The hybrid formulation is obtained by relaxing this constraint to a subset of the fundamental conditions of the problem, typically either the domain equilibrium or the compatibility conditions. The hybrid-mixed formulation results from accepting approximation bases that locally satisfy none of the domain conditions of the problem under analysis. The price paid by relaxing completely the constraints on the domain approximation bases is that atleast three ®elds must be approximated simultaneously. These are the stresses and the displacements in the domain of the element and the displacements on its boundary in the hybrid-mixed stress model for the elastostatic problems considered here. The implied cost is a substantial increase in the number of degrees of freedom of the model, which can only be compensated by using computer friendly approximation bases, i.e. bases that produce sparse solving systems with coef®cients that can be computed at very low cost. The bases that have been studied for the implementation of the hybrid-mixed formulation are orthogonal (Legendre) polynomials and digital functions, namely Walsh functions and wavelets. The experience in these applications is still limited to the elastic, elastodynamic and elastoplastic analyses of stretching and bending plates and of solids [3±5].

*

Corresponding author. Fax: +351-21-849-7650. E-mail address: [email protected] (L.M.S. Castro).

0045-7825/01/$ - see front matter Ó 2001 Elsevier Science B.V. All rights reserved. PII: S 0 0 4 5 - 7 8 2 5 ( 0 0 ) 0 0 3 1 3 - 3

3978

L.M.S. Castro, J.A. Teixeira de Freitas / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3977±3998

The objective of the present paper is to report the use of wavelets as approximation functions in the implementation of the stress model of the hybrid-mixed ®nite element formulation. The approach followed is therefore distinct from that used in the development of the Wavelet Element Method, WEM, as described in [6,7]. The WEM allows for the construction of multiresolution systems and biorthogonal wavelets on general domains. Tensor products involving scaling functions and wavelets de®ned on the interval are used to de®ne the basis in the reference domain that has to be mapped into the di€erent elements in which the original domain was subdivided. The globally continuous biorthogonal basis is built by introducing adequate matching conditions. Although wavelet technology is still in a relatively early stage of development, systems of wavelets are assuming a role of increasing importance in signal and image processing, replacing with remarkable success in several situations the traditional mathematical tools based on Fourier analysis. In 1910, Haar [8] de®ned a set of piecewise constant orthogonal functions forming a complete series with the properties that characterize wavelet systems. In fact, Haar functions are wavelets, although their localisation in Fourier space is bad. The Haar functions have close similarities with the complete set of orthonormal functions introduced by Walsh [9] in 1923. The implementation of the hybrid-mixed stress model using a Walsh approximation basis is discussed in detail in [3,10]. It was only recently that the ideas behind the multiresolution analysis [11,12] theory provided a general framework for the development of wavelet systems. Using wavelets it is possible to represent a given signal or image at different levels of resolution and to compute at marginal cost the information necessary to switch from a coarse level of resolution to higher resolution representations. The main advantage of such systems when compared with the traditional trigonometric functions is that they present a local behaviour not only in the frequency domain, but also in the space/time domain [13]. The development of systems of wavelets is a very active ®eld of research. Important results have been obtained recently both in terms of theoretical development and practical applications. Probably the most successful applications of wavelets lie in image processing and data compression [11]. Some encouraging results are also being reported in numerical analysis [14]. There are several types of wavelets [12,15]. Apart from Haar wavelets, the only system that combines orthogonality and compact support in space domain, meaning that the function is de®ned over a limited and closed interval on the real line, is the one introduced by Daubechies [16]. The lack of regularity and symmetry is the price paid to ensure both properties. Moreover, and as it is shown below, the full use of the orthogonality of wavelet systems in the context of the ®nite element method is still dependent on the implementation of wavelets de®ned on limited and closed intervals [17]. The structure of the paper is as follows. The fundamental equations of elastostatics are recalled ®rst to establish the notation and the terminology in use. The equations of the stress model of the hybridmixed ®nite element formulation are stated next. Only the essential aspects are recalled here, as their derivation, the energy statements they are associated with and the conditions for the existence and uniqueness of the solutions they produce are presented in detail somewhere else [18±20]. The following section addresses the de®nition of wavelet systems, identi®es the techniques used to generate these functions, which do not have analytical expressions, and identi®es their main properties. After describing in detail how wavelets are used to implement the hybrid-mixed stress model, the paper ends with a report on a set of standard tests on convergence in energy, stresses and displacements and sensitivity to mesh distortion and incompressibility. Also reported are the information needed to characterize the structure and the computer time spent to set up and solve the ®nite element governing system.

2. Fundamental relations Let V denote a domain with boundary C. Let Cr and Cu be the complementary regions of the enveloping surface where the tractions tc and the displacements uc are prescribed, respectively (C ˆ Cr [ Cu ; Cr \ Cu ˆ 0). The fundamental relations governing the elastostatic response of plane structures may be summarized as follows:

L.M.S. Castro, J.A. Teixeira de Freitas / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3977±3998

3979

Dr ‡ b ˆ 0 in V ; Nr ˆ tc on Cr ;

…1† …2†

e ˆ D u in V ; u ˆ uc on Cu ;

…3† …4†

ee ˆ fr ‡ eh :

…5†

Vector ufbg lists the components of the displacement {body force} ®eld, and vector rfeg collects the independent components of the stress {strain} tensor. When a geometrically linear behaviour is assumed, the equilibrium (1) and compatibility (3) conditions represent conjugate transformations and the di€erential operators D and D are adjoint and linear. In the static boundary condition (2), matrix N has for entries the components of the unit outward normal associated with the operators listed in the equilibrium di€erential operator, D. The elastic constitutive relation is de®ned by (5), where f is the symmetric ¯exibility matrix of elastic constants characterizing a linear, reciprocal elastic law. Vector eh represents a generalized thermal strain vector.

3. Finite element approximation In the hybrid-mixed stress model used here, the components of the stress and displacement ®elds are independently approximated as follows: r ˆ Sv X

in V ;

…6†

u ˆ U v qv

in V :

…7†

Matrices Sv and Uv collect the approximation functions associated with the generalized stresses and generalized displacements X and qv , respectively. The displacement ®eld on the static boundary, which includes not only the part of the enveloping surface whereon tractions are prescribed but also the inter-element boundaries, is also approximated independently as follows: u ˆ U c qc

on Cr :

…8†

Matrix Uc lists the approximation functions associated with the generalized boundary displacements, qc . The dual transformations of approximations (6)±(8) are used to de®ne the following generalized strains, body forces and boundary tractions: Z e ˆ Stv e dV ; …9† Z …10† Qv ˆ Utv b dV ; Z …11† Qc ˆ Utc tc dCr : Convergence of this hybrid-mixed stress model depends essentially on the completeness of the approximation bases (6)±(8) and on the invariance of the inner product in the ®nite element mapping they imply. Completeness of the approximation bases is ensured by the use of systems of wavelets [13,16,21]. Invariance of the inner product is implied by de®nitions (9)±(11) which ensure that pairs of discrete, dual variables fX; eg, fQv ; qv g and fQc ; qc g dissipate the same energy as the ®elds they represent. Static-kinematic duality and elastic reciprocity in the hybrid-mixed ®nite element model are consequent upon these energy equivalence conditions.

3980

L.M.S. Castro, J.A. Teixeira de Freitas / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3977±3998

4. Finite element equations Di€erent techniques can 2 F Av …a† 6 At 0 …b† 4 v t Ac 0 …c†

be used to derive the resulting elementary governing system: 3 32 3 2 Ac X ec eh 6 7 6 7 0 7 54 qv 5 ˆ 4 Qv 5 0 qc Qc

…12†

Equation (12a) can be interpreted as the Galerkin Sv -weighted residual enforcement of the compatibility condition (3) Z Stv …e D u† dV ˆ 0: The equation above is integrated by parts to preserve static-kinematic duality and the elasticity and kinematic boundary conditions (5) and (4) are directly enforced: Z Z Z Z …N Sv †t u dCr …N Sv †t uc dCu ˆ 0: Stv …fr ‡ eh † dV ‡ …D Sv †t u dV The following expressions are found for the ¯exibility and compatibility matrices associated with the ®nite element approximations (6)±(8), and for the generalized strains associated with the prescribed displacements and residual strains: Z F ˆ Stv f Sv dV ; …13† Z t …14† Av ˆ …DSv † Uv dV ; Z …15† Ac ˆ …NSv †t Uc dCr ; Z …16† ec ˆ …NSv †t uc dCu ; Z eh ˆ Stv eh dV : …17† Similarly, equations (12b) and (12c) can be interpreted as the Galerkin Uv - and Uc -weighted residual enforcement of equilibrium conditions (1) and (2), respectively: Z Z Utv …Dr ‡ b† dV ˆ 0; Utc …Nr tc † dCr ˆ 0: The de®nitions given above for the equilibrium matrices and vectors present in equations (12b) and (12c) are recovered using the stress approximation (6). The global governing system is assembled requiring adjacent elements to share the same boundary displacement approximation (8). The resulting linear system has a structure similar to the one presented in (12). Except for the generalized boundary displacements, which can be shared at most by two neighbouring elements, all the remaining variables are strictly element-dependent. The assembling of the global governing system is straightforward and based on a direct allocation procedure [3,4], without the summations typical of the conventional displacement ®nite element formulation. When linearly independent stress approximation functions are used, the ¯exibility matrix is positivede®nite and system (12) can be condensed to form: K qc ˆ Qc ;

…18†

K ˆ At F 1 A ; Qc ˆ Qc At F 1 Qv ;

…19†



…20†

L.M.S. Castro, J.A. Teixeira de Freitas / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3977±3998

 F ˆ

F

Av

Atv

0



 A ˆ

Ac



0

Qv ˆ



eh

ec Qv

3981

 :

…21†

The condensed governing system (18) can be obtained either at elementary or at a global level. When the Schur complement matrix (19) is obtained at element level, the subsequent assembly procedure is similar to the one used in the conventional displacement formulation of the ®nite element method [3].

5. Wavelet bases 5.1. De®nition of a wavelet system The de®nition of a system of wavelets is based on two functions, the scaling function, /…x†, and the primary wavelet, w…x† [13,21]. In general, these functions do not have analytical de®nitions. Their generation is based on recursive expressions. The scaling function may be de®ned as follows, where ak are the ®lter coef®cients: X /…x† ˆ ak /…2x k†: …22† k2Z

Only a ®nite number of ®lter coecients may be non-zero to ensure compact support. Daubechies' wavelet systems are partitioned into families according to their wavelet number, N, which is equal to one-half of the (even) number of non-zero coef®cients present in de®nition (22). It can be proved [13] that N 1 also corresponds to the maximum degree of the polynomial that can be exactly represented as a linear combination of the scaling function and all its integer translates, /…x k†. Thus, the scaling function for the family with wavelet number N is de®ned by: N /…x†

ˆ

2N X1

ak N /…2x

k†:

…23†

kˆ0

This function has non-zero values only in the closed interval ‰0; 2N 1Š, known as the support of the scaling function. The primary wavelet has the same support and is de®ned as follows: N w…x†

ˆ

2N X1

bk N /…2x

k†;

…24†

kˆ0 k

bk ˆ … 1† a2N

1 k:

…25†

The de®nition of parameters bk is obtained imposing orthogonality between functions /…x† and w…x†. The numerical values of coecients ak are computed using a set of basic conditions derived from the imposition of several constraints involving the properties and behaviour of the set of functions to be generated [13]. Daubechies [16] obtains and presents the values of the wavelet coecients for all the families with N ranging from 2 to 10. Fig. 1 shows the scaling function and the primary wavelet for some of those families. The complete wavelet system is generated by taking translations and dilations on both scaling function and primary wavelet [15,16,21] to yield: N /j;k …x†

ˆ 2j=2 N /…2j x

k†;

…26†

N wj;k …x†

ˆ 2j=2 N w…2j x

k†:

…27†

The dilation and translation operations are illustrated in Fig. 2 for the case of the scaling function 4 /…x†. The support of each of these functions is de®ned by:

3982

L.M.S. Castro, J.A. Teixeira de Freitas / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3977±3998

Fig. 1. Scaling functions and primary wavelets for the families with N ˆ 2; 4; 10.

supp N /j;k …x† ˆ supp N wj;k …x† ˆ ‰k  2 j ; …2N

1 ‡ k†  2 j Š:

…28†

Two-dimensional wavelet systems may be de®ned by products of one-dimensional scaling functions and primary wavelets. The scaling function is then de®ned as the product of the corresponding onedimensional functions. Three di€erent primary wavelets may be obtained in 2D domains. Their de®nitions are: W…1† …x; y† ˆ /…x† w…y†; W…2† …x; y† ˆ w…x† /…y†; W…3† …x; y† ˆ w…x† w…y†:

L.M.S. Castro, J.A. Teixeira de Freitas / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3977±3998

3983

Fig. 2. Illustration of the dilation and translation operations for the case of the functions 4 /0;1 …x† and 4 /1;0 …x†.

Fig. 3. Two-dimensional scaling function and primary wavelets for N ˆ 2.

The plots presented in Fig. 3 show the four two-dimensional functions de®ned above for the family with wavelet number N ˆ 2. The translated and dilated two-dimensional wavelets may be de®ned as follows: Uj;…k1;k2† …x; y† ˆ 2j U…2j x …l† Wj;…k1;k2† …x; y†

ˆ 2j W…l† …2j x

k1 ; 2j y k 1 ; 2j y

k2 † ˆ /j;k1 …x† /j;k2 …y†; k2 †;

l ˆ 1; 2; 3:

3984

L.M.S. Castro, J.A. Teixeira de Freitas / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3977±3998

5.2. Generation of wavelets As in general scaling functions do not have analytical de®nitions, they have to be generated recursively directly from the dilation equation (23). It can be shown [21] that /…k† ˆ 0 when k < 0 or k > 2N 1. Using this property and the dilation equation (23) it is possible to write, LU ˆ U

with Lij ˆ a2i j ;

…29†

where vector U collects the values that the scaling function takes in every integer 0 < k < 2N 1. The generation of function /…x† on those integers is then equivalent to the computation of the eigenvector corresponding to the eigenvalue 1 of matrix L. As the solution is not unique, an additional normalization is necessary. It is usual [21,22] to impose that: 2N X1

/…k† ˆ 1:

…30†

kˆ0

Once the values that the scaling function takes in every integer are obtained, successive application of the dilation equation (23) will lead to the generation of /…x†, ®rst on the half-integers, then on the quarterintegers, and so on. This procedure can be repeated as many times as necessary to ®nd the values of the scaling function at all dyadic points, k  2 n …k; n; 2 Z†, in the interval de®ned by the support of the function. A similar technique is used to determine the derivative of the scaling function: /0 …x† ˆ

2N X1

2ak /0 …2x

k†:

…31†

kˆ0

Following the procedure presented above, it is possible to de®ne 2LU0 ˆ U0 ;

…32†

where U0 collects the values of the derivative of the scaling function on every integer on the support interval. This vector corresponds to the eigenvector associated with the eigenvalue 1=2 of matrix L. An additional condition is again required. It is usual to assume [22] that: 2N X1

k/0 …k† ˆ 1:

…33†

kˆ0

5.3. Properties of orthonormal wavelet bases The following are the most relevant properties of wavelet systems. P1: The wavelets wj;k …x† (with j; k 2 Z) form an orthonormal basis for the L2 …R† space: Z ‡1 wj;k …x† wm;n …x† dx ˆ djm dkn …j; k; m; n 2 Z†: 1

…34†

P2: The functions /j;k …x† (with k 2 Z) form an orthonormal basis of subspace Vj  L2 … R†. The set of nested subspaces Vj constitutes a multiresolution analysis meaning that


 V 2  V 1  V0  V1  V2 


and \ [ Vj ˆ 0; Vj ˆ L2 …R†: j2Z

j2Z

It is possible to verify [15] that Vj ! L2 …R† when j ! 1 and, in consequence of property P2: Z

‡1 1

/j;k …x† /j;m …x† dx ˆ dkm

…j; k; m 2 Z†:

…35†

L.M.S. Castro, J.A. Teixeira de Freitas / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3977±3998

P3: For a given value of dilation parameter j, the functions /j;k …x† and wj;m …x† are orthogonal: Z ‡1 /j;k …x† wj;m …x† dx ˆ 0 …j; k; m 2 Z†: 1

3985

…36†

P4: Every polynomial with a degree less than or equal to N 1 may be exactly represented as a linear combination of the scaling function N /…x† and all its integer translates, N /…x k†, with k 2 Z. 5.4. Approximation of functions Due to property P1, a given function f …x† may be represented on a wavelet basis using the following expansion: XX cj;k wj;k …x†; …37† f …x† ˆ j2Z

k2Z

Z cj;k ˆ

f …x† wj;k …x† dx:

supp N wj;k …x†

…38†

As the wavelets have compact support, when expanding the function f …x† on the interval ‰a; bŠ it is necessary to use the translations and dilations of the primary wavelet that present non-vanishing values inside the domain under consideration. The support de®ned in (28) may be used to identify the functions to use. For a given value of parameter j, the translations to take into account are de®ned by the conditions: k 6 b 2j

1;

k P a 2j ‡ 2…1

N †:

…39†

Expansion (37) involves the use of positive and negative dilations of the primary wavelet, w…x†. The second procedure to expand a function on wavelet basis is the most widely used. It is based on the scaling function and the positive dilations of the primary wavelet. It is known as canonical expansion and may be written as follows, X XX f …x† ˆ dk /…x k† ‡ cj;k wj;k …x†; dk ˆ

j2Z ‡ k2Z

k2Z

Z

supp N /0;k …x†

f …x† /…x

k† dx;

Z cj;k ˆ

supp N wj;k …x†

…40†

f …x† wj;k …x† dx …j; k 2 Z; j P 0†:

The third possibility to represent a given function on a wavelet basis consists in computing the projection of f …x† in the subspace Vj using the basis functions de®ned by the translations of the scaling function with dilation degree j: X dj;k /j;k …x†; j; k 2 Z; …41† f …x† ' Pj f …x† ˆ Z dj;k ˆ

suppN /j;k …x†

k2Z

f …x† /j;k …x† dx:

…42†

6. Implementation of the hybrid-mixed stress model 6.1. Selection of the approximation basis The ®rst problem that arises in the implementation of ®nite element models using wavelet systems is the selection of the functions used in the approximation basis over the ®nite element domain.

3986

L.M.S. Castro, J.A. Teixeira de Freitas / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3977±3998

Eq. (39) de®nes all functions that present non-vanishing values in a limited and closed interval ‰a; bŠ. Besides the functions that are completely de®ned inside this interval, the question arises on how to treat the remaining functions with supports lying partly outside the ®nite element domain. The simplest option consists in building the basis on all functions with non-trivial values in the interval ‰a; bŠ. Consequent upon this option, followed here, orthonormality is destroyed for all functions lying partly outside this interval. Moreover, some of the properties stated before are also lost. Numerical stability problems may occur, a€ecting especially the regions near the boundaries. These problems appear when the truncated functions present small values in the part of the support contained within the domain of the element. To minimize this e€ect, it is assumed here that all the truncated functions considered in the approximation must have a quadratic norm (computed on the interval under consideration) greater than a given tolerance, tol. Unfortunately the value of the tolerance that avoids the numerical instability is problem-dependent. In general, it is assumed that tol ˆ 10 3 . To overcome most of the drawbacks that this truncation may produce, the notion of wavelets on the interval [17] was recently introduced. The functions that intersect any of the edges of the element are so modi®ed as to de®ne a complete system of wavelets on the given interval with exactly the same properties discussed here for the systems de®ned over the real line. The use of these new systems of wavelets as approximation bases of the hybrid-mixed stress ®nite element formulation described here is expected to be reported soon. 6.2. De®nition of the ®nite element approximation In two-dimensional applications, the basic approximations (6) and (7) for the stress and displacement ®elds are written as follows, rˆ uˆ

nr X nr X mˆ1 nˆ1 nu X nu X mˆ1 nˆ1

I3 /j;m …x†/j;n …y† Xmn ;

…43†

I2 /j;m …x†/j;n …y† qvmn ;

…44†

where In is the n  n identity matrix and /j;k …s† represents the kth scaling functions along the direction s and with dilation parameter j. The number of stress and displacement degrees of freedom is av ˆ 3  n2r and bv ˆ 2  n2u , respectively. The weighting vectors are de®ned as Xtmn ˆ fXxx Xyy Xxy gmn and qtvmn ˆ fqvx qvy gmn . The boundary displacement approximation (8) is written in a similar form, uˆ

nc X nˆ1

Ik /j;n …s† qcn ;

…45†

where s is the coordinate system of the boundary whereon either both (k ˆ 2) or one (k ˆ 1) of the displacement components are approximated, depending on the boundary conditions of the element in the mesh. The dimension of the boundary displacement vector qc is bc ˆ nd  nc , where nd de®nes the number of displacement components approximated on the boundary of the element (nd ˆ 2  ns for an element with ns sides whereon both displacement components are unknown). A di€erent level of resolution can be selected for every approximation. To avoid spurious modes it is usual to approximate the stress ®eld using translations of scaling functions with a degree of re®nement j and both the displacements in the domain and on the boundary with scaling functions with dilation parameter j 1. 6.3. Computation of the structural operators Meshes built with non-rectangular elements are easily implemented using the geometric mapping typical of isoparametric elements. Approximations (6)±(8) are now de®ned on the coordinate system of the (square) master element, which is mapped into the desired shape by transformation x ˆ Nng c, where c is the nodal coordinate vector in the global system of reference and matrix Nng lists the shape functions [23].

L.M.S. Castro, J.A. Teixeira de Freitas / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3977±3998

3987

The matrices and vectors in system (12) present the general form (46), where dS is either an area or length element, J the Jacobian matrix of the element mapping and W the product of scaling functions and their derivatives: Z I ˆ W jJj dS: …46† Analytical solutions for the above integral can be obtained when the Jacobian is a constant and W derives from scaling functions intersecting within the domain of integration. Otherwise, numerical integration is used. Newton±C^ otes quadrature rules are applied because the wavelets are de®ned only at dyadic points. The irregularity of the functions requires the use of a large number of control points. 6.4. Solution of the governing system Di€erent strategies are used to solve the alternative descriptions (12) and (18) of the governing system. The large sparse non-condensed system (12) is solved with the Harwell Subroutine Library code MA47 [24] based on the multi-frontal solution technique developed by Du€ and Reid [25] with a pivot sequence derived from the Markowitz criterion [26] to minimize ®ll-in. The condensed system (18) is used when the ®nite element mesh is built on a signi®cant number of geometrically and mechanically identical elements. The Schur complement (19) is formed and its computation is optimized at element level. To minimize the storage requirements and to ensure the eciency of the process, it is essential to avoid the explicit computation of the inverse of matrix F in the de®nition of the triple product in (19). Column j of F 1 A , ckA j , is computed through the solution of the following linear system, F ckA j ˆ colj fA g

…47†

using the factorisation of F obtained by a direct solution procedure based on an algorithm initially introduced by Gustavson [27]. The assembled condensed governing system is then solved using the MA47 package. 6.5. Post-processing operations When wavelet systems are used, it is possible to use fast transform algorithms [28] to decompose or to recompose a given signal or function. However, such techniques could not be applied in the models presented here due to the truncation procedure used. In the numerical models presented here, post-processing operations include the computation of the stress and displacement ®elds in a mesh of 64  64 dyadic points in every element of the mesh using the direct approximations (6) and (7). It also includes the computation of the displacements on 64 dyadic points on every edge of the boundary using the direct approximation (8).

7. Numerical applications The clamped square plate subjected to a uniform transverse top load and the thick cylinder subjected to a uniform internal pressure shown in Figs. 4 and 6 are the two basic problems that support the numerical illustrations reported below. The ®nite element meshes used in these tests are shown in Figs. 5 and 7. The characterisation of the alternative solving systems (12) or (18) in terms of dimension and sparsity and the computer times spent in the pre-processing solution and post-processing phases of the analyses are the aspects addressed ®rst. The second set of results illustrates the convergence of the solution in terms of energy, displacement and stress estimates, and its (in)sensitivity to mesh distortion and incompressibility. No form of smoothing is used in the presentation of the ®nite element solutions. The stress and displacement estimates given below are directly computed from approximations (6)±(8).

3988

L.M.S. Castro, J.A. Teixeira de Freitas / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3977±3998

Fig. 4. Square cantilever plate (m ˆ 0:3).

Mesh A

Mesh B

Mesh C

Fig. 5. Meshes used in the analysis of the square plate.

Fig. 6. Thick cylinder under internal pressure.

7.1. Dimension and sparsity of the solving system The dimension and sparsity of the alternative solving systems (12) and (18) for the hybrid-mixed stress element are illustrated with the cantilever plate problem for the three meshes shown in Fig. 5. Two families of wavelets, N ˆ 5 and N ˆ 10, are used to approximate both the stress and displacement ®elds. The degree of approximation is determined by the dilation parameter, j, of the scaling functions.

L.M.S. Castro, J.A. Teixeira de Freitas / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3977±3998

3989

Fig. 7. Meshes used in the analysis of the thick cylinder.

The number of scaling functions used in the stress approximation, nr , and in the domain and boundary approximations for the displacements, nu and nc , respectively, is presented in Table 1 for mesh A on different dilation parameters, 1 6 j 6 3. The values summarized in Table 2 refer to the non-condensed system (12): ndof and nnz de®ne the total number of degrees of freedom and of non-zero coecients of the symmetric system and g is the sparsity index of the system, de®ned as the ratio between the number of zero coecients and the total number of entries in the governing system. The sparsity index is seen as relatively low for a model based on orthogonal approximation bases. It has been stated above that this is caused by the destruction of the orthogonality of Daubechies' wavelet systems implemented on closed intervals. Parameter g in the same table de®nes the sparsity index that could be obtained if orthogonality was preserved. Orthogonality of the wavelet systems is also destroyed when they are implemented on non-rectangular elements. Consequently, when the meshes shown in Fig. 8 are used, the sparsity indices are of the order of magnitude of those given in Table 2 for Mesh B. Table 1 Number of scaling functions involved in the analysis j…Nˆ5†

1

2

3

j…N ˆ10†

1

2

3

nr nu ˆ nc

6 5

8 6

12 8

nr nu ˆ nc

8 7

10 8

14 10

Table 2 Non-condensed governing system j…Nˆ5†

1

2

3

j…N ˆ10†

1

2

3

Mesh A ndof nnz g g

188 7854 0.5586 0.8320

300 21 162 0.5319 0.8574

608 56 170 0.6973 0.9017

ndof nnz g g

332 25 344 0.5419 0.8205

476 54 910 0.5166 0.8364

848 180 080 0.4999 0.8644

Mesh B nndof nnz g g

732 31 896 0.8818 0.9549

1176 85 016 0.8776 0.9624

2400 219 448 0.9241 0.9745

nndof nnz g g

1300 102 720 0.8789 0.9534

1872 221 560 0.8739 0.9578

3352 717 296 0.8725 0.9655

Mesh C nndof nnz g g

2888 128 544 0.9694 0.9833

4656 343 200 0.9685 0.9904

9536 907 264 0.9801 0.9936

ndof nnz g g

5144 413 568 0.9689 0.9882

7424 890 080 0.9678 0.9894

13 328 2 901 824 0.9674 0.9914

3990

L.M.S. Castro, J.A. Teixeira de Freitas / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3977±3998

Mesh B1

3

1

Mesh B2

4

3

4

2

1

2

Fig. 8. Meshes with non-rectangular elements.

The information summarized in Table 2 also stresses the relatively large dimension of the solving systems of hybrid-mixed stress elements, which should be handled and solved using appropriate large sparse matrix algorithms [26]. An alternative technique consists in condensing system (12) into form (18). The characteristics of the condensed governing system associated with each discretisation are presented in Table 3. There is a clear reduction in the dimension of the ®nal global governing system, Nsys . However, there is also a strong decrease in the sparsity index. 7.2. Computer times All tests presented here are run on an Intel Pentium II machine running Linux. CPU times are measured in seconds. In the results summarized in Table 4, Ttot de®nes the total solution time, and Tsol is the time spent in the solution of the condensed governing system (18). Tpre is the time spent from start-up to setting up the solving system and Tpost de®nes the time spent in the post-processing operations. 7.3. Sensitivity to mesh distortion To illustrate the e€ect of mesh distortion, the square cantilever plate is analysed using the 4-element mesh represented in Fig. 9. The distortion of the mesh is de®ned by parameter e, ranging from 1 to 7. Three di€erent families of wavelets, with N ˆ 6; 8 and 10 are considered. The stress ®eld in the domain is represented at resolution level j ˆ 1 and both displacement ®elds, in the domain and on the static boundary, are independently modelled at resolution j ˆ 0. The graph in Fig. 10 shows, for each tested discretisation, the values obtained for the transverse displacement at point A, wA , normalized to the values obtained with the undistorted mesh, wA …e ˆ 4†. The Table 3 Condensed governing systems j…Nˆ5†

1

2

3

j…Nˆ10†

1

2

3

Mesh A ndof Nsys g

188 30 0.0000

300 36 0.0000

608 48 0.0000

ndof Nsys g

332 42 0.0000

476 48 0.0000

848 60 0.0000

Mesh B ndof Nsys g

732 100 0.5400

1176 120 0.5400

2400 160 0.5400

ndof Nsys g

1300 140 0.5400

1872 160 0.5400

3352 200 0.5400

Mesh C ndof Nsys g

2888 360 0.8426

4656 432 0.8426

9536 576 0.8426

ndof Nsys g

5144 504 0.8426

7424 576 0.8426

13 328 720 0.8426

L.M.S. Castro, J.A. Teixeira de Freitas / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3977±3998

3991

Table 4 CPU times involved in the analysis j…Nˆ5†

1

2

3

j…N ˆ10†

1

2

3

Mesh A Tpre Tsol Tpost Ttot

1.74 3.75 0.03 5.52

1.76 4.12 0.07 5.95

2.18 10.09 0.12 12.39

Tpre Tsol Tpost Ttot

5.01 12.07 0.07 17.15

5.29 15.30 0.11 20.70

5.37 31.19 0.18 36.74

Mesh B Tpre Tsol Tpost Ttot

1.74 3.74 0.16 5.64

1.75 4.12 0.26 6.13

2.19 10.33 0.50 13.02

Tpre Tsol Tpost Ttot

5.02 12.12 0.27 17.71

5.32 15.32 0.40 21.04

5.35 31.29 0.73 37.37

Mesh C Tpre Tsol Tpost Ttot

1.75 3.77 0.64 6.16

1.81 4.32 1.00 7.13

2.22 11.22 1.97 15.41

Tpre Tsol Tpost Ttot

5.03 12.27 1.10 18.40

5.31 15.47 1.56 22.34

5.40 32.95 2.88 41.23

y

3 4

a −α

α =ea/8

1 2

α

A x

a −α

α

Fig. 9. Mesh used to study the e€ect of mesh distortion (E ˆ 1500, m ˆ 0:25).

wA (e) / wA (4) 1,030 1,025

N=6

1,020 N=8 1,015 N=10

1,010 1,005 1,000 0,995 0,990

e

0,985 0

1

2

3

4

5

6

Fig. 10. Displacement sensitivity to mesh distortion (square plate).

7

3992

L.M.S. Castro, J.A. Teixeira de Freitas / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3977±3998

values plotted are directly obtained from the domain approximation (7). It is possible to verify that the solution is quite insensitive to mesh distortion. Moreover, the sequence of results presented in Fig. 10 shows that the richer the ®nite element approximation is, the less sensitive are the estimates to mesh distortion. 7.4. Sensitivity to incompressibility The thick cylinder (re ˆ 2; ri ˆ 10) subjected to a uniform internal pressure is used to illustrate the sensitivity to incompressibility of the hybrid-mixed stress model used here. The results obtained with the single element mesh shown in Fig. 7 are presented in Fig. 11. The graph represents the variation of the strain energy with the Poisson ratio, normalized to the value obtained for m ˆ 0:4. Parameter N 9 represents the number of digits in the extension to the incompressibility limit m ˆ 0:5. For instance, N 9 ˆ 4 identi®es the test on Poisson ratio m ˆ 0:49999. 7.5. Convergence in energy When wavelet systems are used as approximation functions, two alternative p-re®nement strategies may be adopted. The ®rst consists in always using the same degree of re®nement ± or dilation parameter degree, ± while increasing the family wavelet number, N, associated with the selected basis functions. In the second alternative, the wavelet family used in the analysis is always the same while the degree of re®nement is successively increased. Table 5 summarizes the data on the discretisations used in the analysis of the square plate. In this case, p-re®nement is implemented with scaling functions belonging to families with a higher wavelet number, N. Families with N ˆ 4; 5; 6 and 8 are tested. The stress ®elds are approximated with scaling functions with re®nement degree j ˆ 1 and the displacements, both in the domain and on the boundary, are modelled using functions with j ˆ 0.

Fig. 11. Strain energy sensitivity to incompressibility (thick cylinder).

Table 5 Degrees of freedom in the solution of the square clamped plate (j ˆ 1) N

4 5 6 8

Mesh A

Mesh B

Mesh C

av

bv

bc

av

bv

bc

av

bv

bc

75 108 147 243

32 50 72 128

24 30 36 48

300 432 588 972

128 200 288 512

80 100 120 160

1200 1728 2352 3888

512 800 1152 2048

288 360 432 576

L.M.S. Castro, J.A. Teixeira de Freitas / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3977±3998

3993

Table 6 Degrees of freedom in the solution of the square clamped plate (N ˆ 5) j

2 3 4

Mesh A

Mesh B

Mesh C

av

bv

bc

av

bv

bc

av

bv

bc

192 432 1200

72 128 288

36 48 72

768 1728 4800

288 512 1152

120 160 240

3072 6912 19 200

1152 2048 4608

432 576 864

In the second strategy for p-re®nement, the scaling functions are ranged from j ˆ 1 to j ˆ 4 for the same wavelet family N ˆ 5. The degrees of freedom involved in each test on the three alternative ®nite element meshes under analysis are de®ned in Table 6 for 2 6 j 6 4 (see Table 5 for j ˆ 1 and N ˆ 5). The graph in Fig. 12 shows the variation with the total number of degrees of freedom, Ndof ˆ av ‡ bv ‡ bc , of the strain energy, U, normalized to the reference value, U0 , obtained with the discretisation involving the largest number of degrees of freedom (Mesh C, N ˆ 5 and j ˆ 4). Solid lines correspond to the solutions obtained with the approximations based on the use of different families of Daubechies' wavelets (j ˆ 1; 4 6 N 6 8) and dashed lines represent solutions obtained with discretisations using wavelets with N ˆ 5 and at different levels of resolution (N ˆ 5; 1 6 j 6 4). The solution is consistently re®ned either by increasing the number of elements in the mesh, or by augmenting the wavelet number of the scaling functions or the degree of re®nement of the scaling functions. For a given ®nite element mesh, the two alternative forms of p-re®nement produce the same estimate for the strain energy for the same number of degrees of freedom. As it is typical of equilibrium elements applied to force-driven problems, the solution sequences shown in Fig. 12 are monotonic and bound from above the exact value of the strain energy. 7.6. Convergence in displacements The p- and h-re®nements of the displacement solutions found for the clamped square plate are presented in Fig. 13. The weaker approximations produce good estimates for the displacements except at the (singular

U / U0 1,04 (N,j)=(4,1)

1,03 (5,1)

1,02

mesh A

mesh B

mesh C

N = 5 , 1 ≤ j ≤4

mesh A

mesh B

mesh C

j = 1 , 4 ≤N ≤ 8

(6,1)

(5,2)

1,01

(8,1) (5,3)

(5,4)

1,00

log(Ndof) 0,99 2

2,5

3

3,5

4

Fig. 12. Convergence of the strain energy under p-re®nement.

4,5

3994

L.M.S. Castro, J.A. Teixeira de Freitas / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3977±3998 j=1

j=2

j=3

=

j=4

3.75 p a / E

Fig. 13. De¯ected shapes of the clamped square cantilever plate.

stress) corners of the clamped edge. In all other tests, the displacement solutions satisfy with suf®cient accuracy both the kinematic boundary and the inter-element displacement continuity conditions. The boundary displacement solutions presented in Fig. 13 are computed from the domain approximation (7). The thick cylinder problem is used next to illustrate how the solutions from this approximation converge to the solutions computed from the independent boundary displacement approximation (8). The thick cylinder, with re ˆ 2ri ˆ 10 and m ˆ 0:2, is solved for the wavelet family N ˆ 10 in both stresses and displacements and scaling functions with dilation numbers j ˆ 2 and j ˆ 1, respectively. The total number of degrees of freedom is 476 and 952 for the single- and two-element meshes shown in Fig. 7, respectively. The solutions obtained for the horizontal displacement ux at the points indicated in Fig. 7 using the independent domain (V) and boundary (C) approximations are presented in Table 7. These estimates are sensitive to the approximation of the geometry of the cylinder. The results show, however, that re®nement of the mesh induces the convergence of the independent domain and boundary approximations of the displacements. The values presented in Table 7 for the boundary displacement estimate represent the average of the value computed from each boundary connecting at the particular point under assessment. The values of ux …C† at points C and F are not presented because boundary displacement ®eld, ux is not directly approximated on edge CF . 7.7. Convergence in stresses The variation of the polar components of the stress tensor along the interior edge of the thick cylinder (line ABC and 0 6 h 6 p=2 in Fig. 7) is shown in Figs. 14±16 for the single- and two-element mesh solutions. Table 7 Domain and boundary displacement estimates (thick cylinder) Measure point Mesh A Mesh B Exact [29]

ux …V † ux …C† ux …V † ux …C† ux

A

B

C

D

E

F

9.299 8.850 9.214 9.198 9.33(3)

6.458 6.412 6.504 6.496 6.5996

0.058 ± 0.001 ± 0.0

6.416 6.023 6.397 6.396 6.66(6)

4.464 4.528 4.517 4.522 4.7140

0.038 ± 0.000 ± 0.000

L.M.S. Castro, J.A. Teixeira de Freitas / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3977±3998

-0,80

3995

σr

-0,85 -0,90

σ r (exact) = -1

-0,95 -1,00 -1,05

Mesh A

-1,10

Mesh B -1,15

θ (rad)

-1,20 0,00

0,20

0,40

0,60

0,80

1,00

1,20

1,40

1,60

Fig. 14. Variation of rr on the interior edge of the cylinder.

2,40

σθ

2,20

Mesh A Mesh B

2,00

1,80

1,60

σ θ (exact) = 1.6667

1,40

θ (rad) 1,20 0,00

0,20

0,40

0,60

0,80

1,00

1,20

1,40

1,60

Fig. 15. Variation of rh on the interior edge of the cylinder.

The p-re®nement used is the same as the one de®ned above for the displacement estimates. The effect of the approximation of the geometry of the cylinder is still visible. The enforcement of the static boundary condition through Eq. (2) is stronger in the two-element mesh test, which is close to capturing the exact solution, given by rr ˆ 1:0, rh ˆ 1:6667 and rrh ˆ 0:0 [29]. The square clamped plate is used to illustrate the ability of the wavelet approximation to model high gradients in the vicinity of singular stress points. The shear stress solutions obtained for the alternative p-re®nement procedures de®ned above are presented in Figs. 17 and 18. In these ®gures, each row is associated with the same h-re®nement (Meshes A, B and C, respectively) and each column with the same p-re®nement. As in all cases presented above, no smoothing is used in this representation obtained with the graphic package Janela [30]. All tests presented in Figs. 17 and 18 satisfy the boundary condition on the shear stress. The interelement stress continuity condition is also satis®ed with sucient precision in the solution obtained with Meshes B and C.

3996

L.M.S. Castro, J.A. Teixeira de Freitas / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3977±3998

0,15

σ rθ

0,10

σ rθ (exact) = 0

0,05 0,00 -0,05 -0,10

Mesh A -0,15

Mesh B -0,20

θ (rad)

-0,25 -0,30 0,00

0,20

0,40

0,60

0,80

1,00

1,20

1,40

1,60

Fig. 16. Variation of rrh on the interior edge of the cylinder.

N=5

N=6

N=8

N = 10

Mesh C

Mesh B

Mesh A

N=4

-1.3E+00

0.0E+00

Fig. 17. Shear stress ®eld in the cantilever plate; wavelets of di€erent families with j ˆ 1.

The stress oscillation detected in the tests based on the weaker p- or h-re®nement sets is very much affected by the singular ®elds developing at the corners of the clamped edge. These oscillations are associated with a Gibbs-type effect due to truncation in the approximation basis. As the results in Figs. 17 and

L.M.S. Castro, J.A. Teixeira de Freitas / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3977±3998 j=2

j=3

j=4

Mesh C

Mesh B

Mesh A

j=1

3997

-1.3E+00

0.0E+00

Fig. 18. Shear stress ®eld in the cantilever plate; wavelets with di€erent degrees of re®nement and N ˆ 5.

18 show, this local stress instability is solved by increasing the number of elements in the mesh (h-re®nement) or by increasing either the wavelet number N or the resolution parameter j (p-re®nement).

8. Closure The preliminary tests presented here justify a continued investment in the use of wavelets as approximation functions of hybrid-mixed ®nite element models. To exploit fully the orthogonality of the wavelet systems and their gain in sparsity and speed of computation of the ®nite element solving system, it is necessary to focus research on the implementation of wavelets de®ned on the interval [17]. These systems of functions are completely de®ned on a limited and closed interval but present exactly the same properties of the systems de®ned over the real line. The de®nition of a system of wavelets on the interval is expected to minimize the instabilities near the edges and corners of the domain detected in some circumstances. To conveniently exploit the hierarchical nature of the wavelet systems, it is also necessary to use in the approximations the canonical expansion involving the use of both translations of the scaling function and translations and dilations of primary wavelets. The approximations based on the exclusive use of scaling functions at different levels of resolution are not the most effective in terms of adaptivity. The models presented and discussed here result from the implementation of the hybrid-mixed stress model. Preliminary results on the use of wavelets as approximation functions in the alternative hybridmixed displacement model have already been obtained. The structural applications of wavelet systems are also being extended to the elastoplastic analysis of laminar structures [5].

3998

L.M.S. Castro, J.A. Teixeira de Freitas / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3977±3998

Acknowledgements This research work corresponds to a part of the activities of the Structural Analysis Group of Instituto de Engenharia de Estruturas, Territ orio e Construcß~ao, ICIST. It has been partially supported by Fundacß~ao para a Ci^encia e a Tecnologia as part of the research program PRAXIS/2/2.1/CEG/33/94. References [1] J.A. Teixeira de Freitas, J.P. Moitinho de Almeida, E.M.B. Ribeiro Pereira, Non-conventional formulations for the ®nite element method, Struct. Engnrg. Mech. 4 (1996) 655±678. [2] J.A. Teixeira de Freitas, J.P. Moitinho de Almeida, E.M.B. Ribeiro Pereira, Non-conventional formulations for the ®nite element method, in: S.N. Atlury, P.E. O'Donoghue (Eds.), Modelling and Simulation Based Engineering, Palmdale, 1998, pp. 254±259. [3] L.M.S. Castro, Wavelets e Series de Walsh em Elementos Finitos, Ph.D. Thesis, Technical University of Lisbon, 1996. [4] E.M.B.R. Pereira, J.A.T. Freitas, A mixed-hybrid ®nite element model based on orthogonal functions, Int. J. Numer. Methods Engrg. 39 (1996) 1295±1312. [5] L.M.S. Castro, J.A.T. Freitas, Hybrid-mixed ®nite element elastoplastic analysis based on Walsh and wavelet interpolation, in: M. Papadrakakis, G. Bugeda (Eds.), Advanced Computational Methods in Structural Mechanics, CIMNE, 1996, pp. 146±165. [6] C. Canuto, A. Tabacco, K. Urban, The wavelet element method. Part I: construction and analysis, Appl. Comp. Harmonic Anal. 6 (1999) 1±52. [7] C. Canuto, A. Tabacco, K. Urban, The wavelet element method. Part II: realization and additional features in 2D and 3D, Appl. Comp. Harmonic Anal. 8 (2000) 123±165. [8] A. Haar, Z ur theorie der orthogonalen funktionensysteme, Math. Annal. 69 (1910) 331±371. [9] J.L. Walsh, A closed set of normal orthogonal functions, Ann. J. Math. 55 (1923) 5±24. [10] L.M.S. Castro, J.A.T. Freitas, Implementation of Walsh-based hybrid-mixed stress elements, Comput. Struct. (submitted). [11] S.G. Mallat, A theory for multiresolution signal decomposition: the wavelet representation, IEEE Trans. Pattern Anal. Machine Intell. 11 (7) (1989) 674±693. [12] Y. Meyer, Ondelettes et Operateurs, vol. 1, ®rst ed., Hermann, Paris, 1990. [13] J.R. Williams, K. Amaratunga, Introduction to wavelets in engineering, Int. J. Numer. Methods Engrg. 37 (1994) 2365±2388. [14] S. Bertoluzza, G. Naldi, J.C. Ravel, Wavelet methods for the numerical solution of boundary value problems on the interval, in: C.K. Chui, L. Montefusco, L. Puccio (Eds.), Wavelets: Theory, Algorithms and Applications, Academic Press, New York, 1994, pp. 425±448. [15] C.K. Chui, An Introduction to Wavelets, Academic Press, New York, 1992. [16] I. Daubechies, Orthonormal bases of compactly supported wavelets, Commun. Pure Appl. Math. 41 (1988) 909±996. [17] A. Cohen, I. Daubechies, P. Vial, Wavelets on the interval and fast wavelet transforms, Appl. Comput. Harmonic Anal. 1 (1993) 54±81. [18] J.A. Teixeira de Freitas, Duality and symmetry in mixed integral methods of elastostatics, Int. J. Numer. Methods Engrg. 28 (1989) 1161±1179. [19] J.A. Teixeira de Freitas, Mixed and hybrid symmetric formulations for the boundary element method, Eur. J. Mech. A/Solids 9 (1) (1990) 1±20. [20] J.A. Teixeira de Freitas, E.M.B.R. Pereira, L.M.S. Castro, Consistent hybrid-mixed stress and displacement ®nite elements, Internal Report ICIST, Technical University of Lisbon, 2000. [21] G. Strang, Wavelets and dilation equations: a brief introduction, SIAM Rev. 32 (4) (1989) 614±627. [22] B.K. Alpert, Wavelets and other bases for fast numerical linear algebra, in: C.K. Chui (Ed.), Wavelets ± A Tutorial in Theory and Applications, Academic Press, New York, 1992, pp. 181±216. [23] O.C. Zienkiewicz, R.L. Taylor, The Finite Element Method, Volume 1 ± Basic Formulation and Linear Problems, fourth ed., McGraw-Hill, Berkshire, 1989. [24] Harwell Subroutine Library, Advanced Computing Department, Harwell Laboratory, Oxfordshire, 1990. [25] I.S. Du€, J.K. Reid, The multifrontal solution of inde®nite sparse symmetric linear systems, ACM Trans. Math. Softw. 9 (1983) 302±325. [26] I.S. Du€, A.M. Erisman, J.K. Reid, Direct methods for sparse matrices, Monographs on Numerical Analysis, Oxford Science Publications, Oxford, 1986. [27] F.G. Gustavson, Some basic techniques for solving sparse systems of linear equations, in: D.J. Rose, R.A. Willough (Eds.), Sparse Matrices and their Applications, IBM T.J. Watson Center, Yorktown Heights, New York, 1971, pp. 41±52. [28] G. Beylkin, R. Coifman, V. Rokhlin, Fast wavelet transforms and numerical algorithms I, Commun. Pure Appl. Math. 44 (1991) 141±183. [29] S.P. Timoshenko, J.N. Goodier, Theory of Elasticity, third ed., McGraw-Hill, New York, 1982. [30] J.P.B. Moitinho de Almeida, Janela: Uma interface gra®ca destinada a aplicacß~ao em problemas de Mec^anica Computacional, Internal Report, Technical University of Lisbon, 1992.