A numerical integration scheme for special finite elements for the Helmholtz equation

The theory for integrating the element matrices for rectangular, triangular and quadrilateral finite elements for the solution of the Helmholtz equation for very short waves is presented. A numerical integration scheme is developed. Samples of Maple and Fortran code for the evaluation of integration abscissæ and weights are made available. The results are compared with those obtained using large numbers of Gauss–Legendre integration points for a range of testing wave problems. The results demonstrate that the method gives correct results, which gives confidence in the procedures, and show that large savings in computation time can be achieved. Copyright © 2002 John Wiley & Sons, Ltd.


INTRODUCTION
This paper is concerned with recent developments in solving the Helmholtz equation for short wave problems. It is written here for reference as @ 2 @x 2 + @ 2 @y 2 + k 2 =0; ∇ 2 + k 2 =0 (1) k is the wavenumber, given by k =2 = , where is the wave length. is the ÿeld variable of interest, which could be wave elevation, pressure, or an electro-magnetic potential, among many other possibilities.
A weak form of the problem is obtained in the usual way by multiplying Equation (1) by a function , integrating that product over the domain . We assume the domain to be the ÿnite annulus bounded by the closed curves 1 and 2 which would normally be the scattering surface and the outer boundary of the ÿnite element model, respectively, though this is not strictly relevant to the discussion. We then apply the divergence theorem to produce (∇ · ∇ − k 2 ) d = 2 ∇ · n 2 d 2 − 1 ∇ · n 1 d 1 (3) where n i is the outward-directed unit normal vector on i . We obtain a Galerkin ÿnite-element method or more precisely, a Galerkin PUFEM, by taking Our interest is not in the boundary integrals on the right-hand side, but in the integral over the domain . On substituting in the expansions for and , we see that a typical element of the coe cient matrix has the form Next we insert the identity ∇ j =ikê j j into the above equation to obtain the form The hat functions are non-zero over only a few neighbouring ÿnite elements. We can further restrict our attention to the form of the integral over a typical ÿnite element E t , The integral of Equation (7) is highly oscillatory because of the presence of the exponential terms of imaginary argument, which are trigonometrical functions. Since k may be large, the exponential function in Equation (7) may contain many wavelengths. Previous workers, when evaluating the above integral have resorted to very high order Gauss-Legendre integration.
In the examples shown by Laghrouche and Bettess [7] up to 120 by 120 Gauss-Legendre numerical integration is used. So the reduction in the number of active variables, in the system matrix, comes at the price of some computationally intensive numerical integrations over the domain of the elements. The time for these integrations may be comparable with the execution time for the system matrix equation solution. The Gauss-Legendre method assumes that the function to be integrated is of a polynomial form. But it is well known that polynomial representations of trigonometrical functions are not accurate and are also expensive to compute [22]. Indeed it seems completely illogical to develop new types of ÿnite elements specially for short waves, which include a knowledge of wave behaviour, and then to ignore this knowledge, and revert to polynomials, at the integration stage. The challenge then is to develop better integration methods, which incorporate our knowledge of the plane wave functions. This should drastically reduce the time taken to form the element matrices. The following analysis demonstrates how the integrations may be carried out more e ciently, by developing semi-analytical integration rules.
At this point, we introduce the assumption that the local-to-global geometric mapping is a bilinear function of the element co-ordinates, which implies that the edges of the elements are all straight. ‡ For a quadrilateral element, the global cartesian co-ordinates x =(x 1 ; x 2 ) are related to the local element co-ordinates ( ; Á) on E t by where the coe cients a t , b t and c t are vectors of dimension 2. The Á term complicates the evaluation of the integrals. Full details for the case where d t = 0 are given in Reference [23].
Here, we are mainly concerned with triangular or rectangular ÿnite elements, for which d t =0.
In this case, we see that the expression in the square brackets in (7) above is a quadratic polynomial in ( ; Á). (This is because the hat functions are linear and the Jacobian of the geometric mapping is constant.) Consequently, the integral in (7) has the form A t is the area of ÿnite element E t , andÂ is the area ofÊ, the standard rectangular or triangular element. L ' ( ; Á) is an interpolating polynomial, which takes unit value at an integration point, and is zero at all the others. For rectangular and quadrilateral elements it can be constructed from the product of Lagrange polynomials (see Equation (20)). In those cases L ' ( ; Á)=L p ( ) × L q (Á). If the same mapping were retained, but higher-order interpolation were adopted, the upper limits on p and q would, of course, increase. Finally, we introduce and substitute these into (9) to obtain The analysis has now arrived at an integration rule of the following form: where w pq (ÿ; )= Ê e iÿ e i Á L ' ( ; Á) d dÁ (14) ‡ Higher-order mappings could be used, but the integrations would not be exact, which is the case in normal ÿnite elements.
For the rectangle, the evaluation of the integration weights is simpliÿed because the integral of products becomes a product of integrals, (15) Both integrals on the right-hand side are easily evaluated by repeated application of integration by parts. The evaluation of the integration weights for the standard triangle is more complicated. In this case The inner integral can be evaluated through integration by parts. The outer integral, while complicated, can also be evaluated analytically.

