The generalized finite element method

Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

www.elsevier.com/locate/cma

Computational Mechanics Advances

The generalized ®nite element method T. Strouboulis a,b,*, K. Copps a,b,1, I. Babuska a,b b

a Department of Aerospace Engineering, Texas A&M University, College Station, TX 77843-3141, USA Texas Institute for Computational and Applied Mathematics, University of Texas at Austin, Austin, TX 78712, USA

Received 7 July 2000

Abstract This paper describes a pilot design and implementation of the generalized ®nite element method (GFEM), as a direct extension of the standard ®nite element method (SFEM, or FEM), which makes possible the accurate solution of engineering problems in complex domains which may be practically impossible to solve using the FEM. The development of the GFEM is illustrated for the Laplacian in two space dimensions in domains which may include several hundreds of voids, and/or cracks, for which the construction of meshes used by the FEM is practically impossible. The two main capabilities are: (1) It can construct the approximation using meshes which may overlap part, or all, of the domain boundary. (2) It can incorporate into the approximation handbook functions, which are known analytically, or are generated numerically, and approximate well the solution of the boundary value problem in the neighborhood of corner points, voids, cracks, etc. The main tool is a special integration algorithm, which we call the Fast Remeshing approach, which is robust and works for any domain with arbitrary complexity. The incorporation of the handbook functions into the GFEM is done by employing the partition of unity method (PUM). The presented formulations and implementations can be easily extended to the multimaterial medium where the voids are replaced by inclusions of various shapes and sizes, and to the case of the elasticity problem. This work can also be understood as a pilot study for the feasibility and demonstration of the capabilities of the GFEM, which is needed before analogous implementations are attempted in the three-dmensional and nonlinear cases, which are the cases of main interest for future work. Ó 2001 Elsevier Science B.V. All rights reserved. Keywords: Generalized ®nite element method; Hybrid; Meshless; Meshfree; Numerical integration; Adaptive quadrature; Handbook functions

1. Introduction The generalized ®nite element method (GFEM) was introduced in [1,2], and is based on the partition of unity method (PUM) which ®rst appeared in the work of Babuska et al. [3], and subsequently in [4±6], and the papers that originated from these works, e.g., [7±13]. Babuska, and Melenk [4,5,7±9] illustrated the development of the PUM for the Laplacian, the elasticity and the Helmholtz problems, and a general class of elliptic equations, with the main focus being the construction of an approximation which employs the local knowledge about the character of the solution into the approximation. In [4,5,7±9] the PUM was analyzed, and both a priori and a posteriori estimates were derived, and model computations were presented. Oden and Duarte [6,10±13] presented a meshless construction of the PUM called the hp cloud method (hpCM), which employs very general Partitions of Unity called ``clouds'', incorporated special functions which re¯ect the local character of the solution into the construction of the approximation, developed both a priori and a posteriori estimates for the hpCM, and discussed the implementation of the method using object oriented programming (OOP). A complete discussion of the references of various other

*

Corresponding author. E-mail address: [email protected] (T. Strouboulis). 1 Present address: Sandia National Laboratories, PO Box 5800, Albuquerque, NM 87185-0826, USA. 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 1 ) 0 0 1 8 8 - 8

4082

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

relevant engineering works, known under the name meshless methods, may be found in [6], and the recent works of Belytschko [14±18], which employ the PUM for the solution of problems of crack propagation. This paper follows the main ideas presented in [1,2] and focuses on the question of how the GFEM can be used to extend the capability of existing FE codes, to solve problems in domains with complex geometries. The GFEM, as it is understood in [1,2] and this paper, is an extension of the standard ®nite element method (SFEM, or FEM) that is endowed with the following two important capabilities: 1. Domain independent meshes: By this we mean the capability of the GFEM to build the approximation on meshes which may be partially or totally independent of the geometry of the problem domain. This capability of the GFEM is often characterized by the term meshless which could be misleading, because a mesh of elements, which we call the approximation mesh, is indeed used for the construction of the GFEM approximation. The information about the geometry of the domain enters into the formulation in terms of the approximation mesh, and element integration meshes which are used by a sophisticated numerical quadrature algorithm, and can be easily constructed in each element of the approximation mesh, even for very complex domains by employing a special re®nement algorithm. 2. Local enrichment of the approximation by special functions: By this we mean the capability of incorporating any special function, or functions, of interest into the construction of the approximation. An optimal choice of the special functions is the so-called handbook functions, which re¯ect the local character of the solution either in the interior of the domain, or next to a smooth boundary, or in the neighborhood of re-entrant corners and/or multi-material wedge vertices, or in the neighborhood of various inclusions, voids, micro-cracks, etc. The handbook functions are de®ned as exact solutions of handbook problems, which are formulated by using known local information about the differential equation, the local geometry of the boundaries, the boundary conditions, the loads, and the material data and can be determined a priori, i.e., prior to the computation of the approximate solution for the concrete case of interest. The handbook functions are incorporated into the approximation by the PUM, i.e., after they have been multiplied by the mapped bilinear ``hat''-functions de®ned on the approximation mesh employed by the GFEM, to obtain vertex handbook functions. As we will see below, the addition of a few properly selected vertex handbook functions to the basis of the GFEM, can lead to signi®cant improvements in the accuracy of the computed solution. A main contribution of Strouboulis et al. [2], and this work, is the development of the GFEM as an extension of the classical FEM that employs meshes which may not conform to the boundary of the problem domain. Similar ideas for solving problems using meshes which do not conform to the domain geometry, were pursued previously, e.g., by Liszka and Orkisz [19] in the context of ®nite di€erence methods, and more recently by Liszka [20]. The main tool of our approach is a robust numerical integration scheme based on special re®nements of the elements, which practically works well for any mesh and any domain. Another important contribution of Bab uska et al. and Strouboulis et al. [1,2] and this work is the implementation of the GFEM as a combination of the standard FEM and the PUM with the capability to incorporate into the approximation special handbook functions, e.g., corner, void, inclusion, functions, which signi®cantly enhance the accuracy of the GFEM solution. Finally, another important contribution is the development, implementation and validation of a robust linear equation solver for the semi-de®nite systems which result from the application of the GFEM to strongly elliptic partial di€erential equations, e.g., the Laplacian, the elasticity problem, etc. Following Section 1, in Section 2 we discuss the main characteristics of the GFEM, in Section 3 we describe the new numerical integration approach, and in Section 4 we give an illustration of the capabilities of the GFEM for solving boundary value problems for the Laplacian in polygonal domains with simple, or complicated geometry. We give several examples which underline the capability of the GFEM of solving problems by employing GFEM meshes which may overlap the boundary of the problem domain, and hence are simple to construct even for very complex domains which may include hundreds of voids and cracks. The examples also illustrate the e€ect of analytically known corner, crack, and void functions on the accuracy of the GFEM approximation. We also show that the a posteriori estimation and the adaptivity capabilities of the SFEM, can be directly extended to the GFEM. In Section 5 we give an outline of the procedure for the numerical generation and use of handbook solutions in a GFEM approximation, we describe how the implementation was validated, and we give various examples of using the numerically generated handbooks illustrating the potential of the GFEM for achieving high accuracy in problems set in

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4083

arbitrarily complex domains. In Section 6 we give an overview of the architecture of the GFEM computer code, and in Section 7 we discuss the conclusions and give recommendations for further work. 2. The main attributes of the GFEM As we have already mentioned in Section 1, the main capabilities of the GFEM are: 1. The capability of using meshes which are partially, or totally independent of geometry of the domain. 2. The capability of enriching the approximation by any special functions which are of interest. Let us ®rst address the type of meshes which can be employed by the GFEM. The FEM always employs one mesh, which we call the FE mesh, for the construction of the approximation, and the numerical integrations of the sti€ness, load entries, norms, etc. The FE mesh is most often constructed by subdividing the domain into a set of non-overlapping curvilinear triangular and/or quadrilateral subdomains, which are called elements and are required to satisfy various restrictions on the distortion and the connection between neighbors (see [21,22]). Examples of classical FE meshes are shown in Fig. 1(a) for a polygonal domain with straight edges, which we will call Domain I, and in Fig. 1(b), for a curvilinear polygonal domain which we will call Domain II. The GFEM employs two meshes, namely an approximation mesh, which is used for the construction of the basis of the approximation, and a corresponding integration mesh, which is constructed and separately each element of the approximation mesh and is for the numerical evaluation of all the integrals. A pictorial example of an approximation mesh and the corresponding integration mesh for Domain II is shown in Fig. 2. The approximation mesh need only satisfy the requirement that it has to cover completely the problem domain, e.g., as shown in Fig. 2(a), while the integration mesh is obtained by special re®nements of each element of the approximation mesh, which take into account the local geometry of the domain, as shown in Fig. 2(b), and the details shown in Figs. 2(c) and (d). Another example of an approximation mesh which can be used by the GFEM, and the corresponding integration mesh is shown in Fig. 3. Note that the GFEM integration mesh is constructed in one element at a time, using a re®nement algorithm, described in Section 3, which is robust and practically works for any domain, irrespective of its geometric complexity. We will classify the GFEM into three categories according to the relation between the employed GFEM approximation mesh and the geometry of the domain.

Fig. 1. Examples of meshes employed by the classical FEM: (a) A mesh for Domain I. (b) A mesh for Domain II. The above meshes were generated by two uniform nested re®nements of the coarse meshes shown in Figs. 21(a) and (b), respectively.

4084

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 2. An example of the two meshes employed in the construction of the GFEM for Domain II: (a) An approximation mesh obtained by employing uniform re®nements of a square which includes Domain II in its interior. (b) The corresponding integration mesh. (c) and (d) Details of the integration mesh.

GFEM I GFEM II

The approximation mesh is a classical FE mesh for the domain. Examples of GFEM I meshes for Domains I and II are shown in Fig. 1. The approximation mesh is a classical FE mesh for a modi®ed domain which includes the original domain in its interior, and is obtained from the original domain by deleting various parts of the boundary and adding new boundaries. Examples of GFEM II meshes for Domains I and II are shown in Figs. 3 and 4.

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

GFEM III

4085

The approximation mesh is a classical FE mesh for a domain which covers the original domain, and has no common boundaries with it. An example of a simple GFEM III mesh for Domain II is shown in Fig. 2(a).

Fig. 3. Another example of the two meshes employed in the construction of the GFEM: (a) An approximation mesh for Domain II, which is obtained as a coarse mesh of quadrilaterals with straight edges of the polygonal domains obtained by replacing the outer boundary of the domain, by C01 , a polygonal with straight edges which encloses the original boundary. (b) The corresponding integration mesh. (c) and (d) Details of the integration mesh.

4086

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

The second main attribute of the GFEM is the type of functions that can be incorporated into the basis of the approximation. Here, as in [1,2], we are mostly interested in the GFEM obtained by adding special functions to a bi-p FE basis de®ned on an approximation mesh of mapped quadrilaterals, by employing special functions which re¯ect the local character of the solution, e.g., corner functions, functions which satisfy the boundary conditions on a curved edge, crack, or void, etc. As we already mentioned in Section 1, the special functions are added to the basis of the approximation in a conforming way at a set of vertices after they are multiplied by mapped bilinear FE functions corresponding to the vertices. For example, a re-entrant corner

Fig. 4. Examples of GFEM II meshes. (a) A GFEM mesh for Domain I; (b) and (c) GFEM II meshes for Domain II.

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4087

function rkj fj …h†, may be incorporated into the approximation after it has been multiplied in by the vertex ``hat'' functions uX …x; y† for a set of vertices in the neighborhood of the re-entrant corner, e.g., the vertices marked by a triangle in Figs. 5 and 6, to obtain the functions …j†

uX ˆ uX …x; y†rkj fj …h†:

…2:1†

…j†

We will call uX , the jth vertex corner function for the vertex X corresponding to the corner located at the origin of the polar coordinate system …r; h†. The vertices to which the vertex corner functions are added around each corner are determined by specifying the parameter nlayers , which is understood as follows: · For GFEM I, nlayers ˆ 0, speci®es that vertex corner functions, only be used at the vertices at the re-entrant corners, while nlayers ˆ 1, speci®es that vertex corner functions be used at all vertices of the elements connected to a re-entrant corner, and nlayers ˆ 2, speci®es that the vertex corner functions be used at all the vertices of all the elements connected to the vertices obtained for nlayers ˆ 1, etc., as is, for example, shown in Fig. 5.

Fig. 5. An example of a GFEM I mesh with added re-entrant corner functions. This mesh was obtained by starting from the FE mesh shown in Fig. 21(b), by re®ning the elements with a vertex at a re-entrant corner to enforce the rule that any element intersecting a corner may not have vertex special functions from other corners. Triangular symbols indicate that one or more vertex corner functions are associated with that vertex.

4088

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 6. An example of a GFEM III mesh, with added re-entrant corner functions. This mesh was obtained by starting from the uniform mesh of squares shown in Fig. 2(a), by re®ning the elements in the neighborhood of the re-entrant corners to enforce the rule that any element intersecting a corner may not have vertex special functions originating from other corners. The vertices with a triangular symbol have one or more associated vertex corner functions.

· For GFEM II, and GFEM III, nlayers ˆ 0, speci®es that vertex corner functions be used at the four vertices of the elements which overlap a re-entrant corner, while nlayers ˆ 1, and higher, speci®es that the placement of the vertex corner functions continue as in GFEM I. An example is shown in Figs. 6 and 7. The corner functions are one example of analytically known handbook functions used to enrich the basis of the GFEM. Another example of handbook functions, which will be used in the numerical illustrations below, are the circular void functions w2k 1 ˆ R…zk ‡ bz k †; w2k ˆ I…zk ‡ cz k †; k ˆ 1; . . . ; pvoid ; …2:2† p   1, r and h are the polar coordinates associated with the circular void, and a, b and c where z ˆ reIh , I ˆ depend on the type of boundary conditions, e.g., homogeneous Dirichlet or Neumann applied on the boundary of the void. We will also employ the elliptical void functions obtained from the circular void functions by employing an appropriate mapping, and straight crack functions obtained from the elliptical void functions in the special case that the void is ¯at. The circular, elliptical void, and straight crack functions are added to the approximation after they have been multiplied by the vertex hat functions for a speci®ed set of vertices of the approximation mesh. Using the parameter nlayers in a similar way as in the case w0 ˆ a ln r;

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4089

Fig. 7. An example which illustrates the meaning of the parameter nlayers when vertex corner functions are added to the approximation in the neighborhood of the re-entrant corners.

Fig. 8. An example which illustrates the meaning of the parameter nlayers for elliptical voids, and straight crack functions.

4090

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

of the corner functions a set of vertices in the neighborhood of void or crack may be speci®ed, as is, for example, shown in Fig. 8. 3. The integration approach employed in the GFEM 3.1. Element integration algorithms based on the fast remeshing approach In [2] we employed a pixelation approach for the numerical integrations in the GFEM. The idea was to employ a nested subdivision of each element, which is suf®ciently re®ned across the domain boundaries for

Fig. 9. The two approaches for integrating approximation functions and their products by numerical quadrature on the region of domain containing multiple voids and straight cracks covered by a ®nite element mesh of regular squares. (a) The Pixelation approach (see [2]); (b) The Fast Remeshing approach.

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4091

the numerical integration, replacing all boundaries which are overlapped by the GFEM mesh, by a jagged pixel approximation, as is, for example, shown in Fig. 9(a). The Pixelation approach has the disadvantage of being computer memory intensive and slow because of the large number of subdivisions which are needed to control the error due to the approximation of the domain boundaries overlapped by the mesh by jagged polygonal lines. A more ef®cient approach can be obtained by modifying the Pixelation approach near the boundary to employ additional subcells that form a map to the exact geometry of the domain. We will call this modi®ed pixelation approach, the Fast Remeshing approach, because it results in the exact subdivision of the problem domain into a set of curvilinear subcells, as is, for example, shown in Fig. 9(b). Using the integration mesh constructed by the Fast Remeshing approach, the GFEM integrals are computed by employing embedded high-order quadrature rules of degrees 7, 9 and 13, on the cells, and on the subcells which are obtained by the subdivision of the cells intersecting the boundary. The mapping used in the subcells are constructed using the blending function method described in [23, pp. 107±111]. The Fast Remeshing quadrature algorithm can be summarized as follows: Algorithm 1 (Fast Remeshing quadrature). To compute the integral Z I‰fŠ ˆ f:

…3:1†

s

Begin Initial cell division

Assign subcells Initial estimate

Control Process cells

Update

Let ncells ˆ 1, xs;cell ˆ s. 1 Continue to divide all cells crossing a domain boundary until either each cell is completely inside or outside the domain, or the geometry of the domain boundary contained within each cell satis®es one case in the set of cell stopping criteria shown in Fig. 10. For each cell crossing the domain boundary, create a mesh of subcells whose topology is speci®c to the corresponding resulting geometry of the case of stopping criteria. Fig. 11 shows examples of these subcell mesh topologies. For all cells intersecting domain boundaries xks;cell \ X use the seventh degree embedded rule, or other suitable degree rule, in each master subcell and sum them to estimate the value of the integral Ixs;cell and the error Exs;cell over the region. k k For all other xs;cell in the element, use the seventh degree embedded rule in each k master cell and extrapolation to get Ixs;cell and the . P error Exs;cell k k Compute the estimate of the total integral I ˆ I s;cell . x P k Compute the estimate of the error E ˆ Exs;cell . k do while E=jIj2 > erel Find the maximum error in all cells, Emax ˆ maxk …Exs;cell †. k

