SHORT COMMUNICATIONS 405
It usually makes little difference whether moments are extrapolated from values at the Gauss
points or calculated directly at the corners; this result is in agreement with flat plate test cases.
In the four-elementmesh force NZBis considerablyimproved if extrapolation is avoided. However,
the opposite trend is seen in certain of the one-element values of NxB.This effect may be explained
by considering the variation of midsurface strain E,, which is directly proportional to N, because
Poisson’s ratio u = 0. In each element, E, varies quadratically with co-ordinate 4. Extrapolation
imposes a linear variation of E, with 4. Apparently, the beneficial effect of calculating more
accurate strain values at the Gauss points may or may not be overcome by the detrimental effect
of discarding the quadratic variation.
As was the case with flat plates,2there seems no reason to use five internal freedoms per element
rather than three. It is concluded that because internal freedoms are more often helpful than
harmful in plate and shell problems, computer programs should permit the user to employ
three internal freedoms in the eight-node element. As a general rule it seems that all strains at
element corners should be calculated by extrapolation from the Gauss points rather than directly.
REFERENCES
1. 0. C. Zienkiewicz, R. L. Taylor and J. M. Too, ‘Reduced integration technique in general analysis of plates
and shells’, Znt. J. num. Meth. Engng, 3,275-290 (1971).
2. R. D. Cook, ‘More on reduced integration and isoparametric elements’, Znt. J. num. Meth. Engng, 5, 141-142
(1972).
3. S. F. Pawsey and R. W. Clough, ‘Improved numerical integration of thick shell iinite elements’, Znt. J. num.
Meth. Engng, 3, 575-586 (1971).
4. 0. C. Zienkiewicz, The Finite Element Method in Engineering Science, McGraw-Hill, London, 1971.
5. G. R. Cowper, G. M. Lindberg and M. D. Olson, ‘A shallow shell element of triangular shape’, Znt. J. Solids
Strucr. 6, 1133-1156 (1970).
6. D. G. Ashwell and A. B. Sabir, ‘A new cylindrical shell finite element based on simple independent strain
functions’, Znt. J. mech. Sci, 14, 171-183 (1972).
7. F. K. Bogner, R. L. Fox and L. A. Schmit, ‘A cylindrical shell discrete element’, AZAA Jnl, 5, 745-751
(1967).
8. G. Cantin and R. W. Clough, ‘A curved, cylindrical-shell, finite element’, AZAA Jnl, 6, 1057-1062 (1968).
9. A. C. Scordelis and K. S. Lo, ‘Computer analysis of cylindrical shells’, ACZ Jnl, 61, 539-561 (1964).
GAUSSIAN QUADRATURE FORMULAS FOR TRIANGLES
G. R. COWPER
National Research Council of Canada, Ottawa, Ontario, Canada
SUMMARY
Several formulas are presented for the numerical integration of a function over a triangular area.
The formulas are of the Gaussian type and are fully symmetric with respect to the three vertices of
the triangle.
INTRODUCTION
In finite element work, the need sometimes arises to numerically integrate a function over a
triangular area. Suitable formulas have been derived by Silvester; Irons2 and Hammer and
co-worker~.~ The formulas of the latter were compiled in a convenient form by Felippa: and the
compilation has been reproduced, along with Iron’s formulas, in the book by Zienkiewic~.~
The formulas of Silvester have simple coefficients and are fully symmetric with respect to the
three vertices of the triangle, but they are of the Newton-Cotes type and thus are relatively
Received 21 April 1972
406 SHORT COMMUNICATIONS
inefficient compared with the Gaussian formulas. The formulas of Irons, which make use of
successive application of the Gauss and Radau quadrature rules, are highly efficient but have
the aesthetically unappealing feature that the sampling points are not arranged symmetrically
within the triangle. From the point of view of both numerical efficiency and aesthetics the best
formulas are those of Hammer and co-workers which are of the Gaussian type and fully
symmetric with respect to the three vertices of the triangle.
Five different formulas were derived by Hammer and co-workers which have truncation
errors ranging from second- to sixth-order. Unlike the one-dimensional case, the theory of two-
(and three-) dimensional integration is not highly developed. The discovery of additional formulas
remains a matter of ad hoe investigation in each specific case. The purpose of this note is to
present additional formulas of the Hammer and co-workers type which are believed to be new
and which were derived in the course of research on triangular finite elements for shells. The
formulas will be displayed first, followed by some remarks on their derivation.
FORMULAS
The formulas for numerical integration of a functionfover a triangle of area A are all of the form
SI fdA = A i=l
N
x wi f(&,74,w
where &,T,, {, are the area co-ordinates of the i-th sampling point and w, is the weight associated
(1)
with the i-th sampling point. Values of the constants w,,&,yi, Ci for the various formulas are
listed in Table I. For completeness, three of the formulas due to Hammer and co-workers are
also included in the table.* Readers unfamiliar with area co-ordinates may refer to Zienkiewicz’s
book?
Because of the triangular symmetry of the formulas the sampling points occur in groups of
either six, three or one. Thus, if a sampling point has area co-ordinates (a,p, y), none of which
are equal, then there are also five symmetrically disposed points with co-ordinates (a,y , p),
(p, a,y), (p, y, a), (y, a,p) and (y, p, a), all of which have the same weight. A sampling point which
has two equal co-ordinates, for example (a,/?,/?), occurs as a member of a trio of sampling
points with equal weights, the co-ordinates of the other two points of the trio being (p, a,p) and
(p, p, a). A single sampling point with co-ordinates (6, 6, &)may occur at the centroid of the
triangle. To save space the weight and co-ordinates of only one point of a group are given in
Table I. The entry in the column headed ‘Multiplicity’ indicates whether the point belongs to a
group of one, three or six points. The ‘Degree of precision’ shown in the table indicates the
highest degree polynomial which the formula integrates exactly.
REMARKS ON DERIVATION
The formulas are designed to exactly integrate complete polynomials of a given degree. Since
the quadrature-formulas are to be exact for all polynomials of a complete set, one can write (1)
for each such polynomial and the relations so obtained are the equations which determine the
parameters w,,&,q,, Q. The degree of completeness and number of sampling points must of
course be selected so that the number of unknown parameters is consistent with the number of
available equations.
The equations for w,, ti, qi and [,
are highly non-linear and their solution is not straight-
forward, although some steps can be taken to put the equations in a rather manageable form.
* This also affords the opportunity to correct the four-point formula, which is given incorrectly both by
Felippa‘ and Zienkiewic~.~ The ‘severe cancellation errors’ which are said to occur with this formula may
be due simply to use of an incorrect formula.
SHORT COMMUNICATIONS 407
It was found convenient to work in Cartesian co-ordinates originating at the centroid rather
than area co-ordinates, and since the quadrature formulas are independent of the shape of the
triangle it was permissible to use the most convenient shape, namely a unit equilateral triangle.
The requirement of triangular symmetry means that the sampling points can occur only in
symmetric groups of one, three or six points, and imposition of this requirement reduces the
number of unknown co-ordinates to a minimum. The number of independent equations was
found to be correspondingly reduced. The three polynomials of second degree furnished only
one essentially independent equation, as did the polynomials of third, fourth, fifth and seventh
degrees. The seven polynomials of the sixth degree gave rise to two independent equations, while
the first degree equations were identically satisfied as a result of the imposed symmetry. In some
cases it was possible, by judicious manipulation, to solve the equations in closed form. In others,
the solution had to be obtained by trial and error, after first eliminating all but one of the un-
knowns. In all cases the final formulas were checked by using them to integrate a complete set
of polynomials.
Although the derivation was carried out for an equilateral triangle, any triangle can be mapped
onto an equilateral triangle by a linear transformation. Area co-ordinates are invariant under
such mapping, and hence the quadrature formula (1) applies to triangles of any shape.
Table I
Wf 5j 779 5i Mu1tiplicity
3-point formula Degree of precision 2
0.33333 33333 33333 0.66666 66666 66667 0-1666666666 66667 0.16666 66666 66667 3
3-point formula Degree of precision 2
0.33333 33333 33333 0.5oooO 0000000000 0-5oooO 0000000000 0.000000000000000 3
4-point formula Degree of precision 3
- 0.56250 0000000000 0.33333 33333 33333 0.33333 33333 33333 0.33333 33333 33333 1
052083 33333 33333 0.6ooOO 0000000000 0.2oooO 0000000000 0.2oooo 0000000000 3
6-point formula Degree of precision 3
0.16666 66666 66667 0.65902 76223 74092 0.23193 33685 53031 0.10903 90090 72877 6
6-point formula Degree of precision 4
0-10995 17436 55322 0.81684 75729 80459 0.09157 62135 09771 0-09157 62135 09771 3
0.22338 15896 7801 1 0.10810 30181 68070 0.44594 84909 15965 0.44594 84909 15965 3
7-point formula Degree of precision 4
0.37500 0000000000 0.33333 33333 33333 0.33333 33333 33333 0.33333 33333 33333 1
0.10416 66666 66667 0.73671 24989 68435 0.23793 23664 72434 0.02535 51345 59132 6
7-point formula Degree of precision 5
0.22500 0000000000 0.33333 33333 33333 0.33333 33333 33333 0.33333 33333 33333 1
0.12593 91805 44827 0.79742 69853 53087 0.10128 65073 23456 0.10128 65073 23456 3
0.13239 41527 88506 0.47014 20641 05115 0.47014 20641 05115 0-05971 58717 89770 3
9-point formula Degree of precision 5
0.20595 05047 60887 0-1249495032 33232 0.43752 52483 83384 0-4375252483 83384 3
0.06369 14142 86223 0-79711 26518 60071 0.16540 99273 89841 0.03747 74207 50088 6
12-point formula Degree of precision 6
0.05084 49063 70207 0.87382 19710 16996 0.06308 90144 91502 0-06308 90144 91502 3
0.1 1678 62757 26379 0.50142 65096 58179 0.24928 67451 70910 024928 67451 70911 3
0.08285 10756 18374 0.63650 24991 21399 0.31035 24510 33785 0.05314 50498 44816 6
13-point formula Degree of precision 7
-0.14957 00444 67670 0.33333 33333 33333 0.33333 33333 33333 0.33333 33333 33333 1
0.17561 52574 33204 0-4793080678 41923 0.26034 59660 79038 0.26034 59660 79038 3
0.05334 72356 08839 0.86973 97941 95568 0.06513 01029 02216 0.06513 01029 02216 3
0.07711 37608 90257 0.63844 41885 69809 0.31286 54960 04875 0.04869 03 154 253 16 6
408 SHORT COMMUNICATIONS
ACKNOWLEDGEMENT
The assistance of Miss H. A. Tulloch in the preparation of this note is gratefully acknowledged.
REFERENCES
1. P. Silvester, ‘Newton-Cotes quadrature formulae for N-dimensional simplexes’, Proc. 2nd Can. Congr. Appl.
Mech., Waterloo, Canada (1969).
2. B. M. Irons, ‘Engineering applications of numerical integration in stiffness methods’, AIAA JI, 4, 2035-
2037 (1966).
3. P. C. Hammer, 0. J. Marlowe and A. H. Stroud, ‘Numerical integration over simplexes and cones’, Math
Tabl. natn. Res. Coun. Wash. 10, 13G137 (1956).
4. C. A. Felippa, ‘Refined finite element analysis of linear and non-linear two-dimensional structure’, Report
66-22, Dept. Civ. Engng, Univ. of California, Berkeley, Calif. (1966).
5. 0. C. Zienkiewicz, The Finite Element Method in Engineering Science, 2nd edn., McGraw-Hill, London, 1971.
ANALYSIS OF PLATES ON ELASTIC FOUNDATIONS
W . C. CARPENTER
Department of Engineering Science, University of Durham, Co. Durham, England
SUMMARY
The modification required to provide plate on elastic foundation analysis capability with existing
finite element algorithms is presented. As shown by an example, the modification is applicable to
configurations with varying ratios of plate to foundation stiffness.
INTRODUCTION
The finite element method may be adapted to investigate plates on elastic foundations. For a
plate resting on an elastic foundation, an element stiffness matrix is developed and a n example
is presented.
ELEMENT DERIVATION
For small deflections, neglecting the effects of shear distortion, the elastic strain energy of an
element with no foundation, U, is given as4
where
w = The deflection of the element = [PI [a]
[PI= A matrix of algebraic expressions in terms of the co-ordinates of the element
[a]= A matrix of undetermined coefficients
[SP= A matrix of nodal deformation parameters = [C][a]
F]= The element stiffness matrix
[E] = The [Q][a]
generalized strain matrix =
[a]= The generalized stress matrix = [D][E]
Received I March 1973
Revised 26 April 1973