RECTANGULAR FINITE ELEMENTS
In this case the relation between the local, ; Á and global (x 1 ; x 2 ) co-ordinates given in Equation (8) simply becomes where a and b are the dimensions of the rectangular element in the x and y directions, respectively. The resulting integral is given above as Equation (15). As indicated above the two integrations in the x 1 and x 2 directions can now be carried out in sequence. Now consider the simplest case, when p=0. The integral in the direction is simply and integrals of terms containing higher-order polynomials can be found by integrating by parts. A di culty arises since it is possible for ÿ to take the value 0. In this case the integral given in Equation (18) above must be evaluated using l'Hôpital's rule. If ÿ is small then numerical evaluations based on the Equation (18) would also be inaccurate. In such cases it is therefore necessary to adopt a series form for the exponential. This can be integrated term by term. A total of four special cases arises as indicated in Table I.
When the absolute value of ÿ is below some set value, , the function e iÿ , is replaced by the truncated series 19) where N denotes the number of terms retained in the approximation. This can be adjusted, using the value , so that the truncation error in the series is of the order of machine accuracy. The integral of the series in Equation (19) is well behaved as ÿ → 0, and so a seamless transition between the integrated form for the exponential and the series should be achieved. A similar consideration applies in the case → 0, of course.

THEORY OF NUMERICAL INTEGRATION
There are many books and papers devoted to this topic, and we will make no attempt to cover the same ground. The necessary theory can be read in, for example, References [22,24] or [25]. The simple procedure adopted here can be summarized as follows: • Choose a set of integration absciss in two dimensions, distributed over the domain of the ÿnite element. In the present case, the integration absciss were chosen to be equally spaced in the and Á directions, and to include absciss at the ÿnite element vertices. For triangles the ranges were 06 61 − Á and 06Á61, and for the rectangle and quadrilateral −16 61 and −16Á61. Equal spacing was chosen because it was simple. § The number of points selected depends upon highest powers of and Á in the element. A bi-linear rectangular element matrix would have highest powers of 2 Á 2 , by deÿnition. This necessitates three integration points in each direction giving a total of nine points. In general if the highest power of and Á is n, then (n + 1) 2 points will be needed. A bi-quadratic element would need 25 integration points. • The Lagrange polynomials for each integration point are set up using the classical formula for such polynomials, [22,25], in both and Á directions.
The polynomials for each point in the and Á directions are found and multiplied together. ¶ The next step is to integrate each polynomial over the domain of the element. • Integrate analytically the product of the Lagrange polynomial and the wave functions over the domain of the ÿnite element, as seen in Equations (15) and (16) for the rectangle and triangle, respectively, to give the weights. § It may be possible to select optimal locations for the integration points, as in Gauss-Legendre integration. However, this would be almost impossibly di cult and the optimal locations would be di erent for every value of ÿ and . There is a considerable literature on optimal integration absciss which cannot be explored here. See, for example References [22,24,25]. ¶ This is for the rectangle and quadrilateral. The corresponding polynomials are also available for triangles [23].
All these operations are ideal for being performed by a computer algebra program. The integration by parts, in two variables, gives rise to a very large number of terms, which are di cult to evaluate by hand. But because they are routine and standardized they are ideal for computer algebra packages. The existence of special cases when ÿ → 0 and → 0 further complicates the manipulations. The actual manipulations were carried out using Maple and Mathematica.