For each cell that attains Emax , if the cell has subcells, delete the subcells, divide the cell into four new cells and assign a new subcell mesh topologies in those cells crossing the domain boundary; if the cell is completely inside the domain, divide the cell into four cells. Recompute value of global integral I and the error E. end do

The above algorithm is suitable for integrating in domains with arbitrarily complex polygonal boundaries with straight edges, as is, for example, shown in Fig. 12. It is also possible to construct a threedimensional analog of the Fast Remeshing approach, which is the most essential tool which is needed to enable the implementation of the GFEM in three dimensions. 3.2. Numerical quadrature for curved features Let us now give details of the use of the subcell mappings in order to get an exact representation of all domain boundaries, the computation of the boundary integrals, and the veri®cation of the implementation of all numerical integrations.

4092

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 10. The valid stopping criteria for the initial cell division step during the Fast Remeshing algorithm for numerical integration in an element. The domain always lies to the left of the domain edges oriented by the direction of the arrows. Each of these cases, and their permutations created by rotation, also has a corresponding subcell mesh topology, a few examples of which are shown in Fig. 11.

Let us consider the case that the approximation mesh for a curvilinear polygonal domain consists only of one square element which includes the problem domain in its interior, as shown in Fig. 13. Recall, that the Fast Remeshing algorithm, for an element s, ®rst partitions the element into a quadtree set of non-overlapping cells, which we denote xj , obtained from a nested subdivision of the element s, then divides the cells recursively until they satisfy one of a set of stopping criteria. For the numerical evaluation of the integrals on the cells which are completely inside, a standard quadrature rule on the square ^ is used, while for the cells that intersect the domain boundary, namely, xj \ oX 6ˆ ;, the master cell x …x † integration employs a subdivision of xj into a set of non-overlapping subcells, wk j , k ˆ 1; . . . ; nxj , such that

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4093

Fig. 11. The subcell meshes created inside the integration cells xs;cell for the polygonal regions overlapping the domain using the Fast k Remeshing approach. Oriented domain edges are shown with arrows and the subcell topologies are grouped according to the sides of the domain region intersected by the cell.

Fig. 12. A cell/subcell integration mesh created in the initial cell divide step, and subcell assign step, of the algorithm for the fast remesh approach of integration. The algorithm is suited for arbitrarily complex geometry of the problem domain boundary.

xj \ X ˆ

[ k

…x †

wk j :

…3:2†

If the part of the boundary oX inside the cell xj is a curved boundary segment, a mapping is constructed in each subcell by the blending function method (see [23]) in three steps, as depicted schematically in Fig. 14. Step 1: For the purposes of creating the subcell's nodes and edges, the boundary piece of oX cut and lying inside the cell xj is treated as collapsed into a straight segment s0 , as shown in Fig. 14(a).

4094

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 13. The quadtree of integration cells xj inside a single square element s covering the domain X, and the subcells wk inside the cells crossing the boundary oX. Edges of subcells crossing touching the boundary are shown as straight to emphasize that a special blending function mapping is used on these edges in order that the subcell mapping conforms exactly to the geometry of the domain.

Fig. 14. Schematic example of the steps used by the Fast Remeshing algorithm to create the subcell blending function mappings for a triangular region of domain intersected by a cell (a) Step 1; (b) Step 2 and (c) Step 3.

Step 2: The subnode lying on s0 is relocated to the nearest point on the domain boundary oX as shown in Fig. 14(b). The positions of the other subnodes may be adjusted to insure that all internal angles of the subcells remain less than 180°, to insure the invertibility of the Jacobian of the subcell transformation. Step 3: Using the blending function approach, new mappings are created for the subcells which originally ended on the straight edge s0 , as shown schematically in Fig. 14(c). Let us illustrate how we tested the implementation of the Fast Remeshing algorithm, ®rst for a simple example, up to step 1 of the subcell mapping, then in a more complex example, up to the ®nal step 3. We consider integrating the area on a square domain cut by an elliptical void, and meshed by a single element, as shown in Fig. 15(a). In this case the area missing due to the ellipse is pab, where a and b are the major and minor radii. To test the ®rst step of the subcell mapping, we treated the pieces of the boundary of the domain of Fig. 15 cut by integration cells as straight, i.e., without the blending function mapping, and we computed approximations to the area of the domain by employing a sequence of successively uniformly ®ner cells as in Fig. 15(b). In this manner, we expect the approximations to the area to converge to the exact

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4095

Fig. 15. Numerical convergence test of Step 1 of the Fast Remeshing integration algorithm. (a) A square domain cut by an elliptical void and the computed area A0 of the domain using the Step 1 of the Fast Remeshing algorithm, where A ˆ 3:6230088 is the exact area; (b) the sequence of successively ®ner uniform quadtrees of cells of size hcell ; (c) values of the relative error in the computed area.

area. Fig. 15(c) gives the values of the relative error in the computed area; it can be seen that this error is practically negligible even for very coarse cell sizes. For testing the third step of the subcell mapping, we computed the area of Domain II (from Section 2) by employing two types of meshes, the sequence of uniform GFEM I meshes shown in Fig. 16(a), and the uniform GFEM III meshes shown in Fig. 16(b). In the GFEM I meshes, we employed the two-dimensional Gauss rule, of degree 7, on each element, while in the GFEM III meshes we used the full implementation of the Fast Remeshing algorithm in each element, including step 3 of mapping the subcells, with a requested

Fig. 16. Sequences of meshes used to test the full implementation of the Fast Remeshing algorithm. (a) GFEM I, and (b) GFEM III.

4096

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Table 1 Convergence of the computed area on Domain II using Gauss rule of degree 7, and Fast Remeshing erel ˆ 0:1 on GFEM I and GFEM III meshes, respectively Element size, h

Gauss rule (GFEM I), A^

Fast Remeshing (GFEM III), A

 j A^ Aj A^

1 1/2 1/4 1/8 1/16

42.12511623023533 42.12511660277737 42.12511660373347 42.12511660373407 42.12511660373405

42.12511660456875 42.12511660386099 42.12511660373371 42.12511660373401 42.12511660373401

3  10 1  10 2  10 6  10 4  10

9 11 15 16 16

tolerance erel ˆ 0:1. Table 1 reports the values of the relative error in the computed areas. Further testing of the integration algorithm was done by computing the moments of inertia, of ®rst- and higher-order, etc. To further elucidate the role of the subcell mappings in integration, let us trace the contribution of one integration point to a single diagonal entry in an element sti€ness matrix for a typical GFEM basis. The contribution to the diagonal sti€ness matrix entry kjj for the special basis function uX wxj X coming from a single subcell wk is !   Z Z 2  o  2 o xX xX X X u …n; g†wj …x; y† u …n; g†wj …x; y† kjj jwk ˆ ‡ dx dy; …3:3† ox oy wk where uX is the hat function, wxj X is the jth special function associated with vertex X, and the ``physical'' coordinates of the domain are …x; y†; the local coordinate system of the feature giving rise to the special function as …x; y†; the master coordinates of the element as …n; g†; the master coordinates of the subcell as ^ g^†. For the ®rst bracketed term in the integrand, we have …n;  o X o X o u …n; g†wxj X …x; y† ˆ u …n; g†wxj X …x; y† ‡ uX …n; g† wxj X …x; y† …3:4a† ox ox ox !   on o og o o x o o y o x ˆ ‡ uX …n; g†wj X …x; y† ‡ uX …n; g† ‡ wxj X …x; y† ox on ox og ox ox ox o y …3:4b† and similarly for the second term, hence 0"    Z Z    on o og o ^ g^†; Qg …n; ^ g^† wxX Px …n; ^ g^†; Py…n; ^ g^† @ ‡ kjj jwk ˆ uX Qn …n; j ox on ox og ^k w 1 #2 !   ox o o   y o ^ g^†; Qg …n; ^ g^† ^ g^†; Py…n; ^ g^† ‡ uX Qn …n; g; ‡ ‰


Š2 AjJ j dn^d^ ‡ wxj X Px …n; ox ox ox o y

…3:5†

where ‰J Š ˆ

ox oy g on^ o^

ox oy o^ g on^

…3:6†

and Q is the mapping from the master subcell to the master element ^ k 7! si ; Q:w

n ˆ Q…^ n† ˆ Ti 1 …Wk …^ n††;

…3:7†

and P is the mapping from the master subcell to the local coordinates of the feature  l; ^ k 7! X P :w

n ˆ P …^ n† ˆ Rm 1 …Wk …^ n††;

…3:8†

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4097

^ g^†, respectively, W is where the vectors x, x, n, and ^ n, denote the coordinate pairs …x; y†, …x; y†, …n; g†, and …n; the mapping from the master subcell to physical domain x ˆ Wk …^ n†:

^ k 7! wk ; Wk : w

…3:9†

T is the mapping from the master element to physical domain Ti : s^i 7! si ;

x ˆ Ti …n†

…3:10†

and R is the mapping from the local feature coordinates to the physical domain l ! 7 X; Rm : X

x ˆ Rm … x†:

…3:11† R

Let us now also address the computation of the boundary integrals C g/i , in the case that the boundary C crosses an element s in multiple segments sj as is, for example, shown in Fig. 17. The contribution of the element s to fi , an entry in the global load vector, is fi js ˆ

nsegments X jˆ1

Z sj

g…s†/i …s† ds ˆ

nsegments X jˆ1

Z s^j

g…t†/i …t†

ds dt; dt

…3:12†

where t is a parameterization ‰ 1; 1Š on a master segment ^s as shown in Fig. 17(c). The numerical evaluation of this integral is done piecewise by employing a standard one-dimensional quadrature using a suitable order Gauss rule, or if the integrand happens to be highly oscillatory or singular at one or more points, an adaptive quadrature.

Fig. 17. Example of Neumann boundary intersected by an element of a GFEM III mesh.

4098

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 18. The mixed example problem for the wrench shaped domain, Problem O.

Fig. 19. The contours in the relative modulus of the gradient in the overkill solutions for Problem O. (a) The overkill solution for the FEM mesh; (b) the overkill solution for the GFEM II mesh.

3.2.1. The e€ect of the numerical integration tolerance Let us now demonstrate the e€ect of the tolerance employed to control the error in the numerical integrations on the accuracy of the computed solution. We employed the GFEM for solving the Laplacian in the wrench shaped domain perturbed by circular holes and straight cracks as shown in Fig. 18. Figs. 19(a) and (b) show the shades of the computed gradient of the overkill solutions computed by employing a bi-p FE basis of degree p ˆ 6 using void functions of degree pvoid ˆ 3, and re-entrant corner functions of degree pcorner ˆ 1 at nlayers ˆ 0. Fig. 20 depicts the modulus of the gradient using a bi-p FE basis of degree p ˆ 2, and the vertex void and straight crack functions of degree pvoid ˆ 1, and pcrack ˆ 1, at nlayers ˆ 0 on the coarse GFEM II mesh, which are computed using integration tolerances erel ˆ 0:1, erel ˆ 0:05, erel ˆ 0:01, and erel ˆ 0:001 for the requested relative error in the stiffness matrix entries. Note, an integration tolerance erel > 0:05 results in integration errors which pollute the pointwise accuracy of the computed solution. On the other hand, an integration tolerance erel < 0:01 is suf®cient for this problem. The theory of the control of the error in the numerical integration was discussed in [1].

4. The GFEM with analytically known handbook functions Below, we will give illustrative examples which underline the advantages of the GFEM over the traditional FEM, namely its geometrical ¯exibility and superior accuracy.

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4099

Fig. 20. The contours of the relative modulus of the gradient of the computed GFEM II solutions using di€erent relative integration tolerances for the Fast Remeshing algorithm. (a) erel ˆ 0:1; (b) erel ˆ 0:05; (c) erel ˆ 0:01; (d) erel ˆ 0:001.

4.1. Description of the simple example problems First we will consider problems which can also be solved rather easily by the FEM. We will consider the Neumann or the mixed boundary value problem for the Laplacian in either the rectilinear polygonal domain shown in Fig. 21(a), which we will call Domain I, or the curvilinear polygonal domain shown in Fig. 21(b), which we will call Domain II. We will call these domains simple because they can be easily

Fig. 21. The simple domains employed in the illustration of the GFEM and their subdivision into a coarse FE mesh of quadrilaterals. (a) The simple rectilinear domain, which we will call Domain I. (b) The simple curvilinear domain, which we will call Domain II.

4100

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

subdivided into a traditional FE mesh, which is a decomposition of the domain into a set of non-overlapping quadrilaterals, as is, for example, shown in Fig. 21. We will employ the following boundary value problems for the Laplacian in Domains I and II as our simple example problems: Pure Neumann problem: Du ˆ 0 on X; ou ˆ $…2x on ou ˆ0 on

…4:1a†

y†
n

on C1 ;

…4:1b†

on C2 ; C3 ; . . . ; C10 :

…4:1c†

Mixed problem: Du ˆ 0 on X; ou ˆ $…2x on uˆ0

…4:2a†

y†
n

on C1 ;

on C2 ;

ou ˆ0 on

…4:2b† …4:2c†

on C3 ; C4 ; . . . ; C10 :

…4:2d†

Note that in both problems the only loads are due to the Neumann boundary condition ou ˆ $…2x on

y†
n

on C1 :

…4:3†

The domains of the simple example problems together with their boundary conditions are shown in Fig. 22. Below, we will report the values and the convergence of the global energy norm of the error in GFEM approximate solutions, namely s Z def

jjeGFEM jjU ˆ

X

j$eGFEM j2 ;

…4:4†

where def

eGFEM ˆ uEX

uGFEM ;

…4:5†

where uEX denotes the exact solution, and uGFEM the GFEM approximate solution of the boundary value problem. In the results below, we will compute the error in the GFEM solution uGFEM by employing instead of the unknown exact solution uEX , a suciently accurate approximate overkill solution uOVK , which will be chosen such that q q 2 2 2 2 jjuEX jjU jjuGFEM jjU jjuOVK jjU jjuGFEM jjU jjeGFEM jjU ˆ  ; …4:6† jjuEX jjU jjuOVK jjU jjuEX jjU where  means that the equality has been achieved for the ®rst few signi®cant digits. Let us now describe the main features of the overkill solutions for the example problems. Figs. 23 and 24 show a contour plot of the relative modulus of the gradient, and the value of the energy norm of the overkill solution, for Problems A, B, and C, D, respectively. These solutions were computed by employing the meshes which are shown in Figs. 23(c) and 24(c), which were generated by employing three uniform nested re®nements of the initial coarse meshes shown in Figs. 21(a) and (b), respectively, followed by six,

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4101

Fig. 22. The simple domains and the boundary conditions of the example problems. (a) Problem A: The pure Neumann problem on Domain I. (b) Problem B: The mixed problem on Domain I. (c) Problem C: The pure Neumann problem on Domain II. (d) Problem D: The mixed problem on Domain II.

respectively seven, nested re®nements of the elements which are connected to the vertices at the re-entrant corners, and using the bi-p FE basis of degree p ˆ 6, and additional vertex corner functions added at the vertices of the elements connected to the re-entrant corners. The addition of the vertex corner functions to the approximation is discussed below in Section 4.2, and in later sections. In order to determine the number of correct signi®cant digits in jjuOVK jjU , we also computed the GFEM solutions uGFEM for each problem on

4102

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 23. The overkill solutions on Domain I. The contours of the relative modulus of the gradient of the overkill solution for: (a) Problem A; (b) Problem B, and (c) The employed mesh. The overkill solutions were obtained using the bi-p FE basis of degree p ˆ 6 on the mesh with the vertex corner functions added at the vertices of the elements connected to the re-entrant corner vertices.

the same meshes used for the overkill, using the bi-p FE basis of degree p ˆ 1; . . . ; 5, and with the vertex corner functions added at the vertices of the element connected to the re-entrant corner vertices. Fig. 25 shows the p-convergence of the energy norm of the error in the GFEM solutions where we used the p ˆ 6 GFEM overkill solution and Eq. (4.6) for computing the error. Note that we get better than N …2=3†

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4103

Fig. 24. The overkill solutions on Domain II. The contours of the relative modulus of the gradient of the overkill solution for: (a) Problem C; (b) Problem D, and (c) The employed mesh. The overkill solutions were obtained using the bi-p FE basis of degree p ˆ 6 on the mesh with the vertex corner functions added at the vertices of the elements connected to the re-entrant corner vertices.

4104

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 25. p-Convergence of the energy norm of the error for the GFEM solutions constructed using the bi-p FE basis of degree p ˆ 1; . . . ; 5 with the vertex corner functions added at the vertices of the elements which are connected to the re-entrant corner vertices, on the meshes shown in Fig. 23(c) for Problems A, and B, and Fig. 24(c) for Problems C, and D, respectively. Table 2 p-Convergence of jjuGFEM jjU , the energy norm of the GFEM approximation uGFEM computed by employing the bi-p FE basis on the meshes shown in Fig. 23(c) for Problems A, and B, and Fig. 24(c) for Problems C, and D, and, in addition, one layer of vertex corner functions around each re-entrant corner (pcorner ˆ 1, nlayers ˆ 1)a p

1 2 3 4 5 6 %E a

Problem A

Problem B

Problem C

Problem D

N

jjuGFEM jjU

N

jjuGFEM jjU

N

jjuGFEM jjU

N

jjuGFEM jjU

4786 18514 41502 73750 115258 166026

20.572938 20.610058 20.611950 20.612193 20.612297 20.612330

4786 18514 41502 73750 115258 166026

15.613488 15.637612 15.638884 15.639053 15.639121 15.639142

4836 18714 41960 74574 116556 167906

20.800384 20.829377 20.830810 20.830997 20.831079 20.831103

4836 18714 41960 74574 116556 167906

16.093418 16.113976 16.114968 16.115100 16.115157 16.115173

<0.179%

<0.164%

<0.152%

<0.141%

%E is the percentage relative di€erence in the energy norm of the last two solutions in the sequence for each problem.

algebraic convergence. Table 2 reports the values of the energy norm of the computed solution jjuGFEM jjU , versus the degree p, and the number of degrees of freedom for the GFEM approximations of the solutions of Problems A±D. From these results it can be seen that the energy norm of the employed overkill solution has converged to at least ®ve signi®cant digits, and its relative value of the energy norm of the error is of the order 0.1%. The main characteristic of the solutions of Problems A±D, is their singular behavior at the re-entrant corners. More precisely, in the neighborhood of each corner, the exact solution can be written in the form uEX …r; h† ˆ

M X

Aj rkj fj …h† ‡ u~…r; h†;

…4:7†

jˆ1

where r; h is the polar coordinate system associated with the corner shown in Fig. 26(b), kj are the exponents of the algebraic singularity at the corner which depends on the angle and the type of the boundary conditions on the corner edges Ck and Cl , which emanate from the corner vertex, namely

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4105

Fig. 26. Re-entrant corners of Domain I. (a) The enumeration of the vertices at the re-entrant corners. (b) The polar coordinate system and geometrical parameters of a re-entrant corner.

kj ˆ

8 jp > >
if Ck and Cl are both Dirichlet or both Neumann;

> …2j > :

and

( fj …h† ˆ

1†p 2u

…4:8† if Ck and Cl have mixed boundary conditions

cos…kj h†

if Ck and Cl have Neumann boundary conditions;

sin…kj h†

if Ck and=or Cl have Dirichlet boundary conditions;

…4:9†

and u~ is a function which is smoother than the ®rst term and practically negligible in the neighborhood of the corner point with respect to the ®rst term. The coecients Aj are called the intensity factors. For the derivation of the asymptotic expansion (4.7)±(4.9) see [23, pp. 169±173], and the references therein. Table 3 lists the values of the interior angles, the exponents k1 , and the intensity factors A1 for all the re-entrant Table 3 The geometric parameters, exponents k1 , and intensity factors A1 , at the re-entrant corners of Domain I for Problem A Vertex

u (°)

c (°)

k1

A1

Vertex

u (°)

c (°)

k1

A1

4 12 15 16 17 18 19 20 21 22 23 24

270.0 242.1 202.5 225.0 225.0 212.5 202.5 205.0 205.0 212.5 222.5 225.0

0.0 203.1 270.0 112.5 157.5 202.5 235.0 257.5 282.5 307.5 340.0 22.5

0.6667 0.7435 0.8889 0.8000 0.8000 0.8471 0.8889 0.8780 0.8780 0.8471 0.8090 0.8000

0.09583 )0.2916 1.265 )5.268 )4.028 )1.572 0.3163 1.522 4.207 5.178 4.807 2.541