A DESCRIPTION OF THE MAPLE CODE
A web page is available which gives all the Maple codes and fortran codes, used to generate the results in this paper, [26]. The and Á co-ordinates are easily set for each integration point A switch for waves or no waves enables the program to be used simply for Newton-Cotes integrations if wished. In the latter case the ÿ and are set to zero. This is a simple and useful check since results can be compared with tables.
In the Maple code the two Lagrange polynomials for each integration point, were set up using a Maple procedure lagpoly, which simply implemented the classical formula for such polynomials, expressions (20).
After forming the product of each pair of polynomials, it was multiplied by the appropriate harmonic function and integrated analytically between the limits −1 and +1 in the and Á directions, using Maple. This value is then the weight for that integration point. The expressions for the weights in terms of the variables, a, b, ÿ and , are output as Fortran expressions. Maple can automatically generate Fortran or C code to generate these weights. The Fortran option is used here, but the change to C is tri ing.
The Fortran output is organized to be in the form of a subroutine, which can return the weights, on being supplied with the wavenumber, wave directions and element dimensions. This can be compiled directly by standard Fortran compilers. This would be called for a range of values of ÿ and . These in turn would depend upon the values of wavenumber, k, and wave directions Â j .
It might be thought that an integration scheme which has to re-calculate the integration weights for every pair of wave directions would be ine cient. However, for short waves the expense of calculating and storing a large range of integration weights is small compared with the expense of the computations which must be carried out for all terms in the element matrix at each Gauss-Legendre point, in a numerically integrated element. The number of computations is also completely independent of the wavenumber.

INTEGRATION SCHEME LOGIC
Unlike conventional integration schemes, in which the same integration weights are used for all the terms evaluated at a given integration point, in this case each pair of wave directions requires a di erent, complex, weight. There is a signiÿcant overhead of storage and computation in evaluating and storing these integration weights. However, it is still economical compared with conventional Gauss-Legendre integration.
The procedure is as follows. For the ÿnite element all the directions at each node are assembled in a list. A simple example in which nodes 1, 2, 3 and 4 have 3, 2, 5 and 3 directions respectively, is shown below: where Â i j denotes the jth wave direction at the ith node. Let the total number of wave directions be M . In the above example there is a total of 13 wave directions. For each pair of directions, that is for M × (M + 1)=2 cases the integration weights are evaluated by a call to a fortran subroutine. In the example above there are 13 directions and hence 13 × 14=2=91 cases. The active arguments are the dimensions of the rectangular elements, a and b, and the exponents F and G, which are evaluated from Equation (9) from the wave number and wave directions. The integration weights are stored in a large array of dimension M × M × (n + 1) × (n + 1), where n is the highest power of and Á encountered. 5 is an appropriate value, if bi-quadratic elements are used. The weights are evaluated prior to the main integration loops.
Each integration loop is more complex than the conventional Gauss-Legendre-type integration, because as well as the main loops over all the integration points in the and Á directions, it is also necessary to have inner loops over all wave directions. However, the total number of multiplications of terms to evaluate the weighted residual is no larger, it is simply partitioned into smaller sections which use di erent weights. Moreover, some terms in the shape function and shape function derivatives are now present by implication, and not explicitly, so that the formation of the terms corresponding to a single integration point requires less e ort than for Gauss-Legendre integration. And, of course, there are, in general, far fewer integration points. In very general terms the logic is as follows: For all x integration points do For all y integration points do For all i wave directions do For all j wave directions do Form contribution to element matrix at this integration point for these directions Multiply by weight for this point and these directions end do end do end do end do 6. EXAMPLE COMPARISONS Consider a special wave element with 36 wave directions at each node and 9 nodes, using 100 × 100 Gauss Legendre integration. At each Gauss point three inner products of shape functions and shape function derivatives have to be formed to form the contribution to the element matrix. The size of the inner products is number of nodes multiplied by the number of wave directions, that is 36 × 9=324. Hence the number of oating point multiplications is 100 × 100 × 3 × 324 ≈ 10 000 000 In the case of 5 × 5 set of integration points, the total number of oating point multiplications is So provided the number of operations in evaluating the integration weights is less than 9 × 10 7 , the forming of the integration weights is worthwhile. That is we can a ord up to ÿve million oating point multiplications in evaluating the weights for each integration point. With 36 possible directions and 36 × 37=2 ≈ 670 combinations, even 500 oating point multiplications for each of these, is still worth while. In practice it is often possible to remain well below that ÿgure.