25 26 27 28 29 30 31 32 35 36 .. . 62

225.0 225.0 214.9 232.5 240.0 234.9 240.0 232.6 270.0 270.0 .. . 270.0

67.5 145.0 190.0 225.0 277.5 337.5 32.4 92.4 167.5 77.5 .. . 185.0

0.8000 0.8000 0.8373 0.7742 0.7500 0.7662 0.7500 0.7738 0.6667 0.6667 .. . 0.6667

)3.830 )4.979 )2.634 0.2499 3.602 3.709 0.2187 )4.182 )0.5493 )1.425 .. . )0.9698

4106

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

corners for Domain I according to the numbering shown in Fig. 26(a). Here we assumed that the same type of boundary conditions are applied on the edges emanating for each re-entrant corner, as is the case for the boundary conditions employed in Problems A and B. 4.2. The e€ectiveness of the GFEM on simple domains We will now illustrate the e€ectiveness of the three versions of the GFEM in solving the simple example problems described in Section 4.1.

Fig. 27. The sequence of meshes which are used for the uniform h-convergence for the various versions of GFEM on Domain I. (a) GFEM I meshes; (b) GFEM II meshes; (c) GFEM III meshes.

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4107

Example 1 (Uniform h-convergence of the various versions of GFEM). We considered Problem A (the pure Neumann on the rectilinear polygonal domain), and computed the GFEM I, GFEM II, and GFEM III approximate solutions by employing the meshes shown in Figs. 27(a), (b) and (c), respectively, and the bilinear (p ˆ 1) ®nite element basis. For this choice of basis the GFEM I coincides with the classical FEM. Fig. 28(a) shows the graph of the relative error versus the number of degrees of freedom for the three versions of GFEM. From the theoretical analysis we expect that asymptotically jjeGFEM jjU  CN jjuEX jjU

…1=2† min…p;k†

ˆ CN

…1=3†

;

…4:10†

Fig. 28. The uniform h-convergence of the three versions of GFEM for Problem A, with bilinear (p ˆ 1) basis with, and without, corner functions. (a) The log±log graph of the relative error versus the number of degrees of freedom. (b) The values of C, versus the logarithm of the number of degrees of freedom.

4108

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

where p is the degree of the FE basis, N is the number of degrees of freedom and k is the smallest exponent of the algebraic singularity at the re-entrant corners. In this example p ˆ 1, and k ˆ 2=3, and hence we expect the energy norm of the error to converge to zero as N …1=3† . This is con®rmed by computing the values Cˆ

jjeGFEM jjU 1=3 N ; jjuEX jjU

…4:11†

which are shown in Fig. 28(b). Let us also report the e€ect of the vertex corner functions on the convergence of the GFEM. Fig. 28 compares the convergence of the relative value of the energy norm of the error of the GFEM III approximation with and without the vertex corner functions. Here, we employed the GFEM III meshes shown in Fig. 29. These meshes were obtained from the corresponding meshes of Fig. 27(c) by re®ning the elements in the neighborhood of the corners in order to satisfy the rule that any bilinear hat function may not have support over more than one re-entrant corner. Fig. 30 depicts the shades of the relative modulus of the computed gradient  j$uFEM j

1 jXj



Z X

j$uGFEM j ;

obtained using GFEM II and III approximations, without corner functions, and GFEM III approximation with pcorner ˆ 1 of corner functions at nlayers ˆ 0 and 1, with about the same number of degrees of freedom. From this example it can be concluded that: · The theoretical asymptotic rate of convergence of the GFEM for uniform meshes is practically achieved for coarse meshes. · The addition of vertex corner functions and the increase of nlayers shift the convergence curve downward. · Increasing the number of layers from nlayers ˆ 0 to 1 reduces the error further, and shifts the convergence curve downward.

Example 2 (The effect of the Dirichlet conditions on the accuracy of the GFEM). The e€ect of the Dirichlet boundary conditions on the accuracy of the various versions of the GFEM, can be seen from the uniform h-convergence results for Problem B (the mixed boundary value problem on the rectilinear polygonal domain), in which a zero Dirichlet condition is imposed on the polygonal boundary C2 . Here, we will denote by GFEM II(a), the GFEM II approximation which employs the sequence of meshes shown in Fig. 27(b), and by GFEM II(b), the GFEM II approximation which employs the sequence of meshes shown in Fig. 31. For the GFEM I, and the GFEM II(a) and GFEM II(b) meshes, the Dirichlet boundary conditions can be imposed exactly as in the standard FEM, namely, by zeroing the degrees of freedom which are located on C2 , the Dirichlet part of the boundary. In the case that the meshes overlap the Dirichlet boundary, e.g., the meshes shown in Fig. 27(c), the zero Dirichlet condition is imposed by zeroing the degrees of freedom, the support of which intersects the Dirichlet boundary as is, for example, shown in Fig. 32(a). The effect of the approach used to impose the Dirichlet conditions on the accuracy can be seen in Fig. 32(b), which shows the convergence of the energy norm of the error for the various versions of GFEM for Problem B. Note that the GFEM III convergence curve is shifted upward when compared with the GFEM I, GFEM II(a) and GFEM II(b) convergence curves. The accuracy of the GFEM III may be improved by including in the approximation vertex special functions, obtained from the special functions which satisfy the Dirichlet condition, i.e., the corner functions and corresponding functions which satisfy the zero Dirichlet boundary condition on the straight edges.

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4109

Fig. 29. The sequence of meshes employed in the uniform h-convergence of GFEM III with corner functions for Domain I. These meshes are obtained from the corresponding meshes of Fig. 27(c) by enforcing the rule that any vertex corner function may overlap only one re-entrant corner.

From this example it can be concluded that: · The geometrical ¯exibility and accuracy of GFEM may be retained even in the case of mixed boundary conditions by employing meshes which conform to the Dirichlet boundaries. · In the case that the meshes overlap the Dirichlet boundary there is a loss of accuracy which may be improved by including special functions which satisfy the Dirichlet condition.

4110

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 30. Relative modulus of the gradient and energy norm of the computed solution for four GFEM approximations of the solution of Problem A with about the same number of degrees of freedom. (a) GFEM II, p ˆ 1, no corner functions; (b) GFEM III, p ˆ 1, no corner functions; (c) GFEM III, p ˆ 1, pcorner ˆ 1 at nlayers ˆ 0; (d) GFEM III, p ˆ 1, pcorner ˆ 1 at nlayers ˆ 1.

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4111

Fig. 31. The sequence of GFEM II(b) meshes.

Fig. 32. The e€ect of the imposition of the zero Dirichlet condition on the GFEM III. (a) The Dirichlet condition is imposed by zeroing the degrees of freedom for the FE basis functions with support overlapping the boundary, and by adding special vertex functions which satisfy the Dirichlet condition. (b) The uniform h-convergence of the various versions of GFEM for Problem B (mixed problem), with bilinear (p ˆ 1) basis, with, and without, corner functions. The inferior accuracy of GFEM III is due to the approach used for imposing the Dirichlet condition.

4112

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 33. The sequence of meshes which are used for the uniform h-convergence for the various versions of GFEM on Domain II. (a) GFEM I meshes; (b) GFEM II(b) meshes; (c) GFEM II(c) meshes, and (d) GFEM III meshes.

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4113

Fig. 34. The uniform h-convergence of the various versions of GFEM on Domain II. (a) Problem C; (b) Problem D. Note the inferior accuracy of GFEM III for Problem D (mixed problem), which is due to the approach used for imposing the Dirichlet condition.

Example 3 (Uniform h-convergence of the various versions of the GFEM on curvilinear domains). We considered Problem C (the Neumann problem on the curvilinear polygonal domain), and Problem D (the mixed problem on the curvilinear polygonal domain), and computed the approximate solutions using the sequences of the meshes employed for GFEM I±III, as shown in Figs. 33(a), (b) and (c) respectively. Figs. 34(a) and (b) show the h-convergence graphs for Problems C and D, respectively. Note that the results are very similar with the corresponding results obtained in the case of the rectilinear domains. Figs. 35 and 36 show examples of the contours of the modulus of the computed gradient for the GFEM II(b) and (c), respectively, approximations of Problem D using bilinear …p ˆ 1† basis.

4114

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 35. Relative modulus of the gradient computed by a GFEM II approximation of Problem D. For (a) The coarse mesh; and (b) the mesh obtained from one uniform re®nement of the coarse mesh.

Fig. 36. Relative modulus of the gradient computed by a GFEM II approximation of Problem D. For (a) The coarse mesh; and (b) the mesh obtained from one uniform re®nement of the coarse mesh.

Example 4 (Adaptive h-convergence of GFEM III). Let us also show that the adaptive error control capabilities which were developed for the FEM, can be also employed for the GFEM. We computed adaptive GFEM III approximations of the solution of Problem A, by employing a bi-p FE basis, with p ˆ 1 and 2,

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4115

Fig. 37. An illustration of how vertex special functions are assigned to new vertices which are created by the re®nement of the mesh. A vertex special function is assigned to a new vertex of a child element only if it is already assigned to any of the pre-existing vertices of the child element.

with and without vertex corner functions. As initial mesh we employed the coarsest mesh of Fig. 27(c), and, in the case that the corner functions are added, the mesh shown in Fig 29(a). In each case we constructed a sequence of adaptively re®ned meshes by subdividing the elements for which the energy norm of the error is greater than c times the maximum element energy norm of the error, namely c max jjeGFEM jjU…s† 6 jjeGFEM jjU…s† : si

…4:12†

Here we used c ˆ 3=4, and the estimates of the element norms of the error jjeGFEM jjU…s† , obtained by employing an overkill on the same mesh, with the same corner functions and p ‡ 2 ®nite element basis. Each time that an element with corner functions is re®ned, the resulting new vertices are assigned corner functions, only if they happen to belong to an element which already has a vertex with corner functions assigned to it, as is, for example, shown in Fig. 37. The relative error was determined, as before, by employing the approximation q 2 2 jjuOVK jjU jjuGFEM jjU jjeGFEM jjU  jjuOVK jjU jjuEX jjU

…4:13†

using the overkill solution uOVK discussed in the previous section and shown in Fig. 23(a). Figs. 38(a) and 39(a) respectively, show the adaptive h-convergence graphs for p ˆ 1, and 2. In this case, the theoretical asymptotic rate is jjeGFEM jjU  CN jjuEX jjU

…1=2†p

…4:14†

This rate is achieved only for very re®ned meshes, while for coarser meshes a better rate may be obtained, as it can be seen from Figs. 38(b) and 39(b), which graphs the values of Cˆ

jjeGFEM jjU …1=2†p N : jjuEX jjU

…4:15†

This is a typical behavior of h-adaptive convergence, which has also been observed in the h-adaptive convergence of the standard FEM in polygonal domains, especially for elements of degree p greater than 1. Figs. 40 and 41, respectively, show examples of GFEM I and III meshes from the corresponding h-adaptive mesh sequences, and the corresponding GFEM solutions. Fig. 42 compares three h-adaptive meshes for

4116

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 38. The adaptive h-convergence of the three versions of GFEM for Problem A, for p ˆ 1, with, and without, corner functions. (a) The log±log graph of the relative error versus the number of degrees of freedom; (b) The values of C, versus the logarithm of the number of degrees of freedom.

p ˆ 2, from the GFEM I and GFEM III, adaptive sequences. It is interesting to note that the GFEM I meshes are re®ned also away from the corners, while the GFEM III meshes are re®ned primarily at the corners, and are more economical than the GFEM I meshes for suf®ciently high accuracy. This is due to the unsmoothness of the mappings employed in the GFEM I, and is a result of the choice of the employed initial meshes. We can make the following conclusions: · The h-adaptive versions of the GFEM give nearly optimal results, similar to the h-adaptive FEM. · The use of the corner functions leads to preasymptotic rates which are better than the theoretical asymptotic rate.

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4117

Fig. 39. The adaptive h-convergence of the three versions of GFEM for Problem A, for p ˆ 2, with, and without, corner functions. (a) The log±log graph of the relative error versus the number of degrees of freedom. (b) The values of C, versus the logarithm of the number of degrees of freedom.

Example 5 (A posteriori estimation of the energy norm of the error in the GFEM). Let us also give some examples of a posteriori estimation for the GFEM based on a recovered GFEM solution uREC GFEM . Let xX be the vertex patch, which is de®ned as the support of the piecewise bilinear hat function uX corresponding to X the vertex X, as shown in Fig. 43. In each xX , we construct a recovered GFEM solution uREC;x GFEM , in the form X def uREC;x GFEM ˆ

NX X jˆ0

0

aXj wxj X

‡

NX X NX ‡1

aXj wxj X ;

…4:16†

4118

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 40. Final approximate solutions and meshes obtained using the h-adaptive GFEM I, and GFEM III with p ˆ 1, pcorner ˆ 1 at nlayers ˆ 1, for Problem A. (a), and (b): The ®nal GFEM I mesh, and the corresponding solution. (c), and (d): The ®nal GFEM III mesh, and the corresponding solution. Note that the GFEM I mesh is much more re®ned away from the corners than the GFEM III mesh; this is due to the unsmoothness of the mappings for the employed GFEM I mesh. N ‡N 0 ‡1

X X where fwxj X gjˆ1 is the set of basis functions employed in the recovery. For example, for vertices in the interior of the mesh or at a Neumann boundary we employ the NX ˆ 2p0 ‡ 1 harmonic monomials centered at the vertex X, namely

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4119

Fig. 41. Intermediate approximate solutions and meshes obtained using the h-adaptive GFEM I, and GFEM III, with p ˆ 2, and pcorner ˆ 1 at nlayers ˆ 1, for Problem A. (a), and (b) The ®nal GFEM I mesh, and the corresponding solution, (c), and (d) the ®nal GFEM III mesh and the corresponding solution. As in the case p ˆ 1, the GFEM I meshes are more re®ned away from the corners than the corresponding GFEM III meshes. …i†

…i†

…i†

zX †j †; w2j …x; y† ˆ I……z zX †j † …4:17† p for j ˆ 1; 2; 3; . . ., where z ˆ x ‡ Iy, where I ˆ 1, and for vertices in the nlayers neighborhood of the ith corner, we also add the re-entrant corner functions w0 …x; y† ˆ 1;

w2j 1 …x; y† ˆ R……z

4120

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 42. Corresponding meshes from the h-adaptive GFEM I, and GFEM III mesh sequences for p ˆ 2. The re®nement of the GFEM I meshes away from the corners is due to the unsmoothness of the mappings, for the employed GFEM I meshes.

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4121

Fig. 43. Examples of the vertex patches xX in the construction of the recovered solution uREC GFEM .

…i†

…i†

X wxj‡N …r; h† ˆ rkj fj …h†;

…4:18†

for j ˆ 1; 2; 3; . . .. The coecients aXj ; j ˆ 0; . . . ; NX0 , are determined for each vertex X by solving the leastsquares minimization problem 2

X jjuREC;x GFEM

uGFEM jj;xi ˆ min;

…4:19†

where0 jj
jj;xX is a discrete L2 norm in xX . After all the sets of vertex coecients NX faXj gjˆ0 have been recovered, uREC GFEM , the global recovered GFEM solution, is constructed in the form uREC GFEM ˆ

X X 2N…D†

0 uX @

NX ‡NX0

X jˆ1

1 aXj wxj X A;

…4:20†

where N…D† denotes the set of vertices of the mesh D. Employing uREC GFEM in the de®nition of the error instead of uEX , we obtain EREC ˆ

r  X 2 jjuREC u jj GFEM GFEM U…s† ;

…4:21†

s2D

which is an estimator for jjeGFEM jjU , the global energy norm of the error. Employing EREC , we can also compute EREC REC def EREL ˆ q ; jjuGFEM jj2U ‡ …EREC †2

…4:22†

4122

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

which is an estimator for the relative error def

EREL ˆ

jjeGFEM jjU : jjuEX jjU

…4:23†

The accuracy of EREC is measured by the e€ectivity index def

jREC ˆ

EREC jjeGFEM jjU

…4:24†

REC by the e€ectivity index and the accuracy of EREL def

jREC REL ˆ

REC EREL : EREL

…4:25†

Further details about the error estimator and the steps used for its veri®cation are given in Section 4.3. REC Let us now illustrate the accuracy of EREC and EREL . Let us ®rst consider the GFEM approximations for Problems A (the Neumann problem in Domain I), and C (the Neumann problem in Domain II) for the REC sequences of uniform meshes shown in Figs. 27 and 33, respectively. Table 4 shows the values of EREL and EREL versus the number of degrees of freedom for Problem A, in the case that a bilinear uGFEM , is computed using only a bi-p FE basis, with p ˆ 1, and 2. In these computations, the recovered solution uREC GFEM , was constructed by using a basis of p0 ˆ 2p ‡ 1 harmonic polynomials on each patch, i.e., NX ˆ 2p0 ‡ 1 ˆ 4p ‡ 3 harmonic monomials, and one re-entrant corner function at the nlayers ˆ 0 neighborhood of each re-entrant corner. Table 5 reports the same results for Problem B. We note that: REC · In the case of the uniform meshes, the estimators EREC , and EREL , underestimate jjeGFEM jjU and EREL , respectively. This is because the main contribution to the error and the estimate is due to the elements with a vertex at the re-entrant corners, and the element indicator def

gREC ˆ jjuREC s GFEM

uGFEM jjU…s†

…4:26†

underestimates the energy norm of the error jjeGFEM jjU…s† in these elements. Table 4 REC , the exact relative error EREL , and the e€ectivity indices jREC , and jREC The estimated relative error EREL REL , for GFEM I, II, and III, solutions of Problem A, computed using the uniform meshes shown in Fig. 27 and bi-p FE basisa pˆ1

pˆ2

N

REC EREL

EREL

jREC

jREC REL

N

REC EREL

EREL

jREC

jREC REL

GFEM I

212 752 2792

0.1191 0.0802 0.0523

0.1919 0.1217 0.0771

0.614 0.656 0.678

0.621 0.659 0.678

752 2792 10712

0.0801 0.0435 0.0257

0.0867 0.0518 0.0314

0.924 0.839 0.817

0.924 0.840 0.818

GFEM II

107 367 1337

0.1099 0.1077 0.0762

0.2949 0.2021 0.1343

0.358 0.525 0.564

0.373 0.533 0.567

367 1340 5077

0.1087 0.0688 0.0467

0.2013 0.1209 0.0717

0.532 0.566 0.651

0.540 0.570 0.651

GFEM III

158 498 1674

0.1183 0.0948 0.0623

0.2558 0.1698 0.1108

0.450 0.553 0.560

0.462 0.558 0.562

563 1843 6356

0.0915 0.0698 0.0382

0.1645 0.1030 0.0633

0.551 0.676 0.604

0.556 0.678 0.603

a REC The estimators EREC and EREL are based on a recovered solution uREC GFEM , constructed using a 2p ‡ 1 degree harmonic polynomial basis at each vertex, supplemented by one corner function at each re-entrant corner vertex.

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4123

Table 5 The same results as in Table 4, but for Problem C, and its GFEM I, II, and II, solutions computed using the meshes shown in Fig. 33 pˆ1

pˆ2

N

REC EREL

EREL

jREC

jREC REL

N

REC EREL

EREL

jREC

jREC REL

GFEM I

212 752 2792

0.1263 0.0769 0.0485

0.1731 0.1091 0.0693

0.724 0.702 0.699

0.730 0.704 0.700

752 2792 10712

0.0706 0.0399 0.0240

0.0769 0.0469 0.0287

0.917 0.851 0.836

0.917 0.852 0.836

GFEM II

97 321 1140

0.1421 0.1235 0.0854

0.2947 0.2125 0.1389

0.466 0.572 0.611

0.482 0.581 0.615

342 1189 4333

0.1730 0.0812 0.0554

0.2116 0.1282 0.0793

0.811 0.630 0.697

0.803 0.634 0.713

GFEM III

57 163 528

0.1513 0.1119 0.0891

0.3321 0.2391 0.1592

0.435 0.457 0.555

0.456 0.468 0.560

196 583 1961

0.1616 0.0837 0.0646

0.2267 0.1538 0.0966

0.703 0.540 0.667

0.698 0.544 0.669

In Example 4 we gave examples of h-adaptive GFEM computations which were based on the exact element indicators def

gEX s ˆ jjuEX

uGFEM jjU…s† ;

…4:27†

which were obtained via an overkill. Let us now show that the employment of the computed element error indicators gREC leads practically to the same error convergence curves as the exact indicators gEX s s , and that for suciently re®ned meshes, the global e€ectivity index jREC is close to 1. Fig. 44(a) compares the convergence of relative error versus the number of degrees of freedom, for the h-adaptive mesh sequences obtained using the exact indicators gEX and the computed indicators gREC for Problem A with bi-p (p ˆ 1 s s and 2) FE basis, and corner functions. Fig. 44(b) shows the graph of the effectivity index jREC , versus the logarithm of the number of degrees of freedom. We note that: · The h-adaptive mesh sequences obtained, using the exact and the computed indicators, result in practically the same error convergence curves. · As the mesh tends toward being equilibrated with respect to the element error indicators the e€ectivity index jREC , of the estimator EREC , tends toward one. Example 6 (p-Convergence of the GFEM). Let us also give some examples of p-convergence of the GFEM. We considered Problem A (the pure Neumann on the rectilinear polygonal domain), and computed the GFEM I approximate solutions by employing the ®rst two uniform meshes as shown in Fig. 27, with the bip (p ˆ 1; 2; 3; . . . ; 6) ®nite element basis, with the vertex corner functions pcorner ˆ 1 or 2 at nlayers ˆ 0, or 1 added to the approximation. Figs. 45(a) and (b) show the graph of the relative error versus the number degrees of freedom for various choices of the basis. From the theoretical analysis we expect that without the corner functions, asymptotically we have jjeGFEM jjU  CN jjuEX jjU

k

ˆ CN

…2=3†

;

…4:28†

which is achieved, as can be seen in Fig. 45 (see Figs. 46 and 47). We note that: · The addition of the vertex corner functions uX rk1 f1 …h†, corresponding to the ®rst term of the asymptotic expansion at the re-entrant corner vertices shifts the convergence curve of the p-extensions downward. The addition of the vertex corner functions uX rk f …h† corresponding to the second- or higher-order terms of the asymptotic expansions at the re-entrant corners has negligible effect on the global accuracy. · The addition of the vertex corner functions uX rk1 f1 …h†, at nlayers ˆ 1 leads to signi®cant gain in the accuracy in the pre-asymptotic range.

4124

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 44. The adaptive h-convergence of GFEM I and GFEM III on Problem A using the exact and the computed indicators. (a) The log±log graph of the relative error versus the number of degrees of freedom, results for the computed indicators are shown by gray symbols. (b) The effectivity index versus the logarithm of the number of degrees of freedom.

4.3. Veri®cation of the implementation Let us also discuss how we veri®ed the implementation of the GFEM approximation, and the leastsquares recovery error estimator. One way to check the implementation of the GFEM approximation is to employ example problems with exact solutions from the span of the GFEM approximation, for which we expect that uEX  uGFEM : Let us give an example.

…4:29†

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4125

Example 7 (Verification of the implementation of GFEM III). Consider a patch test on the GFEM III mesh shown in Fig. 48, and the GFEM space constructed using the bilinear FE functions, and the elliptical harmonic functions for the voids of degree pvoid ˆ 1 at nlayers ˆ 0. We employed uEX ˆ 2x

y

…4:30†

and applied Neumann boundary conditions at all boundaries, and using an integration tolerance e ˆ 0:0001, we obtained jjuGFEM jjU ˆ 13:4233819485036629;

…4:31†

whereas the exact value is jjuEX jjU ˆ 13:423305730223793:

…4:32†

Note that we have six-digit accuracy in the energy norm, and the loss of accuracy is due to the tolerance employed for the numerical integration. Similarly, to check the implementation of the least-squares recovery, we can employ cases for which 

uREC GFEM  uEX ;

…4:33†



where uREC GFEM is the recovered solution obtained by sampling the values of uEX , and we have uEX



uGFEM ˆ uREC GFEM

uGFEM

…4:34†

and the element e€ectivity indices must be equal to unity in all the elements. For the same exact solution uEX , we should expect that jjuREC uGFEM jjU GFEM  1; jjuEX uGFEM jjU

…4:35†

where uREC GFEM is the recovered solution from the sampled values of uGFEM . Let us now give an example. Example 8 (Verification of the implementation of the recovery estimator). We employed the Neumann boundary value problem for the Laplacian in Domain II, with Neumann boundary conditions compatible with the exact solution uEX ˆ

3 X

…ai R…zi † ‡ bi I…zi ††;

…4:36†

iˆ1

where z ˆ x ‡ Iy, and we arbitrarily chose ai ˆ f4; 3; 1g, bi ˆ f 5; 3; 1g. We computed a GFEM I solution and a GFEM II solution employing the bi-p FE basis with p ˆ 2, while in the least-squares recovery we used the space of harmonic monomials with degree p0 ˆ 3. Figs. 49(a) and (b) show the values of the element effectivity indices ˆ jREC s

gREC s jjuEX jjU …s†

…4:37†

for the computed indicators on the GFEM I and III mesh, respectively. Note that the values are close to 1. When sampling the uEX instead, we obtained e€ectivity indices jREC ˆ 1 to within seven digits in all eles ments. 4.4. Description of the example problems in complex domains Next, we will consider problems which cannot be easily solved by the standard FEM because the domain is complex, and a classical FE mesh is dif®cult, or practically impossible to construct. As examples of domains of increasing complexity we will employ the domains shown in Figs. 50(a)±(c) which may be

4126

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 45. The p-convergence of GFEM I and the effect of the corner functions. The bi-p FE basis was employed with and without vertex corner functions of degree pcorner ˆ 1; 2 at nlayers around the re-entrant corners. (a) The log±log graph of the relative value of the energy norm of the error versus the number of degrees of freedom for the uniform mesh shown in Fig. 47(a). (b) The log±log graph of the relative value of the energy norm of the error versus the number of degrees for the uniform mesh shown in Fig. 47(b).

understood as perturbations of the original Domains I and II, respectively, by the inclusion of a large number of elliptical voids and/or cracks. We will call the domain shown in Fig. 50(a) Domain I0 , the domain shown in Fig. 50(b) Domain I00 , and the domain shown in Fig. 50(c) Domain II0 . We will be interested in solving the same boundary value problems which we solved in the previous section in Domain I0 , Domain I00 , and Domain II0 , with zero Neumann boundary conditions imposed at the boundaries of the elliptical voids and cracks, namely: Pure Neumann problem: Du ˆ 0 on X;

…4:38a†

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4127

Fig. 46. The p-convergence of GFEM I and the effect of the corner functions in a strongly graded mesh. The bi-p FE basis was employed with and without vertex corner functions of degree pcorner ˆ 1; 2 at nlayers ˆ 0; 1 around the re-entrant corners. The mesh has 5 nested re®nements at the re-entrant corners as shown in Fig. 47(c).

ou ˆ $…2x on

y†
n

on C1 ;

…4:38b†

ou ˆ0 on

on C2 ; C3 ; . . . ; C10 ;

…4:38c†

ou ˆ0 on

on C11 ; . . . ; Cnvoids :

…4:38d†

Mixed problem: Du ˆ 0 on X;

…4:39a†

Fig. 47. The uniform meshes used in the p-convergence of GFEM I on Domain I. (a) Mesh with one nested re®nement; and (b) two nested re®nements, and (c) with 5 re®nements at the re-entrant corners.

4128

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

ou ˆ $…2x on

y†
n

on C1 ;

…4:39b†

u ˆ 0 on C2 ; ou ˆ 0 on C3 ; C4 ; . . . ; C10 ; on ou ˆ 0 on C11 ; . . . ; Cnvoids : on

…4:39c† …4:39d† …4:39e†

Here nvoids is the number of internal voids and/or cracks in the employed domain. We will call Problems E and F, respectively, the Neumann, and the mixed problem in Domain I0 , Problems G, and H, the Neumann and mixed problem in Domain I00 , and Problems I and J, the Neumann and mixed problem in Domain II0 . Fig. 51 illustrates the domains and the boundary conditions for these six problems. A solution for Problems E; F; . . . ; J can practically only be determined by employing either a GFEM II, or a GFEM III approximation. For example, for the problems posed in the rectilinear domain we will employ the GFEM I mesh we used for the corresponding problems in the simple domain, which now becomes a GFEM II mesh for the perturbed domain. Figs. 52±54 show a contour plot of the relative modulus of the gradient, and the value of the energy norm of overkill solutions, for Problems E, F, Problems G, H, and Problems I, J, respectively. As in Section 4.1, the overkill solution uOVK was de®ned as the last member of a sequence of GFEM II approximations with bi-p ®nite element basis, from p ˆ 1 to 3 or 4, with vertex corner, void, or crack functions of degrees pcorner ˆ 1, pvoid ˆ 1 or 3, and pcrack ˆ 1, respectively, added to the basis of the approximation. The employed mesh for Problems E, F, and Problems G, H was generated from the coarse FE mesh in the corresponding unperturbed domain, by employing three uniform nested re®nements of the square domain, followed by six nested re®nements of the elements connected to the vertices at the re-entrant corners. Table 6 reports the values of the energy norm of the computed solution jjuGFEM jjU , versus the degree p and the number of degrees of freedom for the GFEM approximations of the four problems E, F, G and H. The mesh employed in the solution of Problems I and J was generated by employing four uniform nested re®nements, followed by six nested re®nements of elements connected to re-entrant corners. Table 7 reports the results for Problems I and J (see Figs. 53±58).

Fig. 48. The GFEM III mesh used in the patch test for uGFEM on Domain I0 .

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4129

Fig. 49. The element e€ectivity indices, jREC , for the GFEM I and GFEM III test for the least squares recovery which samples the s computed GFEM solution. (a) GFEM I; (b) GFEM III.

4130

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 50. The complex domains obtained by perturbing the original simple domains by including a large number of elliptical voids and/or cracks. (a) Domain I0 : The rectilinear domain with 35 elliptical voids. (b) Domain I00 : The rectilinear domain with 158 cracks. (c) Domain II0 : The curvilinear domain with 458 voids and cracks.

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4131

Fig. 51. The complex example problems. (a) Problem E: The pure Neumann problem in Domain I0 . (b) Problem F: The mixed problem in Domain I0 . (c) Problem G: The pure Neumann problem in Domain I00 . (d) Problem H: The mixed problem in Domain I00 . (e) Problem I: The pure Neumann problem in Domain II0 . (f) Problem J: The mixed problem in Domain II0 .