NUMERICAL RESULTS FOR RECTANGULAR ELEMENT
A program was written which would generate the element matrices using the special integration rules developed here, or many conventional Gauss-Legendre integration points. The element matrices obtained in this way were compared with each other. It was expected that the special integration rule would generate the element matrix accurately, irrespective of the wave number, k, whereas the Gauss-Legendre scheme would only be accurate if a su ciently large number of integration points were used. This proved to be the case.
The comparison between the semi-analytical integration and the Gauss-Legendre numerical integration was carried out as follows. The squares of the absolute values of the di erences were summed for every term in the two matrices and then divided by the sum of the squares of the absolute values of the terms in the matrix obtained by semi-analytical integration.
where K is the element matrix, n is the total number of degrees of freedom in the element, (an) denotes the semi-analytical integration, (gl) denotes the Gauss-Legendre integration, and e is the root mean square error. This error is plotted in Figure 1.
The element was square, with length of side 2h. The program was run for a range of kh values between 4 and 140, and for a range of conventional Gauss-Legendre integration points between 1 and 120. A value of kh=140 corresponded to a wavelength, , of approximately =2 h=140 ≈ 0:045h This represents roughly 2=0:045 ≈ 45 wavelengths in the element.
Since the semi-analytical integration should be independent of the wavelength, (apart from the issue of rounding errors at very high values of k), we can predict the behaviour of the results. We should expect the di erence between the semi-analytical integrated value for the element matrix, and the numerically integrated value to be small for small values of k. As k increases we should expect the di erence, or error, to increase, when the Gauss-Legendre scheme no longer has su cient integration points to resolve the wave details. This is exactly what happens. Figure 1 shows a plot of the error against the wavenumber, k, for a range of numbers of integration points. The semi-analytical and numerical errors increase at a value of k, which increases with the number of integration points. For very low numbers of integration points, between 1 and 10, the element matrix is not resolved at all with any acceptable accuracy. For numbers of points greater than 10, the error is around 10 −7 to 10 −6 , until the wavenumber increases to a value at which the number of integration points cannot resolve the wave details. Then as the wave number increases there is a rapid rise in the RMS error, which then plateaus at absurd and unacceptable values of about 10 3 . For example this transition takes place for 100 integration points, for k just below 100. For the results shown, four directions per node were used. Although the results are not shown, the program was also run for a range of di erent wave directions, and also for the special cases of F ≈ 0, G ≈ 0 and F ≈ G ≈ 0. (See Section 2.) No signiÿcant di erences in the results were observed.

TESTS AND TIMING RESULTS FOR RECTANGULAR ELEMENT
The two programs were used to evaluate the element matrices. Of course this is a subjective comparison, since one type of element may be coded more e ciently than the other. The speed-ups varied with the wave length, as shown in Figure 2. The time taken to integrate the element matrix semi-analytically has been found, using timing routines. This is denoted as T sa . The time to integrate the element matrix numerically has been similarly denoted T gl . The ratio of these two times has been termed, for brevity, 'Ratio of analytical to numerical times for element matrix integration', R=T sa =T gl . In Figure 2 this time is plotted as a function of the number of Gauss-Legendre integration points in and Á directions. The plot in Figure 2 demonstrates that the semi-analytical integration gives a considerable beneÿt, especially for shorter wave lengths. The break-even point is for 10 × 10 Gauss-Legendre integration, which corresponds, roughly, to one wavelength in each direction in the special ÿnite element. For 120 × 120 integration, or 12 wavelengths in each direction, the advantage of the semi-analytical integration is a speed-up factor of about 120. This is the sort of zone in which the advantages of the special wave ÿnite elements are particularly felt. It should also be noted that no special e orts have been made to 'tweak' the integration code produced by Maple to be optimal. In general, although the code produced by Maple is adequate, it is not particularly e cient. Variables are often deÿned repeatedly. Also internal symmetries in the calculation of the integration weights have not been exploited. For these, and other reasons, it is felt that the speed-ups shown, are the minimum which can be expected from such techniques. (Fortran code produced by Mathematica was much worse than that produced by Maple.)

TRIANGULAR FINITE ELEMENTS
The integral for the weights for integration over the triangular ÿnite element is given above as Equation (16). and Á can be identiÿed with the area co-ordinates L 2 and L 3 of the triangle, respectively, as is well known, and is shown in Figure 3.
Because of the limits of integration for the triangle, another special case arises corresponding to that mentioned after Equation (21). The case of ÿ ≈ must be considered. Terms of the form e i(ÿ− )Á must then be replaced by the equivalent series. This gives a total of ÿve special cases, listed in Table II. The polynomials for the triangle must be evaluated in terms of local area, or , Á coordinates, using expressions which are in standard textbooks [23], and which will not be repeated here. The integration points are chosen to be regularly spaced in and Á. The simplest formula thus has 1 integration point, at the centre of the triangle (although this would be of little practical interest), the second has three integration points, at the three vertices, the third has in addition the midside points, and so on. The layouts are sketched in Figure 3.