4132

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 52. The overkill solutions for the boundary value problems in Domain I0 . The contours of the relative modulus of the gradient of the overkill solution for: (a) Problem E, and (b) Problem F. (c) The employed mesh.

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4133

Table 6 p-Convergence of jjuGFEM jjU , the energy norm of the GFEM approximation uGFEM computed by employing the bi-p FE basis with degree p ˆ 1; . . . ; 4, and adding a selection of special functions for Problems E, F, G, and Ha p

1 2 3 4 5

Problem E; pcorner ˆ 1; pvoid ˆ 3

Problem G; pcorner ˆ 1; pcrack ˆ 1

Problem H; pcorner ˆ 1; pcrack ˆ 1

N

jjuGFEM jjU

N

jjuGFEM jjU

N

jjuGFEM jjU

N

jjuGFEM jjU

8366 22090 45066 77294 118774

22.205419 22.256259 22.258737 22.259041 22.259145

8366 22090 45066 77294 118774

16.917480 16.956198 16.958051 16.958262 16.958312

7642 21370 44358 76606

21.735537 21.872867 21.881694 21.893280

7642 21370 44358 76606

16.647644 16.783134 16.790549 16.800897

%E a

Problem F; pcorner ˆ 1; pvoid ˆ 3

<0.523%

<0.499%

<3.25%

<3.51%

%E is the percentage relative di€erence in the energy norm of the last two solutions in the sequence for each problem.

We will also consider another example problem which is related to the micro-thermal analysis of real life composites, and was ®rst considered in [24]. We will let Domain III0 and Domain III00 , denote a square domain minus a set of circular voids as shown in Figs. 59(a) and (b), respectively. We will be interested in the solution of the following Neumann problem in Domain III0 and Domain III00 : Pure Neumann problem: Du ˆ 0 on X; ou ˆ $…2x on ou ˆ0 on

y†
n

…4:40a† on C1 ;

…4:40b†

on C2 ; . . . ; Cnvoids ‡1 :

…4:40c†

Here C1 is the outer boundary of the square domain, and C2 ; . . . ; Cnvoids ‡1 are the boundaries of internal voids. We will call the above Neumann problem posed on Domain III0 and Domain III00 , Problems K and L, respectively, as shown in Fig. 60. Figs. 62 and 63 show a contour plot of the relative modulus of the gradient, and the value of the energy norm of the overkill solutions, for Problems K and L, respectively. In these results the overkill solution uOVK was constructed using the GFEM III approximation on a mesh of squares with bi-p ®nite element basis with p ˆ 4, with void functions pvoid ˆ 1 added at nlayers ˆ 0. The meshes, shown in Figs. 62(b) and 63(b), are identical and were generated by employing ®ve nested re®nements of the square domain. Fig. 61 gives an example of a complex classical ®nite element mesh which was used in [24] to obtain comparable accuracy for a similar problem (Figs. 62±64).

Table 7 p-Convergence of jjuGFEM jjU , the energy norm of the GFEM approximation uGFEM computed by employing the bi-p FE basis with degree p ˆ 1; 2; 3 and adding a selection of special functions for Problems I, and Ja p

1 2 3 %E a

Problem I; pcorner ˆ 1; pvoid ˆ 1; pcrack ˆ 1

Problem J; pcorner ˆ 1; pvoid ˆ 1; pcrack ˆ 1

N

jjuGFEM jjU

N

jjuGFEM jjU

25723 59999 116951

23.647289 23.750562 23.772250

25675 59951 116903

18.305110 18.408457 18.413077

<4.27%

%E is the error in the energy norm of the p ˆ 2 result with respect to same with p ˆ 3.

<2.22%

4134

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 53. The overkill solutions for the boundary value problems in Domain I00 . The contours of the relative modulus of the gradient of the overkill solution for: (a) Problem G, and (b) Problem H. (c) The employed mesh.

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4135

Fig. 54. The overkill solutions for the boundary value problems in Domain II00 . The contours of the relative modulus of the gradient of the overkill solution for: (a) Problem I; (b) closeup of Problem I; (c) Problem J; (d) close up of Problem J; (e) The employed mesh.

4136

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 55. The integration mesh employed in the evaluation of all the integrals in the overkill solution for Problem E or F on Domain I0 . (a) The integration mesh for the entire domain; (b) and (c) details of the integration mesh.

Table 8 reports the values of the energy norm of the computed solution jjuGFEM jjU , versus the degree p, and the number of degrees of freedom for the GFEM approximations of Problem K and L. From these results it can be seen that the energy norm of the employed overkill solution has converged to at least four signi®cant digits, and its relative value of the energy norm of the error is of the order 0.3%. 4.5. Illustration of the GFEM on complex domains We will now illustrate the e€ectiveness of the GFEM for solving boundary value problems in complex domains, like the domains given in Section 4.4.

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4137

Fig. 56. The integration mesh employed in the evaluation of all the integrals in the overkill solution for Problem G or H on Domain I00 . (a) The integration mesh for the entire domain; (b) and (c) details of the integration mesh.

Example 9 (The geometrical flexibility of GFEM II and GFEM III, and their utility in obtaining accurate fixed mesh solutions). We considered Problem E (the pure Neumann problem in Domain I0 ) and computed four solutions using three di€erent GFEM II meshes, and one GFEM III mesh, using the bi-p FE basis of degree p ˆ 2, with added vertex void functions pvoid ˆ 1, at nlayers ˆ 0. Fig. 65 shows the employed meshes along with the modulus of the computed GFEM gradient. Note that the mesh of Fig. 65(a) overlaps the elliptical voids, the mesh of Fig. 65(b) overlaps the elliptical and the square voids, the mesh of Fig. 65(c) overlaps all boundaries except C1 , whereas the mesh of Fig. 65(d) overlaps the entire domain. All four meshes have elements roughly of the same size at the corners, and similar relative values of the energy norm of the error: 3.48%, 6.43%, 2.88%, and 3.70%, respectively.

4138

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 57. The integration mesh employed in the evaluation of all the integrals in the overkill solution for Problem I or J on Domain II0 . (a) The integration mesh for the entire domain; (b) and (c) details of the integration mesh.

We next used a GFEM III mesh of nested square elements to solve Problems I and J, employing a bi-p FE basis of degree p ˆ 2, with the addition of vertex void functions of degree pvoid ˆ 1, vertex crack functions of degree pcrack ˆ 1, and vertex corner functions of degree pcorner ˆ 1 functions, all at nlayers ˆ 0. The integration mesh used in the evaluation of all integrals is shown in Figs. 66. Fig. 67(a), (b) and (c), (d) depict the shades of the relative modulus of the computed gradient, along with the number of degrees of freedom, and the energy norm of the computed solution for Problems I and J, respectively. Fig. 68 shows

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4139

Fig. 58. p-Convergence of the energy norm of the error for the GFEM solutions for Problem A on the unperturbed Domain I versus for Problem E on the perturbed Domain I. The Problem A approximation was constructed using the bi-p FE basis of degree p ˆ 1; . . . ; 5 with the vertex corner functions added at the vertices of the elements which are connected to the re-entrant corner vertices at nlayers ˆ 1 on the mesh shown in Fig. 23(c). While the Problem E approximation was constructed using the bi-p FE basis of degree p ˆ 1; . . . ; 4 with the vertex corner functions of degree pcorner ˆ 1 and the elliptical void function of degree pvoid ˆ 1 added at nlayers ˆ 0 on the mesh shown in Fig. 52(c).

the shades of the computed gradient using a GFEM II mesh obtained from the re®nement of the coarse mesh shown in Fig. 4(c). Note that the results for Problem I, which are obtained using GFEM II and GFEM III, are of comparable accuracy, however, for Problem J the GFEM III solution is less accurate than the GFEM II solution, because (as we discussed above) of the approach used for the application of the Dirichlet boundary condition in GFEM III. The relative value of the energy norm of the error for the GFEM III computed solutions of Problems I and J is 8.71% and 19.24%, respectively, while the relative value of the energy norm of the error for the GFEM II computed solutions of Problems I and J is 10.5% and 11.1%. From this example we note: · GFEM II and GFEM III can be used to compute the solution of problems in very complex domains with good engineering accuracy, e.g., 5% accuracy in the energy norm. · In a domain perturbed by multiple voids and or cracks, a mesh designed to ®t the unperturbed domain can still be used. If the perturbation changes, i.e., topology of the internal structure changes, the GFEM is ¯exible enough that remeshing may be unnecessary. · As with the previous results on the unperturbed domain, GFEM II or GFEM III can be used to obtain approximations with the same order accuracy. Example 10 (The effect of the vertex corner, void, and crack functions on the accuracy of the GFEM). We considered various ®xed mesh GFEM approximations on Domains I0 and II0 , we used elements of low bi-p degree, p ˆ 1 or 2, and added the vertex void and vertex crack functions of degree pvoid ˆ 1 or 2 and pcrack ˆ 1 or 2 at nlayers ˆ 0, and we examined the effect of these special functions on the error in areas where the gradient of the solution takes on high values. First, we examined Problem E (the Neumann problem on the Domain I0 which includes 35 elliptical voids), when a ®xed GFEM III mesh of squares is employed with a bi-p FE basis of degree p ˆ 2, with three combinations of void and corner functions. In order to analyze the results, we also computed a separate

4140

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 59. Two square domains perturbed by arrays of circular voids. (a) Domain III00 : a square domain encompassing 393 voids; (b) Domain III0 : a square domain encompassing 429 voids; (c) Enhanced micrograph of an actual ®ber reinforced composite and the frames used for obtaining Domain III00 and Domain III0 (from [24]).