QUADRILATERAL FINITE ELEMENTS
In the case of the bi-linear quadrilateral the Jacobian of the mapping from local to global coordinates is no longer constant. The integrations will therefore only be approximate. However, it is possible to follow the same procedures as before. An additional term to those of Equations (11) is introduced, After some manipulation and the assumption of a constant Jacobian the expressions for the weights become e iÿ e i Á e i Á L ' d dÁ (24) In the integral of Equation (24), the di culties arise from the presence of the Á term. However, the integrals, while complicated, can be carried out. In the quadrilateral, there are three constants, ÿ, and , and a total of 12 special cases arise, corresponding to these values being zero or non-zero, and to ÿ or being approximately equal to , which again causes a problem. When is non-zero, the weights involve the exponential integral, Ei(x). This function is not supplied in Fortran, and so it was evaluated in our code in two separate ways. The ÿrst method was using the code from Numerical Recipes [27]. The second was by using code written by Bulirsch [28] and modiÿed by Clark [29]. This code uses Chebychev coe cients to evaluate the Cosine and Sine integrals. Polynomials are evaluated using Clenshaw's algorithm. The values for Ei(x) obtained from the two codes were checked against each other and against values evaluated to arbitrary precision by Maple. Because 12 separate special cases arise, details are given elsewhere [30].

Evanescent modes
In order to test the rectangular element integration schemes on a realistic and non-trivial problem the so-called evanescent mode problem was solved in the interior of a rectangular domain. While the comparisons given in Figure 1 are useful, because of the ill-conditioned nature of the matrices generated by these elements it is essential to test the element matrices on a real problem. The solution to this problem is given by Morse and Feshback [31], who demonstrate that the Helmholtz equation can be separated into two equations, so that Equation (1) becomes 1 X where k is the wavenumber, given by k =2 = , where is the wavelength, and =X (x)Y (y). After some manipulation this leads to a solution of the form where is an arbitrary constant. For values of smaller than k, the conventional plane wave propagating in a direction deÿned by the two x and y components is obtained. But if ¿k, the behaviour in the y direction is exponential, and the apparent wave length in the x direction is given by 2 = , and is therefore less than the true wavelength, given by 2 =k. Since the solution in this case has a di erent wavelength from that assumed in the special element shape functions it is a searching test of the approach. The solution can be further generalized by rotating it through an angle ÿ. The most general form is therefore where in the usual notation, with ÿ being the angle of rotation of the wave ÿeld s=x cos ÿ + y sin ÿ; t = − x sin ÿ + y cos ÿ The boundary condition applied on all four boundaries is the Robin boundary condition [7] @ @n + ik =g (29) where g is the known solution, derived from Equations (28) and (29) above, and n is the outward normal. For the numerical solutions, the domain is 06x612, 06y61:2. The problem is ÿrst solved for the parameters k =1, =1:5 and ÿ =30 • . For this case a mesh of 4 × 4 equal sized rectangular elements were used, with 8 wave directions at each node. Figure 4 shows the real part of the analytical solution for the evanescent mode. The real and imaginary potential along the centreline of the domain, y =0:6 is plotted in Figure 5, for analytical, and the numerically integrated and the semi-analytically integrated rectangular elements. Figure 6 shows the errors in the ÿnite element solutions along the same line.
The same evanescent modes problem is next solved for the parameters k =3, =5 and ÿ =10 • . The same element mesh was used as before but with 24 wave directions at each node. Figure 7 shows the real part of the analytical solution for the evanescent mode. For the numerical solutions Robin boundary conditions are used. The real and imaginary potential along the centreline of the domain, y =0:6 is plotted in Figure 8, for analytical, and the numerically integrated and the semi-analytically integrated rectangular elements. Figure 9 shows the errors in the ÿnite element solutions along the same line. The results are encouraging. Aside from the generally acceptable accuracy of the results, they demonstrate that the results from the numerically integrated elements and the semi-analytically integrated elements are virtually identical. This gives conÿdence in the scheme for the rectangular element.