overkill solution on the same mesh, and used it to obtain the pointwise modulus of the error in the gradient of the computed GFEM III solutions, as shown in Figs. 69(a)±(f). The overkill used the bi-p FE basis of p ˆ 5, with vertex void functions pvoid ˆ 3 and vertex corner function of pcorner ˆ 1. The effect of adding the void functions can better be seen in the area where the gradients of the solution are high. Here, the quadratic p ˆ 2 approximation without any pvoid functions has a relative error in the gradient of at least more than 50% in the area immediately surrounding the voids. The addition of the void functions, with pvoid ˆ 1 at nlayers ˆ 0, reduces the local error around the voids to a practically negligible amount. The error is further reduced by the addition of the re-entrant corner functions, with pcorner ˆ 1 at nlayers ˆ 0, as shown in Figs. 69(e) and (f). Second, we examine Problem K (the voids in the square Domain III0 , here we employed a uniform mesh of squares and we applied two approximations, the ®rst with the bi-p basis p ˆ 2 without any vertex void functions, and the second with the bi-p basis p ˆ 2 with pvoid ˆ 1 at nlayers ˆ 0. The overkill approximation used to compute the error in these two examples was computed on the same mesh and used the bi-p p ˆ 3 basis and pvoid ˆ 2 functions at nlayers ˆ 0. The contours of the modulus of the error in the gradient for the

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4141

Fig. 60. Two complex example problems on domains perturbed by circular voids. (a) Problem K: the Neumann problem on Domain III0 ; (b) Problem L: the Neumann problem on Domain III00 .

Fig. 61. An example of a classical FE mesh used to discretize a square subdomain of the actual composite; from [24]. The mesh uses  20000 elements. Note the complexity of the mesh and the distortion of the elements.

4142

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 62. The overkill solution for Problem K. (a) The contours of the relative modulus of the gradient and the energy norm of the overkill solution; (b) detail of the contours; and (c) the employed mesh. The overkill solution was determined using the bi-p basis of degree p ˆ 4, with pvoid ˆ 1 vertex void at nlayers ˆ 0.

two approximations are shown in Figs. 70(a) and (b) along with the number of degrees of freedom used in the approximation and the resulting energy norm of the solutions. Third, we examined the accuracy of GFEM II approximations of Problem H (the mixed problem on the Domain I00 which includes 158 cracks), when a ®xed GFEM II mesh that respects all domain boundaries except for the cracks with re®nements at the re-entrant corners is employed as shown in Fig. 71(c). We employed two approximations, the ®rst using a bilinear FE basis (p ˆ 1) with added vertex crack functions of degree pcrack ˆ 1 at nlayers ˆ 0; the second using a biquadratic FE basis

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4143

Fig. 63. The overkill solution for Problem L. (a) The contours of the relative modulus of the gradient and the energy norm of the overkill solution. (b) The employed mesh. The overkill solution was determined using the bi-p basis of degree p ˆ 4, with pvoid ˆ 1 vertex void at nlayers ˆ 0.

(p ˆ 2), with added vertex crack functions of degree pcrack ˆ 2 at nlayers ˆ 0. Both approximations also included the pcorner ˆ 1 functions at nlayers ˆ 0 around the re-entrant corners. The overkill approximation used to compute the error in these two examples was computed on the same mesh with an additional four nested re®nements of elements at the re-entrant corners and used the bi-p ‡ 2 FE basis with added vertex crack functions of degree pcrack ‡ 1, and also pcorner ˆ 1 functions at nlayers ˆ 0. The contours of the modulus of the error in the gradient for the two approximations are shown in Figs. 71(a) and (b).

4144

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 64. The integration mesh employed in the evaluation of all the integrals in the overkill solution for Problem K on Domain III0 . (a) The integration mesh for the entire domain; (b) and (c) details of the integration mesh.

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4145

Table 8 p-Convergence of jjuGFEM jjU , the energy norm of the GFEM approximation uGFEM for Problems K, and L, computed by employing the bi-p FE basis with degree p ˆ 1; . . . ; 4 and adding vertex void functions, pvoid ˆ 1, at nlayers ˆ 0a p

Problem K; pvoid ˆ 1

1 2 3 4

N

jjuGFEM jjU

N

jjuGFEM jjU

14583 26978 47523 76218

550.26542 550.68523 550.70847 550.71119

13435 25832 46383 75088

532.02068 532.57787 532.59266 532.59460

%E a

Problem L; pvoid ˆ 1

<0.270%

<0.314%

%E is the percent error in the energy norm of the p ˆ 3 solution with respect to the p ˆ 4 solution.

We conclude: · The addition of the vertex corner, void and crack functions dramatically improves the global and the pointwise accuracy of the GFEM. · As in the standard FEM, the error in an area of interest can be partitioned into local and pollution error. Employing the idea of goal adaptive meshes presented in [25], it is possible to construct an adaptive algorithm for adding the special functions such that we achieve high accuracy only in the areas of interest. Example 11 (The local accuracy of the recovery error estimator). We considered Problem K and the bi-p p ˆ 1, pvoid ˆ 1 approximation and used the recovery error estimator with the functions p0 ˆ 2 and pvoid ˆ 2. Fig. 72(a) shows the pointwise estimated relative modulus of the error of the computed gradient  j$eREC GFEM j

1 jXj



Z X

j$uEX j

and (b) shows the exact error in the computed gradient. In order to obtain the exact error of the computed gradient $eGFEM , we used an overkill solution with the bi-p basis of degree p ˆ 3, with vertex void functions of degree pvoid ˆ 3 at nlayers ˆ 0. We conclude that in the cases that the pollution is negligible, as is the case presented in this example, the recovery estimator may be interpreted as the error. Example 12 (The h and p-convergence of GFEM in complex domains). We considered Problem E (the pure Neumann problem in the rectilinear domain shown in Fig. 50), and computed GFEM III approximate solution by employing only the bilinear p ˆ 1 ®nite element basis on the approximation mesh, by adding special corner and ellipse functions, and by employing a uniform or an adaptive mesh. For approximations that do, or do not use, the re-entrant corner functions, the initial meshes are shown in Figs. 73(a) and (b), respectively. The convergence of the relative error for the uniform and the adaptive mesh sequences is shown in Fig. 74. Fig. 75 shows two adaptive meshes with almost the same relative error, in the case that no special functions are used, and in the case that both corner and void functions are employed. We conclude that the adaptive h, p, and also hp, GFEM can be used to achieve high accuracy for problems in complex domains.

4146

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 65. Four di€erent GFEM approximations of the same problem illustrating the ¯exibility of meshing for complex domains. All four approximations use p ˆ 2, in addition, the vertex void functions of degree pvoid ˆ 1, at nlayers ˆ 0. The shades of the computed gradient are shown along with the corresponding mesh in each case. (a) The employed GFEM II mesh overlaps the elliptical void boundaries; (b) The employed GFEM II mesh overlaps both the elliptical and the square voids; (c) The employed GFEM II mesh overlaps all domain boundaries except C1 ; (d) The employed GFEM III mesh overlaps all domain boundaries.

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4147

Fig. 66. The integration mesh employed in the evaluation of all the integrals in the GFEM III approximations for Problem I or J. (a) The integration mesh for the entire domain; (b) and (c) details of the integration mesh.

4148

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 67. GFEM III solutions of the boundary value problems on Domain II0 , which includes 458 voids and cracks. GFEM III solution using a mesh of squares for (a)±(b) Problem I; and (c)±(d) Problem J.

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4149

Fig. 68. GFEM II solutions of the boundary value problems on Domain II0 , which includes 458 voids and cracks. Both GFEM II meshes were obtained by using four nested uniform re®nements and the approximations used the bilinear FE basis (p ˆ 1), with added vertex void functions pvoid ˆ 1 and vertex crack functions pcrack ˆ 1 at nlayers ˆ 0. (a)±(b) Problem I; and (c)±(d) Problem J.

4150

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 69. The e€ect of the vertex corner, void, and crack functions on the accuracy of the GFEM. Bi-quadratic (p ˆ 2) GFEM III approximation (a) without any added special functions; (c)±(d) with added vertex void functions of degree pvoid ˆ 1 at nlayers ˆ 0; (e)±(f) with added vertex void, and corner functions of degrees pvoid ˆ 2, and pcorner ˆ 1 respectively, both at nlayers ˆ 0. Note that although there is relatively little di€erence in the global accuracy between the three cases, there is a dramatic di€erence in the local accuracy.

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4151

Fig. 70. Two approximations for Problem K and their resulting contours in the relative modulus of error in the computed gradient in a local area. The biquadratic (p ˆ 2) GFEM II approximation (a) without any vertex functions; (b) with added vertex crack functions of degrees pcrack ˆ 1 at nlayers ˆ 0. (c) area of interest inside Domain III0 .

4152

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 71. Two approximations on a GFEM II mesh for Problem H and their resulting contours in the relative modulus of error in the computed gradient in a local area. The bi-linear (p ˆ 1) GFEM III approximation (a) with added vertex crack functions of degree pcrack ˆ 1 at nlayers ˆ 0; (b) with added vertex crack functions of degrees pcrack ˆ 2 at nlayers ˆ 0; (c) the approximation mesh and the boundary window of the local area.

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4153

Fig. 72. Comparison of the estimated error versus the exact error for the GFEM II approximation of Problem K. (a) The contours of the estimated error in the computed gradient; (b) contours of the exact error in the computed gradient; (c) area of interest.

4154

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 73. The initial meshes used in the uniform and adaptive h-convergence on the complex domain. (a) Initial mesh employed in the case that no corner functions are employed; (b) Initial mesh employed when corner functions are used.

Fig. 74. The uniform and adaptive h-convergence of the GFEM III on the rectilinear domain perturbed by multiple elliptical voids.

5. The GFEM with numerically generated handbook functions In this section, we discuss the implementation and give illustrative examples of the GFEM with numerically generated handbook functions. 5.1. Incorporation of the numerically generated handbook functions Let us outline the major steps in using the numerically generated handbook solutions in the GFEM approximation. First, let us assume a suitable set, or library, of handbook solutions has been pre-computed, each member of the set being stored in the form of the:

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4155

Fig. 75. Two intermediate approximations with almost the same relative error. (a) p ˆ 1 without any special functions; (b) p ˆ 1 with the addition of pcorner ˆ 1 and pvoid ˆ 1, both at nlayers ˆ 0.

1. Geometry of the handbook domain: The handbook problem domain consists of an outer boundary surrounding one or more interior boundaries. The handbook problem domain must be chosen suf®ciently large such that the effect of the outer boundary on the usable part of the handbook solutions is negligible, such that it encompasses the required support of the handbook functions in all of their occurences. 2. Approximation parameters: The GFEM approximations include a small set of parameters which specify the mesh, and the GFEM basis used in the numerical solution of the handbook problems. The approximation parameters can be chosen adaptively to produce suf®ciently accurate handbook solutions. 3. Matrix factorization: The factorization of stiffness matrix corresponding to the mesh and and the GFEM basis. This factorization is used to generate the handbook functions at run-time by applying a series of Neumann boundary conditions, corresponding to $R…zi †
n and $I…zi †
n, i ˆ 1; . . . ; q, where q is order of the handbook functions to be generated. This process is trivially fast for numerically generated handbook solutions with number of degrees of freedom N less than 10 000. For handbook solutions using with N  10 000, it is more ef®cient to pre-compute and store the solution coef®cients for each handbook function. The computation of a GFEM solution utilizing the library of precomputed handbook functions proceeds as follows: Algorithm 2. Solve a boundary value problem using the GFEM and a library of numerically generated handbook solutions. De®ne main problem Place handbooks

Determine the domain geometry, boundary conditions, and mesh for the boundary value problem. Determine, from the description of the domain geometry, which handbooks are needed, and for each occurrence of each handbook, its position, scaling, and rotation, and the degree q of the handbook functions.

4156

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 76. An example used in the veri®cation of the implementation of the numerically generated handbook solutions. (a) The domain and employed GFEM mesh (h ˆ 1=2 element size); (b) the handbook domain and mesh (h ˆ 1=4); (c) placement of the handbook mesh on top of the domain.

Generate handbooks Create PUM

For each handbook used, compute the handbook functions up to the required degree q. Assign the appropriate handbook functions at vertices of the mesh at nlayers around each feature.

5.2. Validation of the implementation of the numerically generated handbook functions The implementation of the GFEM with numerically generated handbook functions is rather complex, and it is necessary to design several example problems for checking its correctness. Let us give an example of such a veri®cation problem. Example 13 (Numerically generated handbook functions for an elliptical void feature). Here we constructed an example which checks the implementation of numerically generated handbook functions by taking advantage of the known analytical expressions of the harmonic handbook functions for elliptical voids. We employed the Neumann boundary value problem for the Laplacian in the square domain shown in Fig. 76(a), with Neumann boundary conditions corresponding to the exact solution. Process

Assemble the discrete system of GFEM equations by integrating sti€ness and load vector on elements of the GFEM mesh. Drop integration points onto the appropriate (translated/scaled/rotated) handbook mesh for evaluation of the handbook functions as necessary, employing the inverse mapping from physical coordinates to master elements of the handbook mesh. Solve the system of GFEM equations. Also use the handbook solutions in the basis for the least-square recovery estimator, if needed, at a higher order than that used in the assemble/solve step.

Postprocess

uEX ˆ

2 X

ak R…zk ‡ z k † ‡ bk I…zk

z k †;

…5:1†

kˆ1

where z is the conformally mapped coordinate from the ellipse to the unit circle, ak ˆ f1; 1g, and bk ˆ f1; 1g, and uEX satis®es the homogeneous Neumann boundary condition $uEX
n ˆ 0, on the boundary of the elliptical void. For this problem, we computed a GFEM solution uGFEM by employing a bi-p FE basis of degree p ˆ 1 on the mesh shown in Fig. 76(b) with numerically generated handbook functions of degree q ˆ 2 for the void at all vertices of the mesh. The handbook solutions, u~GFEM;i ; i ˆ 1; . . . ; 4, for q ˆ 2, were computed by imposing boundary conditions for the series of four handbook problems are the pure Neumann boundary conditions

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4157

Table 9 A series of uGFEM solutions for the main problem resulting from using bi-p FE basis with q ˆ 2 numerical handbook solutions for the elliptical void. The computed value of the energy norm jjuGFEM jjU for various combinations of meshsizes for the handbook problem h and meshsizes for the main problem h are shown Handbook problem

Global problem

Meshsize, h

Degrees of freedom, N

Meshsize, h

Degrees of freedom, N

jjuGFEM jjU

1 1 1/2 1/2 1/4 1/4

20 20 45 45 125 125

1 1/2 1 1/2 1 1/2

20 45 20 45 20 45

31.854018731881 31.854018770962 31.854018730056 31.854018730121 31.854018730056 31.854018729899

jjuEX jjU ˆ

31.854018730346

gi ˆ

o  ai R…zi ‡ z i † ‡ bi I…zi on

z i †



 on C

…5:2†

with a1 ˆ f1; 0g, and b1 ˆ f0; 0g; a2 ˆ f0; 1g, and b2 ˆ f0; 0g; a3 ˆ f0; 0g, and b3 ˆ f1; 0g; a4 ˆ f0; 0g, and b4 ˆ f0; 1g, and by employing the handbook domain and the mesh shown in Fig. 76(c), using the bip FE basis of degree p ˆ 1, with the analytically known special void functions pvoid ˆ 2, added at all vertices of the handbook mesh. Table 9 gives the values of the energy of the computed solution jjuGFEM jjU for various choices of the mesh size for the meshes for the global problem and the handbook problems. As expected, for the employed choice of approximation parameters we obtain jjuGFEM jjU ˆ jjuEX jjU with the number of correct signi®cant digits depending on integration tolerance and round-off error. 5.3. Examples illustrating the use of numerically generated handbook functions We will now give several examples of GFEM solutions which employed numerically generated handbook functions. Example 14 (GFEM solutions of the simple Neumann Problems using numerically generated square void functions). We computed eight handbook functions for the square void by solving the Neumann problems for the Laplacian in the master domain shown in Fig. 77(a), with Neumann boundary conditions corresponding to the real and the imaginary part of zi , for i ˆ 1; . . . ; 4, namely ou ˆ $…R…zi ††
n on

…5:3a†

ou ˆ $…I…zi ††
n on

…5:3b†

and

applied on the outer boundary, and zero Neumann boundary condition, ou=on, applied on the inner boundaries. The solutions were computed using the p-method with a variable p degree approximation on the GFEM I mesh with degree p1 employed in all the elements except those adjacent to the outer boundary of the master square in which the higher degree p2 ˆ p1 ‡ 2 was employed, as shown in Fig. 78. In addition the vertex corner functions of degree pcorner ˆ 1 were added to the approximation at nlayers ˆ 1 around each re-entrant corner. Figs. 79 and 80 show the contour plots of the relative modulus of the gradient, and the value of the energy norm for the pairs of the numerically generated handbook functions computed by

4158

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 77. The boundary value problems employed in the numerical generation of the handbook functions for a square void. (a) The square void in the employed master domain; (b) and (c) the employed Neumann boundary conditions for q ˆ 1, and q ˆ 2, respectively.

Fig. 78. The mesh, and the assignment of degree p of the elements, used to solve numerically the handbook problems for the square void in a master square domain. The depicted mesh was generated from the coarse mesh of 20 elements, the boundaries of which are shown with the thick lines, by re®ning the elements with a vertex at a re-entrant corner 5 consecutive times.

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4159

Fig. 79. The shades of the computed gradient for the numerically generated handbook functions for the square void. (a) i ˆ 1, and (b) i ˆ 2.

employing p1 ˆ 4, for i ˆ 1; 2 and 3, 4, respectively, while Table 10 reports the p convergence of the energy norm of the computed handbook solutions, for p1 ˆ 1; . . . ; 4. From these results it can be seen that the energy norm of the handbook solutions obtained using p1 ˆ 4, has converged to at least six digits, and the relative value of the energy norm of the error for the solution computed using p1 ˆ 3, is of the order of 1%.

4160

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Table 10 p-Convergence of the energy norm of the GFEM I approximation of the handbook functions for the square void computed in the master domain and the mesh shown in Fig. 78, obtained by employing the variable bi-p FE basis with added vertex corner functions of degree pcorner ˆ 1 at nlayers ˆ 0 around each re-entrant cornera p1

1 2 3 4 %E a

N

316 904 1868 3208

jjuGFEM jjU R…z†

I…z†

R…z2 †

I…z2 †

R…z3 †

I…z3 †

R…z4 †

I…z4 †

32.135082 32.137038 32.137110 32.137123

32.135082 32.137038 32.137110 32.137123

418.05196 418.05597 418.05597 418.05597

418.05196 418.05597 418.05597 418.05597

3230.9671 3230.9686 3230.9686 3230.9686

3230.9671 3230.9686 3230.9686 3230.9686

18087.847 18089.654 18089.655 18089.655

18087.847 18089.654 18089.655 18089.655

<0.090%

<0.090%

<0.0068%

<0.0034%

<10 4 %

<10 4 %

<10 5 %

<10 5 %

%E is the percentage relative di€erence in the energy norm of the last two solutions in the sequence for each problem.

Further, we constructed GFEM approximations of the solution of Problems A and C, which were de®ned in Part I, on a coarse uniform mesh, using the bi-p FE basis of degree p ˆ 1; 2; 3; employing the numerically generated handbook functions of degree q ˆ 1, at nlayers ˆ 0, around each square void. Fig. 81 illustrates the placement of the handbook meshes over the coarse GFEM II mesh for Domain I, while Fig. 82 indicates the elements over which the void functions are integrated. The placement of the vertex square void functions used in these solutions are shown in Fig. 83. Figs. 84 and 85 show the contour plots of the relative modulus of the gradient, and the value of the energy norm of the solution for the GFEM II and GFEM III approximations of Problem A, respectively, and Fig. 86 shows the corresponding results for the GFEM III approximation of Problem C. Note that for the Problem A the GFEM with numerically computed square void functions achieves accuracies of 13%, 4%, and 2%, respectively, by employing p ˆ 1, 2, and 3, which translates into only 221, 481, and 893, degrees of freedom for the global GFEM approximation, respectively. Similar results were obtained for GFEM III approximations of Problem A and C. Let us note that the dramatic improvement in the accuracy measured by the relative value of the energy norm of the error is because the main contribution to the energy is near the re-entrant corners of the square voids. The main conclusion from this example is that: · The constructed implementation of the GFEM with numerically generated handbook functions can be used for general complex domains, and with a judicious choice of the handbook functions a dramatic improvement in the accuracy can be obtained. Example 15 (GFEM solutions of the simple mixed problems using both the analytically known reentrant corner functions and numerically generated square void functions). We constructed GFEM approximations of the solution of Problems B and D (see Fig. 22), using the bi-p FE basis of degree p ˆ 1±3 with the numerical handbook functions of degree q ˆ 1 for each square void at nlayers ˆ 0 around each void, and the re-entrant corner functions of degree pcorner ˆ 1, at nlayers ˆ 0, around each re-entrant corner of the domain not part of any of the square voids. The employed mesh and the placement of the vertex square void functions are shown in Fig. 87. Figs. 88 and 89 show the contour plots of the relative modulus of the gradient, and the value of the energy norm of the solution for the GFEM II and alternate GFEM II approximations of Problem B, respectively, and Fig. 90 shows the corresponding results for the GFEM II approximation of Problem D. Note that for the Problem B we have achieved accuracies of 10%, 2.1%, and 0.3%, by employing p ˆ 1, 2 and 3, respectively, by employing only 307, 744, and 1449 degrees of freedom in the global GFEM

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4161

Fig. 80. The shades of the computed gradient for the numerically generated handbook functions for the square void. (a) i ˆ 3, and (b) i ˆ 4.

approximation, respectively. Again, similar results were obtained for other approximations of Problem B and D. Once more, we conclude that: · By employing a judicious choice of the employed numerically generated and analytically known handbook functions, we can obtain good accuracy using rather coarse meshes.

4162

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 81. The handbook problem master squares and meshes, shown in light gray, placed on top of the Domain I.

Example 16 (The h-adaptive control of the error in a GFEM approximation using numerically generated handbook functions). We considered Problem A and computed GFEM II and GFEM III, h-adaptive sequence of approximate solutions by employing only the bilinear p ˆ 1, FE basis on the approximation mesh with added numerically generated vertex square void functions of degree psquare ˆ 1, at nlayers ˆ 0, around each square void. We employed the least-squares recovery estimator, using a vertex patch basis of the harmonic polynomials of degree p0 ˆ 2, and also the square void functions of degree psquare ˆ 2 at nlayers ˆ 0, and we used this indicator to drive the feedback h-adaptivity, starting from the coarse GFEM II and GFEM III meshes shown in Fig. 91(a), and (b), respectively. Fig. 92(a) shows the convergence of the relative error for the adaptive mesh sequences, and Fig. 92(b) shows the corresponding values of the global effectivity index. Fig. 93 shows two intermediate meshes obtained during the sequence with almost the same relative error.

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4163

Fig. 82. The handbook problem meshes placed on top of the GFEM II mesh employed for the solution of Problem A. Here, the outer layer of elements in each handbook mesh was omitted for clarity. The elements over which the square void functions are supported are shown as shaded.

Fig. 83. The vertex locations of the numerically generated square void functions on the coarse GFEM meshes used to solve the simple Neumann problems, for (a) The GFEM II approximation employed for Problem A; (b) The GFEM III approximation employed for Problem A; (c) The GFEM III approximation employed for Problem C.

4164

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 84. Relative modulus of the gradient computed by a GFEM II approximation of Problem A. These approximations were computed by employing the bi-p FE basis of degree (a) p ˆ 1, (b) p ˆ 2, and (c) p ˆ 3, with the numerically generated square void functions of degree q ˆ 1, added at the vertices shown in Fig. 83(a).

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4165

Fig. 85. Relative modulus of the gradient computed by a GFEM III approximation of Problem A. These approximations were computed by employing the bi-p FE basis of degree (a) p ˆ 1, (b) p ˆ 2, and (c) p ˆ 3, with the numerically generated square void functions of degree q ˆ 1, added at the vertices shown in Fig. 83(b).

4166

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 86. Relative modulus of the gradient computed by a GFEM III approximation of Problem C. These approximations were computed by employing the bi-p FE basis of degree (a) p ˆ 1, (b) p ˆ 2, and (c) p ˆ 3, with the numerically generated square void functions of degree q ˆ 1, added at the vertices shown in Fig. 83(c).

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4167

Fig. 87. The vertex locations of the numerically generated square void functions, and the vertex corner functions on the coarse GFEM meshes used to solve the simple mixed problems. (a) The GFEM II mesh employed for Problem B; (b) Another GFEM II mesh employed for Problem B; (c) The GFEM II mesh employed for Problem D.

We note that: · The capabilities of a posteriori estimation and adaptive control of the error, which have been developed for the standard FEM, can be directly extended to the most general version of the GFEM with numerically generated handbook functions. Example 17 (The effect of higher degree numerically generated handbook functions on the error). In this example, we illustrate the e€ect of higher degree numerically generated handbook functions on the accuracy of the approximation in solving Problem C using a GFEM III approximation, and the pointwise e€ectivity of the GFEM least-squares recovery. In this example, in addition to the numerically generated functions for the square voids, we also generated and used handbook functions for the two polygonal voids that form boundaries C2 and C3 for use in solving Problem C. Using the p-method, as we did for the square void problem, with the meshes shown in Figs. 94 and 95, we solved a series of eight Neumann problems for the Laplacian in the surrounding square master domain with Neumann boundary conditions corresponding to the real and the imaginary part of zi , for i ˆ 1; . . . ; 4. Again, we used the variable p degree approximation with degree p1 ˆ 1; . . . ; 4; employed in all the elements except in the two layers adjacent to the outer boundary where we employed degree p2 ˆ p1 ‡ 2, and the vertex corner functions of degree pcorner ˆ 1, at nlayers ˆ 1 around each re-entrant corner. Figs. 96 and 97 show the contour plots of the relative modulus of the gradient, and the value of the energy norm of the pairs of the numerically generated handbook functions obtained by employing p1 ˆ 4. Tables 11 and 12 report the values of the energy norm of the computed solutions jjuGFEM jjU , for p1 ˆ 1; . . . ; 4, for the handbook problems for C2 and C3 . From these results it can be seen that the energy norm of the handbook solutions obtained using p1 ˆ 4 has converged to at least six digits. We computed GFEM III approximations of Problem C using the numerically generated square void functions, and the two polygonal void functions for boundaries C2 and C3 placed as shown in Figs. 98 and 99. We used the bi-p FE functions of uniform degree p ˆ 1, or 2, and added the numerically generated functions of degree q ˆ 1, or 2 at nlayers ˆ 0, around the square and polygonal voids. We also computed the pointwise error, by subtracting these solutions from an overkill solution, which was computed on the same mesh with one additional uniform re®nement and six nested re®nements of

4168

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 88. Relative modulus of the gradient computed by a GFEM II approximation of Problem B. We employed the bi-p FE basis of degree (a) p ˆ 1, (b) p ˆ 2, and (c) p ˆ 3, with the re-entrant corner functions of degree pcorner ˆ 1, and the numerically generated square void functions of degree q ˆ 1, added at the vertices of the mesh as shown in Fig. 87(a).

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4169

Fig. 89. Relative modulus of the gradient computed by an alternate GFEM II approximation of Problem B. We employed the bi-p FE basis of degree (a) p ˆ 1, (b) p ˆ 2, and (c) p ˆ 3, with the re-entrant corner functions of degree pcorner ˆ 1, and the numerically generated square void functions of degree q ˆ 1, added at the vertices of the mesh as shown in Fig. 87(b).

4170

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 90. Relative modulus of the gradient computed by a GFEM II approximation of Problem D. We employed the bi-p FE basis of degree (a) p ˆ 1, (b) p ˆ 2, and (c) p ˆ 3, with the re-entrant corner functions of degree pcorner ˆ 1, and the numerically generated square void functions of degree q ˆ 1, added at the vertices of the mesh as shown in Fig. 87(c).

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4171

Fig. 91. The initial coarse meshes used in the adaptive h-convergence employing the numerically generated vertex square void functions. (a) Initial GFEM II mesh with 76 elements; and (b) initial GFEM III mesh with 256 elements. Table 11 p-Convergence of jjuGFEM jjU , the energy norm of the GFEM I approximation of the handbook functions for the polygonal void (C2 of Problem C) computed in the master domain and the mesh shown in Fig. 94 by employing the variable bi-p FE basis with added vertex corner functions of degree pcorner ˆ 1 at nlayers ˆ 0a p1

1 2 3 4

N

780 2060 4100 6900

%E a

jjuGFEM jjU R…z†

I…z†

R…z2 †

I…z2 †

14.228930 14.232471 14.232522 14.232523

14.238259 14.241907 14.241958 14.241959

77.263590 77.273518 77.273543 77.273544

77.263114 77.273017 77.273042 77.273042

<0.0313%

<0.0348%

<0.0107%

<0.0098%

%E is the percentage relative di€erence in the energy norm of the last two solutions in the sequence for each problem.

Table 12 p-Convergence of jjuGFEM jjU , the energy norm of the GFEM I approximation of the handbook functions for the polygonal void (C3 of Problem C) computed in the master domain and the mesh shown in Fig. 95 by employing the variable bi-p FE basis with added vertex corner functions of degree pcorner ˆ 1 at nlayers ˆ 0a p1

1 2 3 4 %E a

N

546 1442 2870 4830

jjuGFEM jjU R…z†

I…z†

R…z2 †

I…z2 †

10.200325 10.203318 10.203431 10.203432

10.184354 10.187199 10.187286 10.187287

59.280167 59.293479 59.293854 59.293858

59.300688 59.312855 59.313226 59.313230

<0.0441%

<0.0441%

<0.0368%

<0.0367%

%E is the percentage relative di€erence in the energy norm of the last two solutions in the sequence for each problem.

4172

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 92. The adaptive h-convergence of GFEM II and GFEM III on Problem A using the computed indicators, showing the effect of the numerically generated square void functions versus the vertex corner functions. (a) The log±log graph of the relative error versus the number of degrees of freedom. (b) The effectivity index versus the logarithm of the number of degrees of freedom when numerically generated square void functions are used.

elements intersecting the re-entrant corners; for this overkill approximation we used the bi-p FE basis with degree p ‡ k, where k ˆ 4, with re-entrant corner functions of degree pcorner ˆ 1 at nlayers ˆ 0 added around each re-entrant corner. Figs. 100(a)±(d) show the contours of the relative modulus of the gradient of the error for combinations of bi-p FE basis of degree p ˆ 1 and 2 with numerically generated handbook functions of degree q ˆ 1 and 2. Note that the effect on the error of the higher degree q ˆ 2 versus the q ˆ 1 handbook functions while using the same degree FE bi-p basis is negligible, however, there is dramatic improvement in the accuracy using the handbook functions from the bilinear p ˆ 1 to the biquadratic p ˆ 2 case. In order to illustrate the use of the handbook solutions in the a posteriori estimation of the error, we computed the estimated error eREC GFEM for the ®rst uGFEM above, namely the one with p ˆ 1, q ˆ 1, using the harmonic monomial basis of degree p ˆ 2, and all the q ˆ 2 handbook solutions. Figs. 101(a)±(d) show the

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4173

Fig. 93. Two intermediate approximations of the h-adaptive sequences of meshes using bi-p FE basis of degree p ˆ 1 and the numerically generated vertex square void functions, psquare ˆ 1, at nlayers ˆ 0. (a) GFEM II mesh; and (b) GFEM III mesh.

4174

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 94. The mesh, and the assignment of degree p of the elements, used to numerically solve the handbook problems for the polygonal void (C2 of Problem C) in a master square domain. The mesh was generated from the coarse mesh of 20 elements, the boundaries of which are shown with the thick lines, by re®ning the elements with a vertex at a re-entrant corner 5 consecutive times.

Fig. 95. The mesh, and the assignment of degree p of the elements, used to numerically solve the handbook problems for the polygonal void (C3 of Problem C) in a master square domain. The mesh was generated from the coarse mesh of 14 elements, the boundaries of which are shown with the thick lines, by re®ning the elements with a vertex at a re-entrant corner 5 consecutive times.

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4175

Fig. 96. The shades of the computed gradient for the numerically generated handbook functions for the polygonal void (C2 ), i ˆ 1, and i ˆ 2.

4176

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 97. The shades of the computed gradient for the numerically generated handbook functions for the polygonal void (C3 ), i ˆ 1; and i ˆ 2:

contours of the relative modulus of the gradient of the estimated error for the recovery and the global e€ectivity index of the estimator jREC . Note that the least-squares recovery estimator incorporating the handbook solutions is accurate, and that the contour plots of the estimated error are similar to the exact error previously shown in Fig. 100.

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4177

Fig. 98. The handbook problem master squares and meshes, shown in light gray, placed on top of Domain II.

In conclusion, we note that: · Our implementation of the GFEM with numerically generated handbook functions is robust enough to be used for many combinations of numerically generated and analytically known handbook functions, including the capability of error estimation. · We cannot say a priori which handbook functions will be most e€ective when added to the approximation of the solution and its error. Example 18 (The versatility of the GFEM with numerically generated handbook functions). In this example, we further underline the versatility of the numerically generated handbooks for solving boundary value problems with complex geometries. As an example we considered the case where several bifurcated cracks may be added at arbitrary locations in the problem domain.

4178

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 99. The handbook problem meshes placed on top of the GFEM III mesh employed for the solution of Problem C. Here, the outer layers of elements in each handbook mesh was omitted for clarity.

Table 13 p-Convergence of jjuGFEM jjU , the energy norm of the GFEM I approximation of the handbook functions for the bifurcated crack computed in the master domain and the mesh shown in Fig. 103 by employing the variable bi-p FE basis with added vertex corner functions of degree pcorner ˆ 1 at nlayers ˆ 0a p1

N

jjuGFEM jjU R…z†

1 2 3 4

733 1590 2867 4564

5.0113269 5.0118357 5.0118581 5.0118628

%E a

<0.137%

I…z† 5.0765512 5.0798867 5.0800088 5.0800637 <0.464%

R…z2 †

I…z2 †

10.206048 10.206562 10.206564 10.206564

10.207923 10.208039 10.208046 10.208050

<0.0287%

<0.0822%

%E is the percentage relative di€erence in the energy norm of the last two solutions in the sequence for each problem.

We will be interested in solving a Neumann problem on Domain I with containing four large bifurcated cracks as shown in Fig. 102(a), which we will call Problem M, with zero Neumann boundary conditions imposed on the internal boundaries of the cracks, namely: Problem M: Du ˆ 0 on ou ˆ $…2x on ou ˆ 0 on on ou ˆ 0 on on

X; y†
n

…5:4a† on C1 ;

…5:4b†

C2 ; C3 ; . . . ; C10 ;

…5:4c†

C11 ; . . . ; Cncracks :

…5:4d†

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4179

Fig. 100. Relative modulus of the gradient of the error in four GFEM III approximations of Problem C showing the e€ect of the degree of the handbook functions, we used square void functions, and polygonal void functions for C2 and C3 . The approximations using (a) p ˆ 1, q ˆ 1, and (b) p ˆ 1, q ˆ 2, (c) p ˆ 2, q ˆ 1, and (d) p ˆ 2, q ˆ 2.

4180

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 101. Relative modulus of the gradient of the estimated error eREC GFEM in the four GFEM III approximations of Problem C showing the e€ect of the degree of the handbook functions, we used square void functions, and polygonal void functions for C2 and C3 . The approximations using (a) p ˆ 1, q ˆ 1, and (b) p ˆ 1, q ˆ 2, (c) p ˆ 2, q ˆ 1, and (d) p ˆ 2, q ˆ 2.

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4181

Fig. 102. Problem M, a modi®ed form of Domain I including four bifurcated cracks at various locations. (a) The domain and applied Neumann boundary conditions; (b) contours of the modulus of the computed overkill gradient.

Fig. 103. The mesh, and the assignment of degree p of the elements, used to numerically solve the handbook problems for the bifurcated crack in a master square domain. The mesh was generated from the coarse mesh of 15 elements, the boundaries of which are shown with the thick lines, by re®ning the elements with a vertex at a re-entrant corner 5 consecutive times.

4182

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 104. The shades of the computed gradient for the numerically generated handbook functions for the bifurcated crack, i ˆ 1; and i ˆ 2:

We computed four handbook functions (real and imaginary parts of the sequence i ˆ 1, 2, as in the three previous examples) for a bifurcated crack on the mesh on the master square as shown in Fig. 103 (see also Fig. 104). Table 13 reports the values of the energy norm of the computed solutions kuGFEM kU , for p ˆ 1; . . . ; 4, to the handbook problems for the bifurcated crack. Figs. 105 and 106 illustrate the placement of the handbook meshes over the coarse GFEM II mesh for Problem M. We computed an overkill solution

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4183

Fig. 105. The handbook problem master squares and meshes for the bifurcated crack, shown in light gray, placed on top of the modi®ed Domain I.

on a GFEM II mesh with three uniform nested re®nements using a bi-p FE basis of degree p ˆ 2 with handbook functions for the bifurcated crack of degree q ˆ 2 at nlayers ˆ 2 around these features and reentrant corner functions of degree pcrack ˆ 1 at nlayers ˆ 2 around all the re-entrant corners other than the cracks. Fig. 102(b) shows this overkill mesh and the energy norm and contours of the modulus of the gradient of the overkill solution. Using the overkill solution we computed the energy norm of the error in two GFEM approximations that used the bi-p FE basis p ˆ 1, and 2, respectively; both approximations used the re-entrant corner functions of degree pcrack ˆ 1 at nlayers ˆ 0 around the re-entrant corners, and the numerically generated bifurcated crack functions of degree q ˆ 1 at nlayers ˆ 0 around each bifurcated crack. Figs. 107(a) and (b) show the contours of the modulus of the computed gradient of these two p ˆ 1 and 2 solutions, the error in the energy norm was 19% and 6%, respectively. Figs. 107(c) and (d) show the contours of error eGFEM in the computed gradient with respect to the overkill solution. Once more, we see that the GFEM with numerically generated handbook functions mixed with a bi-p FE basis of degree p ˆ 1 or 2 can be used to achieve engineering accuracy using a minimal number of degrees of freedom.

4184

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 106. The handbook problem meshes placed on top of the domain employed for the solution of Problem M. Here, the outer layer of elements in each handbook mesh was omitted for clarity.

We conclude: · The use of the numerically generated handbook for a complex geometrical domain features can provide accurate solutions when it becomes dicult to construct a traditional ®nite element mesh that conforms to the problem domain. 6. Computer implementation of the generalized ®nite element method The GFEM coupled with a library of numerically generated handbook functions involves the manipulation of multiple load cases, multiple meshes and domains, and the many mapping on and between the elements and domain of each. In order for the software implementation of such a complex system to exhibit maintainability and allow productive work to be accomplished, an object model design of the software architecture is necessary. As a means to help manage such complexity, the object model addresses the organization of the major conceptual entities, for example, see [26]. In this section, we outline how we used the object model to organize and implement the GFEM, mainly the parts of our code directly applicable to the numerically generated handbook functions, and also the analytically known handbook functions. Our code is called GFEM2D and it was written speci®cally to adhere to the Fortran 95 standard, which is an object-based language, rather than a fully object-oriented language such as C++, in the sense that it allows user-de®ned classes or types, and with the notion of modules it allows data hiding or encapsulation, the association of functions with types, and function overloading. Fortran 95 does not, however, allow inheritance or strict polymorphism, see [26,27]. We selected Fortran for its above-mentioned object-based extensions, since they were added to a language we were already familiar with, and also because of its portability, its powerful array expression capability, and its maturity, which results in its speed of compilation and execution.

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4185

Fig. 107. GFEM II approximations of Problem M using bi-p FE functions of degree p ˆ 1 and p ˆ 2 and the numerically generated handbook function for the bifurcated crack at nlayers ˆ 0. The relative modulus of the gradient of the computed solutions for the cases (a) p ˆ 1, and (b) p ˆ 2, and also the relative modulus of the error in the gradient with respect to the overkill solution for the same cases (c) p ˆ 1, and (d) p ˆ 2.

4186

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 108. The data structures relating to the de®nition of the approximation on a mesh of elements. The main type is the mesh, de®ned in MESH_DATA_MODULE, which has relationships with types in the DOMAIN_MODULE (see Fig. 109), PUM_DATA_MODULE (see Fig. 112), BASIS_MODULE, and QUANTITY_MODULE.

The numerically generated functions in the GFEM is the use of numerical solutions of one or more boundary value problems in the approximating basis of another main problem. This is necessarily made easier to implement by treating each conceptual entity of the FEM and GFEM as a data object with associated method, including, but not limited to the approximation mesh of elements, the problem domain and the individual geometrical features making up the domain, the partition of unity (PU) and special families of vertex functions, and the element integration mesh. We will examine these major components and how they are implemented as classes for objects in our code GFEM2D. The largest composition in GFEM2D is our class representing the mesh, MESH_T. Fig. 108 shows a uni®ed modeling language (UML), see [28], class diagram of MESH_T and its relationship with other major classes, and their containing modules of GFEM2D. The MESH_T responsibilities are to:

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4187

1. collect and catalog all connections of the standard FE entities, including elements, nodes; 2. maintain the geometric relationship of mesh entities to a speci®c problem domain, X, the class DOMAIN_T, and using the element mappings SUPER_ELEMENT_T, these attendant classes are explained below; 3. keep track of the approximation basis, including the standard bi-p FE basis functions in BASIS_T, and also the PUM special function basis (both analytically known and numerically generated handbooks) associated with speci®c vertices of the mesh in PUM_MESH_T. An object of MESH_T is created, or generated, on top of a problem DOMAIN_T by the group of parameters in a MESH_GENERATE_T. The main purpose in the life of a mesh is to control the error in a group of engineering quantities of interest, to this end, the MESH_T collaborates with multiple instances of the class QUANTITY_OF_INTEREST. We will not detail this latter class in this paper. For a single execution run of GFEM2D, one object of the type MESH_T exists for each boundary value problem, i.e., one MESH_T for the main problem, and one MESH_T for each numerically generated handbook problem that is used in the main problem. The numerically generated handbook meshes are read from storage on disk; note only one MESH_T exists for each handbook, even if the handbook will be used in multiple occurrences. All the MESH_T objects coexist concurrently in memory as the code proceeds to integrate and solve the GFEM system of equations. The main procedures of GFEM2D all basically take a MESH_T as an input argument, making it easy to perform tasks without regard for whether we are operating on a handbook mesh/domain or on the main problem mesh/domain. Let us also give details of the object model comprising the problem domain. Fig. 109 shows the UML class diagram of the data structures de®ning the domain entities, including what we term the super-elements, super-edges, the features (which are geometric sections of the domain boundary), and their relationships. We will de®ne these objects further below. The domain type owns both an array of super-elements and an array of super-edges, the elements(:) and edges(:), respectively. The domain type also keeps track of the re-entrant corners, or singularities by storing an array of domain singularity type. Additionally, the domain stores the name of a set of materials, the MATERIAL_SET_T, and the name of a set of geometric features, the FEATURES_T. The class DOMAIN_T acts as a stand-in for the mathematical domain X of the boundary value problem is the two-dimensional region on which the partial di€erential equation is solved, and is represented by a data structure which collects together the geometric entities used to de®ne the domain, its boundary conditions, and a coarse mesh of elements on top of X. The fundamental entities used to de®ne a domain are the super-nodes and super-edges. The supernodes and super-edges are grouped together in such a way as to de®ne both the boundary of the domain (along with boundary conditions), and the coarsest starting mesh for the problem. The starting mesh of coarse elements is called the set of super-elements. Groupings of the super-edges are called features. Examples of de®ning the domain using these entities and their interactions are explained in this section. · A super-node is a point-like entity with position being its only attribute. A super-node serves the purpose of acting as a vertex for two or more adjoining super-edges; additionally a super-node may form the vertex of one or multiple quadrilateral super-elements. · A super-edge is a line or curve-like entity in two dimensions. It de®nes one of two types of boundaries: (1) the boundary of the domain or (2) the shared boundary of two super-elements. In the ®rst case ± (1) the boundary of the domain ± the super-edge may or may not also form the boundary of a single super-element. Super-edges occur singly or grouped sets. Within sets, the super-edges are connected end-to-end into a closed loop by shared super-nodes. · A super-element is the coarsest ®nite element possible, i.e., an unre®ned element. A super-element is a four sided region in two dimensions. A quadrilateral super-element is always de®ned by four super-nodes at its vertices and four surrounding super-edges, these nodes and edges may or may not be shared by other super-elements. The set of super-elements acts as the initial ®nite element mesh, they also act as the mapping from master ®nite element coordinates to physical coordinates. Material properties are constant inside the super-element. · A feature is a domain or material boundary not coincident with the boundary of super-elements. A feature is de®ned by a set of one or more super-edges. The set may be one single super-edge, or the set may

4188

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 109. The data types relating to the de®nition of the domain. These include mainly the types of the DOMAIN_MODULE and their relationship with the types in the FEATURES_DATA_MODULE and the MATERIAL_MODULE.

be multiple grouped super-edges forming a loop. There are di€erent kinds of features, each has a name in the code as shown in Fig. 110. Each super-edge has an orientation, with respect to the two super-nodes; the super-edge points from the ®rst to the second super-node. If grouped into a closed loop, all super-edge orientations continue either clockwise or counter-clockwise. The domain is de®ned to lie always to the left of such looped super-edges. Thus the outermost boundary must necessarily be oriented counter-clockwise, further, counter-clockwise oriented feature loops always enclose the domain, while clockwise-oriented loops (and also feature ellipses) encompass a void in the domain. See again Fig. 110. Non-grouped (single) super-edges act either as the boundary between two super-elements (with super-nodes at both ends), or as the only de®ning edge of a feature. Only in this very last role may the super-edge exist without super-nodes at its end points. Fig. 111 show an enumeration of super-nodes and geometric features de®ning a domain surrounded by a single

c

Fig. 111. The computer representation of a typical GFEM III domain showing the enumeration of super-nodes and geometric features. This con®guration of super-nodes, super-edges, super-elements, and features results from reading an input ®le for the boundary problem description. Note the code's enumeration of the super-nodes and the features respectively; the outer boundary is feature 36, and the two largest polygonal voids are features 37 and 38.

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4189

Fig. 110. The kinds of features used to de®ne boundaries in the domain and their parameter names used as integer constants in the Fortran code.

4190

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

Fig. 112. The data types relating to the de®nition of the special functions at element vertices. These are the main types contained in the PUM_DATA_MODULE.

Fig. 113. The data types relating to the de®nition of the integration cell mesh. These are the basic types contained in the ICELL_MODULE and the two mesh types in ICCELL_MESH_MODULE.

square super-element, thus serving as the start of a GFEM III approximation mesh. Finer meshes of elements would be created inside this one super-element. The class which is important for the necessary book-keeping of di€erent families of special vertex functions arising from various geometric features is the PUM_MESH_T class. Fig. 112 shows the UML class diagram for the PUM_MESH_T class and its composition and collaborators. The PUM_MESH_T is always associated with only one MESH_T, but one PUM_MESH_T may be used for the basis of the approximation and another PUM_MESH_T will be used to describe the special functions in the recovery error estimator for an approximation. The parameters describing the degrees of special vertex functions: pvoid , pcorner , pcrack , or q for numerically generated functions, the number of layers nlayers is given in PUM_MESH_READ_T; these parameters are used as input for creating the PUM_MESH_T associated with a given MESH_T for a problem. For each element of the approximation MESH_T, we create a PUM_ELEMENT_T describing the number of

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4191

special vertex functions with support in the element and their required integration tolerances during the assembly of the element equations. Similarly, for each mesh vertex (node), we have an associated PUM_VERTEX_T which describes: 1. the families, or kinds, of special functions associated with the associated vertex; 2. the number of functions in each family, or equivalently, the degree of the functions; 3. the geometric feature giving rise to the family of functions. The last major group of classes in the GFEM2D design, is the integration mesh, which is a mesh of quadrature cells we denote as ICELL_MESH_T. The details of the elementwise Fast Remeshing algorithm for numerical integration were previously given in Section 3.1. Fig. 113 shows the UML class diagram for ICELL_MESH_T and the hierarchy of the other attendant classes of objects. The classes ICELL_T, IGROUP_T, and INODE_T actually form the nested tree of the integration mesh of quadrature cells. The IQUAD_MESH does the actual work of performing Gauss like quadrature rules of a speci®c order inside each integration ICELL_T object. As we discussed in Section 3.1, the cells crossing the domain boundary are given special attention and are subdivided and mapped appropriately to conform exactly to the domain geometry, which is performed by the IBOUNDARY_CELL_T and its component subcells, the ISUBCELL_T. 7. Conclusions In this paper we gave an example of the design, implementation, and illustration of the performance of a partition of unity (PU)-based generalized ®nite element method (GFEM). Here, we presented the GFEM for the Laplacian in two dimensions; however, the implementation can be directly extended to more general elliptic problems in two dimensions, like for example, the orthotropic Laplacian and the elasticity problems in heterogeneous media, and with considerably more e€ort to the same problems in three dimensions. Let us recall the main attributes of the method: Geometrical flexibility. The GFEM bypasses the main dif®culty in the application of the classical FEM in engineering problems, namely the generation of a mesh as a partition of the problem domain into triangles and/or quadrilaterals. The GFEM is capable of constructing approximate solutions for the boundary value problem of interest, by employing meshes which are allowed to overlap part or all of the problem domain boundary, for example, meshes of squares generated by the nested re®nement of one square which contains the problem domain in its interior. Hybrid capability. The GFEM builds the approximation by employing a standard FE basis on the employed mesh, enriched by special or handbook functions which are incorporated into the approximation by employing the the PUM, using a PU corresponding to the employed FE mesh. This construction allows us to incorporate into the approximation local handbook or special bases which can have great impact on the achieved accuracy, while building the implementations of the method using existing FE infrastructure. The main achievements of this work, which were essential in bringing this paper to a successful conclusion are: · The implementation and the analysis of the robustness of a direct Gauss solver for the semi-de®nite, or nearly semi-de®nite systems of linear equations employed by the various PUM and/or GFEM approximations. · The design and implementation of a robust numerical integration algorithm with adaptive control of the integration error, for integrating the entries of the sti€ness matrix, the load vector, the various norms of the solution, etc., in elements which either overlap the domain boundary and/or include special functions with singular, or nearly singular behavior. 7.1. Summary of conclusions We illustrated the capabilities of the GFEM by solving boundary value problems for the Laplacian in curvilinear polygonal domains, which may also include a large number of voids, and cracks which make the construction of a classical FE mesh for the problem domain very dicult and in some cases practically impossible. From the presented illustrative examples we can make the following conclusions:

4192

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

1. The GFEM makes possible the numerical solution with high accuracy of boundary value problems in domains that cannot be meshed and solved by the classical FEM. 2. The incorporation into the approximation of analytically known or numerically computed handbook functions can greatly enhance the local and global accuracy of the computed solution. 3. The capabilities of a posteriori error estimation and adaptive control of the error which were developed for the classical FEM, can be extended to the GFEM. 4. The geometrical ¯exibility of the GFEM is enhanced, because special functions for arbitrarily complex feature geometries can be created. A mesh for a complex domain feature may be created once, used to solve a series of small problems, and the resulting solutions pasted into a larger problem, whose mesh need not conform to the feature geometry at all. The majority of the cost of engineering analysis on complex part geometries is spent in human time creating an approximation mesh, the GFEM can mitigate this cost, or even make a solution possible in settings where it was previously impossible to create satisfactory FE meshes. 5. The GFEM is enabled with the capability of using the most optimal basis functions. As we have seen, the use of the numerically created ``handbook'' functions allowed us to obtain engineering accuracies (i.e., <5% error measured in the energy norm) on complex boundary value problems ± normally requiring tens of thousands of degrees of freedom using the standard FE method ± with only hundreds of degrees of freedom. 7.2. Recommendations for further work There is great potential for extending the work on the GFEM in the following directions: · Quadrature optimization. The GFEM, especially when using the numerically generated handbook functions, moves the bulk of the computational work of the analysis from solving the algebraic system of equations, to integrating the products of the basis functions. Thus, there is great potential for further optimization of the GFEM by focusing on the numerical quadrature and the evaluation of the numerical handbook problems. The numerical quadrature in the GFEM is a perfect candidate for parallelization because the quadrature is performed element by element independently, requiring no communication across elements, the only communication is between the basis functions used and the shared domain geometry. · The GFEM for steady-state and nonlinear time-dependent elliptic problems. The GFEM can be extended to provide a powerful computational method for obtaining highly accurate solutions of nonlinear steadystate and time-dependent problems, including guaranteed a posterior estimation for the outputs of interest. The main tool in this case, is the employment of parameterized numerically generated solution dependent handbook functions. Possibilities being currently explored by other investigators are, for example, the extension of the GFEM to solid mechanics problems with large strains and the problem of crack propagation, to name a few. Three-dimensional problems. The utility of the GFEM with numerically generated handbook functions would obviously be advantageous in solving boundary problems posed in complex three-dimensional domains, where a priori expansions of solutions are not generally known for even most simple feature geometries, e.g., in elasticity. The main work in bringing the GFEM into the three-dimensional realm is in the element remeshing algorithm for the numerical integration.

Acknowledgements The authors would like to thank Dr. Luise Couchman of the Oce of Naval Research for supporting this research. Kevin Copps was primarily supported by AASERT grants from the Oce of Naval Research and the Army Research Oce. The support of Dr. Richard Lau of the ONR and Dr. Kailasam Iyer of the ARO is gratefully appreciated. T. Strouboulis and K. Copps were supported by the ONR under Grant N00014-99-1-0726; K. Copps was also supported by ARO grant DAAG55-97-1-0186. I. Babuska was supported by the NSF under grant DMS-9802367.

T. Strouboulis et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4081±4193

4193

References [1] I. Babuska, T. Strouboulis, K. Copps, The design and analysis of the generalized ®nite element method, Comput. Methods Appl. Mech. Engrg. 181 (1) (2000) 43±69. [2] T. Strouboulis, K. Copps, I. Babuska, The generalized ®nite element method: an example of its implementation and illustration of its performance, Int. J. Numer. Methods Engrg. 47 (8) (2000) 1401±1417. [3] I. Babuska, G. Caloz, J. Osborn, Special ®nite element methods in a class of second order elliptic problems with rough coecients, SIAM J. Numer. Anal. 31 (1994) 945±981. [4] J.M. Melenk, Finite element methods with harmonic shape functions for solving laplace's equation, M.S. thesis, University of Maryland, College Park, MD, 1992. [5] J.M. Melenk, On generalized ®nite element methods, Ph.D. dissertation, University of Maryland, College Park, MD, 1995 (Advisor: I. Babuska). [6] A. Duarte, The hp cloud method, Ph.D. dissertation, University of Texas at Austin, Austin, TX, 1996 (Advisor: J.T. Oden). [7] I. Babuska, J.M. Melenk, The partition of unity ®nite element method: basic theory and applications, Comput. Methods Appl. Mech. Engrg. 139 (1996) 289±314. [8] I. Babuska, J.M. Melenk, The partition of unity method, Int. J. Numer. Methods Engrg. 40 (1997) 727±758. [9] J.M. Melenk, I. Babuska, Approximation with harmonic and generalized harmonic polynomials in the partition of unity method, Comput. Assist. Mech. Engrg. Sci. 4 (1997) 607±632. [10] C.A.M. Duarte, J.T. Oden, An hp adaptive method using clouds, Comput. Methods Appl. Mech. Engrg. 139 (1996) 237±262. [11] C.A.M. Duarte, J.T. Oden, Hp Clouds ± an hp meshless method, Numer. Methods Partial Di€erential Equations 12 (1996) 673± 705. [12] T.J. Liszka, C.A.M. Duarte, W.W. Tworzydlo, Hp-meshless cloud method, Comput. Methods Appl. Mech. Engrg. 139 (1996) 263±288. [13] J.T. Oden, C.A. Duarte, O.C. Zienkiewicz, A new cloud-based hp ®nite element method, Comput. Methods Appl. Mech. Engrg. 153 (1998) 117±126. [14] T. Belytschko, T. Black, Elastic crack growth in ®nite elements with minimal remeshing, Int. J. Numer. Methods Engrg. 45 (1999) 601±620. [15] N. Moes, J. Dolbow, T. Belytschko, A ®nite element method for crack growth without remeshing, Int. J. Numer. Methods Engrg. 46 (1) (1999) 131±150. [16] N. Sukumar, N. Moes, B. Moran, T. Belytschko, Extended ®nite element method for three-dimensional crack modeling, Int. J. Numer. Methods Engrg. 48 (2000) 1549±1570. [17] J. Dolbow, N. Moes, T. Belytschko, Modeling fracture in Mindlin±Reissner plates with the extended ®nite element method, Int. J. Solids Struct. 37 (2000) 7161±7183. [18] C. Daux, N. Moes, J. Dolbow, N. Sukumar, T. Belytschko, Arbitrary branched and intersecting cracks with the extended ®nite element method, Int. J. Numer. Methods Engrg. 48 (2000) 1741±1760. [19] T. Liszka, J. Orkisz, The ®nite di€erence method at arbitrary irregular grids and its application in applied mechanics, Comput. Struct. 11 (1980) 83±95. [20] T. Liszka, An interpolation method for an irregular net of nodes, Int. J. Numer. Methods Engrg. 20 (1984) 1599±1612. [21] I. Babuska, T. Strouboulis, The Finite Element Method and its Reliability, Oxford University Press, London, 2001, to appear. [22] S.C. Brenner, L.R. Scott, in: The Mathematical Theory of Finite Element Methods, Springer, New York, 1994. [23] B.A. Szabo, I. Babuska, in: Finite Element Analysis, Wiley, New York, 1991. [24] I. Babuska, B. Andersson, P.J. Smith, K. Levin, Damage analysis of ®ber composites. Part I: statistical analysis on ®ber scale, Technical Report, FFA TN 1998-15, Flygtekniska F ors oksanstalten (Aeronautical Research Institute of Sweden), 1998. [25] I. Babuska, T. Strouboulis, D.K. Datta, K. Copps, S.K. Gangaraj, A posteriori estimation and adaptive control of the error in the quantity of interest. Part I : a posteriori estimation of the error in the Mises stress and the stress-intensity factors, Comput. Methods Appl. Mech. Engrg. 181 (2000) 261±294. [26] G. Booch, in: Object Oriented Design with Applications, Benjamin/Cummings, New York, 1991. [27] S.J. Chapman, in: Fortran 90/95 for Scientists and Engineers, McGraw-Hill, New York, 1997. [28] J. Rumbaugh, I. Jacobson, G. Booch, in: The Uni®ed Modeling Language Reference Manual, Addison-Wesley, New York, 1998.