Circular cylinder di raction
The problem considered consists of a rigid circular cylinder of radius, a=1. The ÿnite element mesh extends from r=a=1 to 5. The mesh is regular, with 36 sectors, each subtending 10 • . Each sector is modelled using 8 linear triangular elements. (Linear in the sense that the  underlying shape function is linear.) The wavenumber, k is chosen so that ka=4. The mesh is shown in Figure 10. Figure 11 shows the results for plane waves di racted by the cylinder, for ka=4, where k is the wave number and a is the cylinder radius. All the waves considered are incident from the negative x direction, that is from Â =180 • and r → ∞. They progress towards Â =0 • and r → ∞. All incident waves have unit amplitude. The numerically integrated elements use a Gauss-Legendre scheme involving 8 integration points in two directions. (This unsymmetrical scheme is not entirely satisfactory, but we had no access to higher triangular  integration schemes.) The elements on the boundary, used to implement the boundary conditions were integrated using 40 Gauss-Legendre integration points, and so the order of quadrature should not in uence the results, as it is so much higher than the minimum requirement. Plane damper boundary conditions were used, since they are simple to implement, and the present investigation is to determine the accuracy of the elements. The authors are aware that much more accurate techniques are available. Both numerically integrated and semianalytically integrated elements use 8 equally spaced wave directions, i.e. 0 • , 45 • ; : : : ; 315 • . Figure 12 shows the errors in the real and imaginary components of the complex wave potential as a function of the angle around the cylinder. The plots are of the di erence between the  analytical expression and the numerical values. The results demonstrate that there are virtually no di erences between the values obtained by the analytically integrated and numerically integrated elements. Figure 13 shows the results for plane waves di racted by the cylinder, for ka=16, where k is the wave number and a is the cylinder radius. The numerically integrated elements use a Gauss-Legendre scheme involving 24 integration points in two directions. The same arrangement as before served for boundary elements and plane damper boundary conditions were again used. Both numerically integrated and semi-analytically integrated elements use 18 equally spaced wave directions, for this shorter wave problem i.e. 0 • ; 20 • : : : 340 • . Figure 14 shows the errors in the real and imaginary components of the complex wave   Table III. The computations were carried out on a standard PC. In the case of the numerically integrated 24 × 24 integration point elements, the time to form the element matrices, was 899 s. That for the semi-analytically integrated elements was 60 s. This latter time was about the same as the system matrix solution time. Too much should not be read into one set of ÿgures, but they do indicate that the time for the formation of the element matrices can be brought broadly into line with the time for the solution of the system matrix, which indicates that the method should be competitive. The ÿgure for the semi-analytical integrations should be independent of wave length. If the criterion of 10 nodes per wavelength, for wave modelling, is used a regular mesh would require 101 859 nodes, and an unstructured mesh 48 892 nodes. This assumes that pollution would not play a role, which might well be optimistic in a problem with over 12 wavelengths in the radial direction. In the present model, with 180 nodes, and 18 directions per node, the total number of degrees of freedom is 3240. This gives a factor of reduction of number of degrees of freedom for the unstructured mesh of roughly 15, and for the structured mesh roughly 31. It is likely that the number of nodes in the present mesh could be reduced if an iso-parametric version of these elements (which could model more accurately the curved cylinder), were developed, since the nodal spacing circumferentially is much closer than in the radial direction.

CONCLUSIONS
A method for integrating semi-analytically the element matrices arising from the new schemes of special ÿnite elements for the solution of Helmholtz equation in two dimensions has been presented. The results show that the integration schemes give results which are hardly distinguishable from those obtained using very large numbers of conventional Gauss-Legendre integration points. As expected the element matrix integrations appear to be independent of wave number, as shown by testing. For larger wavenumbers the semi-analytical integrations lead to large savings in execution times, by factors which increase with increasing wavenumber. Preliminary results for some classical wave problems demonstrate that these savings seem to be achieved without any loss of accuracy. At any rate for the elements of simpler geometry, semi-analytical integration seems to be a very promising way forward. There remain the challenges of higher-order parametric elements and the extension to three dimensions. In three dimensions the beneÿts of this approach should be even more marked, if it can be got to work. As the semi-analytical integration codes are still experimental it is likely that they could be improved upon. The other approach, that of Ortiz [16], could well be even more e cient. Time and experience will tell which method is the most e ective, but it would be premature to make a judgement at the moment.
The PostScript ÿle containing this paper, and the Maple codes to determine the integration weights, the corresponding generated Fortran integration weight codes and the test programs which generated the element matrix comparisons shown in Section 9, are all freely available on the web site [26].