- Email: [email protected]
Gauss-Jacobi quadrature rules for n-simplex with applications to finite element methods
COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING NORTH-HOLLAND
GAUSSJACOBI
QUADRATURE
~PLI~A~ONS
TO ~
40 (1983) 293-307
RULES FOR n-SIMPLEX WITH FLUENT GOODS
Masayuki OKABE Mi~~~i Mining & Smelting Company Ltd., Niho~~~h~-MMro~~hi,
Chuo-Ku, Tokyo 103, Japan
Received 21 January 1983 Revised manuscript received 8 April 1983 For numerical integrations over the n-simplex, the Gauss-Jacobi quadrature operators are presented definitively. It is demonstrated in finite-element applications that the Gauss-Jacobi quadratures could be highly effective under the i~parametri~ mapping as well as under the semi-radial singul~ity mapping.
1. Introduction
Numerical integration plays an increasingly important role in numerical techniques such as the finite element and boundary element methods. A variety of finite element interpolations are known over the n-cube [l-S] as well as over the n-simplex [6-8], but we have some difficulty in approximate integrations to keep the required accuracy with necessary minimum numerical efforts. For the n-cube the product-Gauss rules [3,9] are fully explored. Here the operator is constructed as a product of Gauss-Legendre approximate integral operators defined on every axis. It is, indeed, optimal in the case of distorted cubes, although some redundancy exists for undistorted cubes. For the n-simplex, on the other hand, we have no complete devices. Only in two dimensions the so-called Gauss-Radau rules on triangles are fully known also in a product form [lo]. Gauss-Radau sampling points are located asymmetrically with respect to triangular vertices, and the same order accuracy could be attained with less sampling points. In fact, Cowper [ll] and Moan [12] simultaneously developed more effective formulas with symmetrical points but only in lower order accuracy. This paper concentrates on applications of Jacobi polynomials [13] to numerical integrations over the n-simplex. We first treat the isoparametric mapping, and Gauss-Jacobi quadrature operators are developed as natural generalization of Gauss-Radau ones. Then we define the semi-radial singularity mapping [8,14-16] in II dimensions, and the other kind of Gauss-Jacobi quadrature operators are presented.
2. Definite integral operators in one dimension
For the l-simplex operator I* by ~5-7825/83/$3.~
of (-1,l)
in the &-system, we define the analytical definite integral
@ 1983, Elsevier Science Publishers B.V. bosh-Holland)
M. Okabe, Gauss-Jucobi
294
quadrature for n-simplex
(1) Here @ is an appropriate real number (p > -1) and f denotes a sufficiently smooth function. Let Jc denote the approximate definite integral operator by
(2) Here tf denotes the 5, coordinate of a sampling point i, ~~~~) is the f-value at i and Wi designates the corresponding weight (i = 1, . . _, m). We then aim to integrate any polynomials up to 2m - 1 such that Jc [([$I = Ip [(&y] , PR~P~S~~O~
(3)
j = 0, . . . ,2nz - 1 .
2.1. Suppose that the system
m-i C ak/(P + k + I) = -2m +pCm/(/3+ m + I) ,
1 = 1,
. .
. >m
(4)
k=O
has the unique sohitions ak, and suppose that a polynomial by PY!%SI)= Z,+,Ct&
m-1 + ‘$I)>”f z: ak&l
+
(5)
6))”
k=O
has m-independent zero points. Here ~~+pl.$,= r(2m + /3 + l)/T(m + l)S(m + p + 1) with the gamma function r. Then (2) and (3) can be solved uniquely with respect to the weights Wi by selecting those zero points as [f. PROOF. Noting
(3) can explicitly be written as
2 {f(1+ [;))‘w,
= 2/(p+j+1),
j=O ,...,
2m-1.
(7)
i=l
Above-mentioned of the form
procedures
then guarantee
~{~(I+g;))“‘W=2/(p+j+I+2), i=l
with respect to Wi.
the identity of m systems with I = 0 to
j=O,...,m-1
m - 1
(81
295
M. Okabe, Gauss-Jacobi quadrature for n-simplex
3. Gauss-Jacobi quadrature operator in one dimension The polynomial Py(Q1) of (5) is known as the Jacobi polynomial in m degrees, see Appendix’ A. Let the monomials ((&)m)~=O orthogonalize :hrough the weight {i(l + &))B in (-1, l), then we have the Jacobi polynomials {P~(&))~=o. Here the orthogonal condition is written as (9) with nonnegative THEOREM satisfy (3).
delta S,,.
integers m and 12,and with the Kronecker
quadrature operator JB, can be uniquely defined so as to
3.1. The Gauss-Jacobi
PROOF. It is easy to prove that the Jacobi polynomial points in (-1, l), see Appendix B.
P?(&)
has m-independent
zero
We show the Gauss-Jacobi quadrature rules explicitly for some integer p in Tables 1 to 4. Zero points of PtB(&) and P?@(&) are plotted in Fig. 1 against @. REMARK
3.2. To obtain the formulas in closed form, we first make the Jacobi polynomial
-1
13* 0
5 1
2
3
4
5
6
Fig. 1. Zero points of Py8(&) and P$@(&) against j3.
M. Okabe, Gauss-Jacobi
296
Table 1 Gauss-Jacobi quadrature operators Ji which are identical to Gauss-Legendre ones. Computations of each abscissa and weight are done in double precision
m
Abscissa
1
0.0
2
3
4
5
51
Weight
W,
m
-0.57735027 0.57735027
1.0 1.0
1
0.33333333
-0.77459667 0.0 0.77459667
0.55555556 0.88888889 0.55555556
2
-0.28989795 0.68989795
0.36391724 0.63608276
-0.86113631 -0.33998104 0.33998104 0.86113631
0.34785485 0.65214515 0.65214515 0.34785484
-0.57531892 0.18106627 0.82282408
0.13965396 0.45848221 0.40186383
-0.90617985 -0.53846931 0.0 0.53846931 0.90617985
0.23692689 0.47862867 0.56888889 0.47862867 0.23692689
-0.72048027 -0.16718086 0.44631397 0.88579161
0.06236194 0.25969510 0.40692914 0.27101383
-0.80292983 -0.39092855 0.12405038 0.60397316 0.92038029
0.03149583 0.14781774 0.29277397 0.33434928 0.19356318
m
Abscissa
1
0.5
3
4
5
Table 2 Gauss-Jacobi quadrature operators Jk which are employed in Gauss-Radau rules. These operators can be applied also to the twodimensional semi-radial singularity mapping elements in linear fracture mechanics of A = l/2 as shown later
2.0
Table 3 Gauss-Jacobi quadrature operators J,?, which are used with Jz and J,!, in numerical integrations on tetrahedrons
2
quadrature for n-simplex
6;
Weight
Wi
3
4
5
Abscissa
5;
Weight
W,
1.0
Table 4 Gauss-Jacobi quadrature operators 54 which can be applied to the three-dimensional semiradial singularity mapping elements in the case of h = l/2 as shown later
m
Abscissa
0.66666667
1
0.6
0.5
-0.08830369 0.75497035
0.20157176 0.46509490
2
0.05971587 0.79742699
0.13381050 0.36618950
-0.41000442 0.30599247 0.85401195
0.05990141 0.29249254 0.31427272
-0.27347074 0.39762254 0.87584820
0.03295812 0.20919980 0.25784209
-0.59170284 -0.03409459 0.52279852 0.90299890
0.02070448 0.13726777 0.28691758 0.22177683
-0.47704442 0.07169289 0.58056646 0.91569416
0.00931673 0.08508345 0.21800874 0.18759108
-0.70210843 -0.26866695 0.22022723 0.65303936 0.93084212
0.00822765 0.06411120 0.17840032 0.25239792 0.16352957
-0.60757599 -0.16579958 0.29714001 0.69121030 0.93887141
0.00304138 0.03391465 0.12088991 0.20063301 0.14152106
3
4
5
51
Weight
Wi
M. Okabe, Gauss-Jacobi
explicit through the recurrence
297
quadrature for n-simplex
rule of the form
2(m + l)(m + P + 1)(2m + P)E?+&)
-(2m +p + 1Wm
+ P@m
+ P +2)&- ~"E!(&)
+ 2m (m + P)(2m + p + 2)PT_,(&) = 0 ,
m> 0.
(10)
Here E%)
= I,
W(51)
= w
(11)
+ 2) 51- PI *
Then we obtain its zero points by the Newton iteration referring to (I - &)DP?(&)
= (m + P + l)U?!(&)
(12)
- P??“(&))
where D = d/d&. REMARK
3.3. For the polynomial
f(51) = 2
f(!f~)h(tl)
+
f(tl) up to 2m - 1 degrees, we thus have
PZ(51)
mil
where cj denotes appropriate Wi =
coefficients. The weights Wi can then be obtained as
P[L*(&)] , i = 1,.
Here each interpolation Li(ei) =
(13)
Cj(5J
j=O
i=l
. . , m.
(14
function Li(t,) should be the polynomial
up to 2m - 1 satisfying
Sij a
(1%
It is especially remarked that the above-mentioned conditions allow unlimited selection of Li(tl)* If we adopt the Lagrangian interpolation polynomials with respect to 5: (i = 1, . . . , m) as usual, then naturally (15) is equivalent to (8) with 1= 0. REMARK 3.4. The Jacobi polynomial Pz(&) is identical Hence JZ is called the Gauss-Legendre quadrature operator.
to the Legendre
polynomial.
4. Isoparametric mapping between n-simplex and n-cube with collapse Consider an n-simplex in the global Cartesian system (x1, . . . , x,). Then the volume given by V, = det[“$ is1
“2 xj-l*q]/n! j=1
.
V,
is
(16)
298
M. Okabe, Gauss-Jacobi
quadrature for n-simplex
Here xi denotes the Xj coordinate of vertex i with xb = 1, and ei is the (n + 1)st unit vector having only the ith component of unity. The volume coordinate @k related to appex k can be defined as ~k=det[“~~(SkiX,-,+(l_4i)X;,)eie,]/n!V., i=l ;=I
k=l,...,n+l
(17)
where x0 = 1. Next we consider an n-cube normalized to (-1,l) in the local parametric system (tl, . . . , 5”). The subcube in II - 1 dimensions of [I = -1 is collapsed into node 1, and the subcube in n - 2 dimensions of (I = 1 and & = - 1 into node 2. In the same way, we collapse the subcube in it - 1 dimensions of [j = 1 (j < 1) and 4 = - 1 into node 1 up to I = n - 1. The residual two nodes with & = -1 and & = 1 are numbered II and II + 1, respectively. The collapsing procedures in two and three dimensions are, for clarity, illustrated in Figs. 2 and 3, respectively. Then the volume coordinates ok of (17) can be rewritten in terms of ,$I to & such that w1=
1(1- 61))
,i,k = { ‘fi
+(l - &) ,
$(I+
j=l
wn+1
=
fi
$(I
+
k = 2,
. . . , II,
(18)
5;.) -
i=l
THEOREM
4.1. The n-simplex
1
in (x1, . . . , x,)
to the (-1, 1) n-cube in
2
Fig. 2. Normalized 2-cube where node 1D is collapsed to node 1.
can be normalized
Fig. 3. Normalized 3-cube with collapse of nodes ID, ID’ and 1D” into node 1, and node 2D into node 2.
299
M. Okabe, Gauss-Jacobiquadrature for n-simplex
(51,* - * 5”) through the parametric representation 9
n+l
Xi=
J$X$Oky j=l,**.,
(19)
?Z.
k=l
Here
wk
is given by (18).
THEOREM
4.2. The Jacobian inherent to (19) with (18) can be writtenas
PROOF. It is clear in (18) that the summation
of all volume coordinates
takes unity;
(21) Then (19) can be expressed as n+l
xi-x$=
C (x!-xf)~~,
n.
j= l,...,
(2%
k=2
Differentiation
of (22) with respect to & to & thus gives (20).
REMARK 4.3. In usual finite-element applications, polynomial trial functions in terms of & to & are used. To the linear element, for example, volume coordinates of (18) are directly applied. Then the mapping function space and the trial function space have the same interpolation basis (@k), and hence the parametric representation of (19) with (18) could appropriately be called the isoparametric mapping.
5. Gauss-Jacobi quadrature rules for n-simplex under the isoparametric mapping Over the n-simplex we define the analytical integral operator I by (23) and the approximate
integral operator J,,, by
Jmlf(W,.*.,O n+1)1=v,
5 f(d,
i=l
. . . , d+,)W .
(24)
Here the number of sampling points is denoted by AX We aim to integrate ~lynomi~s up to 2m - 1 in terms of x1 to x, exactly by J,,,. Then
300
M. Okabe, Gauss-Jacobi
quadrature for n-simplex
noting (18) and (19) J,,, should realize
Jrn[fi(tj)k]=I[fi((j)k]T j=l
(25)
k=0,...,2m-1.
j=l
We restrict here f to monomials up to 2m - 1 in terms of volume coordinates. seperated in variables such that f(W1, .
* * 3 %+I)
=
fi
f;.(tj>
Then f can be
-
(26)
j=l
THEOREM 5.1. For f of (26), let the Gauss-Jacobi quadrature operator J,,, be constructed as a product of one-dimensional operators JE (p = 0, . . . , n - 1) by n-1
Jm [f
(w, . . . , ontl>l = (n! K/2”)
n J!Z[f n-dSn-d
.
(27)
p=o
Then J,,, realizes (25). REMARK
5.2.
In the Gauss-Jacobi quadrature rules, coordinates of sampling points are given in a tensor product of zero points of P?([,+) through p from 0 to n - 1. The weights are also given in a product of those for J!. Thus the number M of sampling points to realize (25) can be expressed as M=m”.
(28)
REMARK
5.3.
REMARK
5.4. The Gauss-Jacobi
REMARK
5.5. For
In order to construct the quadrature rules in the form of (24) sampling points should be written in volume coordinates by the use of (18). quadrature [lo], and hence called the Gauss-Radau ones.
rules for 2-simplex are established
by Radau
f of (26), the analytical integral operator I by (23) can also be separated in
variables such that
I[f(@l, . . . , @,+I)]= (n !VJ2”),@I”-‘[f;.Gfj>l . Here I”-’ is given by (1). Notice that fi in (29) could be nonpolynomials. (19) we can derive the formula of the form
(29) Then noting (18) and
(30) Here lj denotes the real number satisfying lj > -1. For nonnegative integers lj, (30) is well known especially in two and three dimensions
[3].
M. Okabe, Gauss-Jacobi
301
quadrature for n-simplex
6. Semi-radial singularity mapping between n-simplex and n-cube with collapse With the collapsing procedure alternately expressed as 01
=
1 - {$(1+ &)}lIA)
ok
=
I.%1
+
51))~‘~
previously argued, let the volume coordinates
02
=
G(1+
SW
tj)}!i(l -
(l@l%l+
[k)
O,Cl = {4(l+ 51)}l” { fi $(I +
[j)}
4(1-
7
52)
ok of (17) be
,
k = 3,. . . , n,
(31)
-
j=2
Here A is a positive constant below unity (0 < h < 1). THEOREM
6.1. The n-simplex in (x1, . . . , x,)
(51, * - * , &) through THEOREM
the parametric
representation
can be normalized to the (- 1,l) of (19) but with ok of (31).
n-cube in
6.2. The Jacobian inherent to (19) with (31) can be written us
(JI = {n! Vd2”AXd(l+
t~)}~~-l fi {f(l +
&j)}“-j
.
(32)
j=2
THEOREM r=
6.3. Let r denote the Euclidean 5
1 j=l
(Xj-
radius in (x1, . . . , x,) centered at node 1 by
112 ,
Xi)”
and let p be the nondimensional
(33) radius by
pA = $(1+ 61) . Then the nondimensional is identical to O(r). PROOF.
Substituting
(34) radius p is proportional to the Euclidean
radius r in the sense that O(p)
(34) into (31), (19) and (33) yield “2(x!-xj)M:
(35)
k=2
where M;
=
f(l -
52))
k=3,...,n,
M:={~~f(l+6)};(1-~k),
M ,*+I= fi 6(1+ j-2
5j)
*
(36)
302
M. Okabe, Gauss-Jacobi
quadrature for n-simplex
Evidently the MZ functions (k > 1) of (36) have no singularities and hence O(r) is identical to O(p). We then approximate
over the domain of interest,
a quantity # over n-simplex by fb” in the trial function form of
n+l
chh=
c &Nk.
(37)
k=l
Here & denotes the #-value at node k (& = +t), and Nk is the shape function taking the value unity at node k and zero at other nodes. We adopt the linear interpolation basis for the normalized n-cube but with collapse in the form of (18) such that N, = $(I - !5),
Nk = {
‘fi$(I+ tj)}$(l5)
7
k=2,...,n,
j=l
N n+l
=
fi
t(l+
(j)
(38)
.
j=1
THEOREM 4.4. Under the parametric representation of (19) with (31) the first derivatives of the trial function 4” of (37) with (38) hold O(r *-I) singularities within the vicinity of node 1. PROOF, Noting that the nondimensional radius p takes the value zero at node 1 and unity at all other nodes, it is easy to verify that our trial function of (37) with (38) can reproduce the ph term such that n+l
p* = c Nk.
(39)
k=2
The first derivatives of C;bhcan then be expected to have O(pAi’) singularities. In evaluating the singularities, it suffices to examine
since the shape functions Nk of (39) satisfy n+l
2 Nk = 1.
(41)
k=l
For simplicity but without any loss of generality, we place node 1 at the origin and node k + 1 on the &-axis distant by Xk from the origin (k = 1,. . . , n). Then the Jacobian transformation relations can explicitly be written as
303
M. Okabe, Gauss-Jacobi quadrature for n-simplex
and d[j/aXi
=
P_‘( (I-
[j)
i
Sil
-
(1 +
tj)Si,
j-l}/X
l=j
j=3,***,?2,
(jfi$(I+ !$I)} 7
1=2
i--l )..‘,
n.
(43)
It is clear in (42) and (43) that a&/L&, has the p”-’ term and @$/8Xi(j > 1) has the p-l term. However, aNJag (j > 1, k > 1) takes a value of strictly zero at l1 = -1, and hence it has a factor of $(l + &), i.e. ph. It is further noted that Nk (k > 2) has the term of nf1; $(l+ 6,). Since aNk/atj = 0 (j > k), we can conclude that aNk/aXi is a product of the p”-’ term and another polynomial. The proportionality of p and r in Theorem 6.2 then guarantees that the first derivatives of 4” hold O(rA-‘) singularities as expected. REMARK 6.5. The parametric representation of (19) with (31) is thus termed the semiradial singularity mapping. In two dimensions, systematic de~vation of (31) is given in 1141. REMARK 6.6. Under the semi-radial singularity mapping, any higher order polynomial interpolations can be applied to the trial function approximations, which naturally correspond to the Taylor expansion with respect to p*.
7. Gauss-Jacobi quadrature rules for n4mplex In finite-element formulations tegration may appear in the form
Then noting Theorem for the integrand by
under the semi-radial singularity mapping
with semi-radial
singularity
mapping
elements,
the in-
6.4, it suffices to consider the analytical integral operator I of (23) only
(45) Noting further Theorem
I[f(%
6.2, our analytical integral operator I for f of (45) can be written as
* - * W,+l)] = {IZ!V~2”A}11+‘“-2”A[fi(51)1fi In-jlfi(Si)l * 7
We now construct the approximate THEOREM
j-2
(46)
integral operator J, so as to satisfy (25) for I of (46).
7.1. For I of (46) with (45), let the Gauss-Jacobi
quadrature operator J,,, be
304
M. Okabe, Gauss-Jacobi quadrature
for n-simplex
constructed in the product form as
(47) Then J, realizes (25). REPARK
7.2. In the Gauss-Jacobi quadratures under the semi-radial the number of sampling points to realize (25) are also given by (28).
singularity mapping,
REMARK
7.3. In two dimensions, .!, by (47) under the semi-radial singularity mapping can be identified with that by (27) under the isoparametric mapping (except for A). The Gauss-Radau rule is thus always adequate for arbitrary A, which is of great significance in practical applications. REMARK
7.4. The semi-radial
singularity mapping is designed so that it is identical to the mapping at the limit of A = 1. Thus naturally substitution of A = 1 into (47) gives
isoparametric (27). REMARK
7.5. Let & in (31) and (38) be replaced
a1 + &>y= 4(1+
by [I through (48)
5,) .
Then (37) represents the non~lynomi~ trial function under the traditional polynomial mapping in terms of & and & to &. Here the Gauss-Jacobi quadrature operator .T, of (47) should be used by transforming the abscissa 5: into gl also through (48). For details of this alternative approach in fracture mechanics, see [17-211. REPARK
7.6. In the preceding
formulations,
we tacitly assume n > 1. In the case of n = 1,
(47) is simply written as Jm[f(w, @*)I = (w2~)J!?“[f1(51)1
*
(49)
Thus the semi-radial singularity Gauss-Jacobi quadratures for the l-simplex are valid only for h > $. Otherwise, the energy in the form of (44) could never be bounded.
8. Concluding
remarks
We have thus established the Gauss-Jacobi quadrature rules for the n-simplex as natural generalization of the Gauss-Radau ones for the 2-simplex. The Gauss-Jacobi quadrature operators for the n-simplex are given in a product of one-dimensional operators, derivation of which depends strongly on collapse of normalized n-cube. Evidently the Gauss-Jacobi quadrature is not optimal. The same order accuracy may be kept with less sampling points as is realized in two dimensions [lo, 111. Furthermore, it is
M. Okabe, Gauss-Jacobi quadrature
for n-simplex
305
aesthetically unpleasing since the sampling points are asymmetric with respect to vertices in the global Cartesian system. Some unpreferable bias may be produced due to asymmetry in approximate integrations. In finite-element applications, however, exact results can occasionally be obtained even through approximate quadrature operators as are demonstrated in this paper. The GaussJacobi quadratures could then be highly effective. It is further emphasized in pursueing much more desirable formulas that the number M of sampling points by (28) could be an upper bound to realize (25). The classical problem to find the optimal number, which should be at least lower than M of (28), has never been completely solved even for the Z-simplex [9]. Appendix A. Coeflkients of Jacobi polynomials We prove here that the polynomial of (5) with coefficients & of (4) is identified with the Jacobi polynomial. We express the Jacobi polynomial Pz(&) in the form of (5) such that
P”,B(&)=
k$o ak{i(l +&)I” .
The coefficient of the highest term in P?(&)
a,
= 2m*&m/2m
The orthogonal
condition
(A.1) is known as
.
64.2)
(9) further gives
W(1+ 51)/2YP?(5,)] = 0 , j = 0, * * . , m - 1 *
(A.3)
Substituting (A.l) into (A.3) and practising analytical integrations by the use of (6), we have (4). Appendix B. Zero points of Jacobi ~lynomi~s Unique existence of Gauss-Jacobi quadrature operators is guaranteed by the fact that the Jacobi ~lynomial ~~(~~) has ~-independent zero points in f-l, 13, to which we present a brief proof. The orthogonal condition (9) gives
~~~(~*~P~(~1)] = 0,
j = 0,. . . , m - 1.
W)
Let & (Jo = 1, . . . , n) denote n independent points in [-1, 1) where P”,B(&) changes its sign. Then evidently (n”,,l (& - ~~)~~~(~~) has no sign change in [-1,1). Suppose IZ< m, then (B.1) yields
03.2)
M. Okabe,
306
Gauss-Jacobi
quadrature
for n-simplex
Since the weight {i(l + &)}@is positive in [-1, l), (B.2) gives a contradiction that P”,p(,$,) should be thoroughly zero. Thus naturally II = m, and hence the Jacobi polynomial I??([,) has mindependent zero points in [-1, 1). Noting further P?(l)
= 1)
(B-3)
all zero points 5’1”are in E-1, 11.
Acknowledgment
The author is grateful to Prof. Norio Takenaka and suggestions.
of Nihon University for his useful comments
References [1] J.H. Argyris, K.E. Buch, I. Fried, G. Mareczek and D.W. Scharpf, Some new elements for matrix displacement methods, Proc. 2nd Conf. Matrix Methods in Structural Mechanics, Air Force Inst. Tech., Wright Patterson Air Force Base, OH, 1968. [2] J.H. Argyris and I. Fried, The LUMINA element for the matrix displacement method, Aero. J. 72 (1968) 514-517. [3] O.C. Zienkiewicz, The Finite Element Method (McGraw-Hill, London, 3rd ed., 1977). [4] M. Okabe, Y. Yamada and I. Nishiguchi, Reconsideration of rectangular Lagrange families with hierarchyranking bases, Comput. Meths. Appl. Mech. Engrg. 23 (1980) 369-390. [S] M. Okabe, Generalized Lagrange family for the cube with special reference to infinite-domain interpolations, Comput. Meths. Appl. Mech. Engrg. 38 (1983) 153-168. [6] J.H. Argyris, Triangular elements with linearly varying strain for the matrix displacement method, J. Roy. Aero. Sot. 69 (1965) 711-713. [7] J.H. Argyris, I. Fried and D.W. Scharpf, The TET20 and TEA8 elements for the matrix displacement method, Aero. J. 72 (1968) 618-625. [8] M. Okabe and N. Kikuchi, Some general Lagrange interpolations over simplex finite elements with reference to derivative singularities, Comput. Meths. Appl. Mech. Engrg. 28 (1981) l-25. [9] G. Strang and G.J. Fix, An Analysis of the Finite Element Method (Prentice-Hall, Englewood Cliffs, NJ, 1973). [lo] Radau, J. de Math. (1880) 283. [ll] G.R. Cowper, Gaussian quadrature formulas for triangles, Internat. J. Numer. Meths. Engrg. 7 (1973) 405-408. [12] T. Moan, Experiences with orthogonal polynomials and “best” numerical integration formulas on a triangle: with particular reference to finite element applications, Z. Angew. Math. Mech. 54 (1974) 501-508. [13] W. Magnus and F. Oberhettinger, Formeln und St&e fur die Speziellen Funktionen der Mathematischen Physik (Springer-Verlag, Berlin, 1948). [14] M. Okabe, Fundamental theory of the semi-radial singularity mapping with applications to fracture mechanics. Comput. Meths. Appl. Mech. Engrg. 26 (1981) 53-73. [15] M. Okabe, Quadrilateral element with semi-radial singularity mapping, Internat. J. Numer. Meths. Engrg. 17 (1981) 104.5-1067. [16] M. Okabe and N. Kikuchi, Semi-radial singularity mapping theory for line singularities in fracture mechanics, Comput. Meths. Appl. Mech. Engrg. 36 (1983) 257-276.
M. Okabe, Gauss-Jacobi [17] W.S. Blackbum,
[18] [19] [20] [21]
quadrature for n-simplex
307
Calculation of stress intensity factors, in: J.R. Whiteman, ed., The Mathematics of Finite Elements and Applications (Academic Press, London, 1973) 327-336. J.E. Akin, The generation of elements with singularities, Intemat. J. Numer. Meths. Engrg. 10 (1976) 12491259. J.E. Akin, Elements for problems with line singularities, in: J.R. Whiteman, ed., The Mathematics of Finite Elements and Applications, Vol. 3 (Academic Press, London, 1978) 65-75. M. Stern, Families of consistent conforming elements with singular derivative fields, Internat. J. Numer. Meths. Engrg. 14 (1979) 409-421. T.J.R. Hughes and J.E. Akin, Techniques for developing ‘special’ finite element shape functions with particular reference to singularities, Internat. J. Numer. Meths. Engrg. 15 (1980) 733-751.
GAUSSJACOBI
QUADRATURE
~PLI~A~ONS
TO ~
40 (1983) 293-307
RULES FOR n-SIMPLEX WITH FLUENT GOODS
Masayuki OKABE Mi~~~i Mining & Smelting Company Ltd., Niho~~~h~-MMro~~hi,
Chuo-Ku, Tokyo 103, Japan
Received 21 January 1983 Revised manuscript received 8 April 1983 For numerical integrations over the n-simplex, the Gauss-Jacobi quadrature operators are presented definitively. It is demonstrated in finite-element applications that the Gauss-Jacobi quadratures could be highly effective under the i~parametri~ mapping as well as under the semi-radial singul~ity mapping.
1. Introduction
Numerical integration plays an increasingly important role in numerical techniques such as the finite element and boundary element methods. A variety of finite element interpolations are known over the n-cube [l-S] as well as over the n-simplex [6-8], but we have some difficulty in approximate integrations to keep the required accuracy with necessary minimum numerical efforts. For the n-cube the product-Gauss rules [3,9] are fully explored. Here the operator is constructed as a product of Gauss-Legendre approximate integral operators defined on every axis. It is, indeed, optimal in the case of distorted cubes, although some redundancy exists for undistorted cubes. For the n-simplex, on the other hand, we have no complete devices. Only in two dimensions the so-called Gauss-Radau rules on triangles are fully known also in a product form [lo]. Gauss-Radau sampling points are located asymmetrically with respect to triangular vertices, and the same order accuracy could be attained with less sampling points. In fact, Cowper [ll] and Moan [12] simultaneously developed more effective formulas with symmetrical points but only in lower order accuracy. This paper concentrates on applications of Jacobi polynomials [13] to numerical integrations over the n-simplex. We first treat the isoparametric mapping, and Gauss-Jacobi quadrature operators are developed as natural generalization of Gauss-Radau ones. Then we define the semi-radial singularity mapping [8,14-16] in II dimensions, and the other kind of Gauss-Jacobi quadrature operators are presented.
2. Definite integral operators in one dimension
For the l-simplex operator I* by ~5-7825/83/$3.~
of (-1,l)
in the &-system, we define the analytical definite integral
@ 1983, Elsevier Science Publishers B.V. bosh-Holland)
M. Okabe, Gauss-Jucobi
294
quadrature for n-simplex
(1) Here @ is an appropriate real number (p > -1) and f denotes a sufficiently smooth function. Let Jc denote the approximate definite integral operator by
(2) Here tf denotes the 5, coordinate of a sampling point i, ~~~~) is the f-value at i and Wi designates the corresponding weight (i = 1, . . _, m). We then aim to integrate any polynomials up to 2m - 1 such that Jc [([$I = Ip [(&y] , PR~P~S~~O~
(3)
j = 0, . . . ,2nz - 1 .
2.1. Suppose that the system
m-i C ak/(P + k + I) = -2m +pCm/(/3+ m + I) ,
1 = 1,
. .
. >m
(4)
k=O
has the unique sohitions ak, and suppose that a polynomial by PY!%SI)= Z,+,Ct&
m-1 + ‘$I)>”f z: ak&l
+
(5)
6))”
k=O
has m-independent zero points. Here ~~+pl.$,= r(2m + /3 + l)/T(m + l)S(m + p + 1) with the gamma function r. Then (2) and (3) can be solved uniquely with respect to the weights Wi by selecting those zero points as [f. PROOF. Noting
(3) can explicitly be written as
2 {f(1+ [;))‘w,
= 2/(p+j+1),
j=O ,...,
2m-1.
(7)
i=l
Above-mentioned of the form
procedures
then guarantee
~{~(I+g;))“‘W=2/(p+j+I+2), i=l
with respect to Wi.
the identity of m systems with I = 0 to
j=O,...,m-1
m - 1
(81
295
M. Okabe, Gauss-Jacobi quadrature for n-simplex
3. Gauss-Jacobi quadrature operator in one dimension The polynomial Py(Q1) of (5) is known as the Jacobi polynomial in m degrees, see Appendix’ A. Let the monomials ((&)m)~=O orthogonalize :hrough the weight {i(l + &))B in (-1, l), then we have the Jacobi polynomials {P~(&))~=o. Here the orthogonal condition is written as (9) with nonnegative THEOREM satisfy (3).
delta S,,.
integers m and 12,and with the Kronecker
quadrature operator JB, can be uniquely defined so as to
3.1. The Gauss-Jacobi
PROOF. It is easy to prove that the Jacobi polynomial points in (-1, l), see Appendix B.
P?(&)
has m-independent
zero
We show the Gauss-Jacobi quadrature rules explicitly for some integer p in Tables 1 to 4. Zero points of PtB(&) and P?@(&) are plotted in Fig. 1 against @. REMARK
3.2. To obtain the formulas in closed form, we first make the Jacobi polynomial
-1
13* 0
5 1
2
3
4
5
6
Fig. 1. Zero points of Py8(&) and P$@(&) against j3.
M. Okabe, Gauss-Jacobi
296
Table 1 Gauss-Jacobi quadrature operators Ji which are identical to Gauss-Legendre ones. Computations of each abscissa and weight are done in double precision
m
Abscissa
1
0.0
2
3
4
5
51
Weight
W,
m
-0.57735027 0.57735027
1.0 1.0
1
0.33333333
-0.77459667 0.0 0.77459667
0.55555556 0.88888889 0.55555556
2
-0.28989795 0.68989795
0.36391724 0.63608276
-0.86113631 -0.33998104 0.33998104 0.86113631
0.34785485 0.65214515 0.65214515 0.34785484
-0.57531892 0.18106627 0.82282408
0.13965396 0.45848221 0.40186383
-0.90617985 -0.53846931 0.0 0.53846931 0.90617985
0.23692689 0.47862867 0.56888889 0.47862867 0.23692689
-0.72048027 -0.16718086 0.44631397 0.88579161
0.06236194 0.25969510 0.40692914 0.27101383
-0.80292983 -0.39092855 0.12405038 0.60397316 0.92038029
0.03149583 0.14781774 0.29277397 0.33434928 0.19356318
m
Abscissa
1
0.5
3
4
5
Table 2 Gauss-Jacobi quadrature operators Jk which are employed in Gauss-Radau rules. These operators can be applied also to the twodimensional semi-radial singularity mapping elements in linear fracture mechanics of A = l/2 as shown later
2.0
Table 3 Gauss-Jacobi quadrature operators J,?, which are used with Jz and J,!, in numerical integrations on tetrahedrons
2
quadrature for n-simplex
6;
Weight
Wi
3
4
5
Abscissa
5;
Weight
W,
1.0
Table 4 Gauss-Jacobi quadrature operators 54 which can be applied to the three-dimensional semiradial singularity mapping elements in the case of h = l/2 as shown later
m
Abscissa
0.66666667
1
0.6
0.5
-0.08830369 0.75497035
0.20157176 0.46509490
2
0.05971587 0.79742699
0.13381050 0.36618950
-0.41000442 0.30599247 0.85401195
0.05990141 0.29249254 0.31427272
-0.27347074 0.39762254 0.87584820
0.03295812 0.20919980 0.25784209
-0.59170284 -0.03409459 0.52279852 0.90299890
0.02070448 0.13726777 0.28691758 0.22177683
-0.47704442 0.07169289 0.58056646 0.91569416
0.00931673 0.08508345 0.21800874 0.18759108
-0.70210843 -0.26866695 0.22022723 0.65303936 0.93084212
0.00822765 0.06411120 0.17840032 0.25239792 0.16352957
-0.60757599 -0.16579958 0.29714001 0.69121030 0.93887141
0.00304138 0.03391465 0.12088991 0.20063301 0.14152106
3
4
5
51
Weight
Wi
M. Okabe, Gauss-Jacobi
explicit through the recurrence
297
quadrature for n-simplex
rule of the form
2(m + l)(m + P + 1)(2m + P)E?+&)
-(2m +p + 1Wm
+ P@m
+ P +2)&- ~"E!(&)
+ 2m (m + P)(2m + p + 2)PT_,(&) = 0 ,
m> 0.
(10)
Here E%)
= I,
W(51)
= w
(11)
+ 2) 51- PI *
Then we obtain its zero points by the Newton iteration referring to (I - &)DP?(&)
= (m + P + l)U?!(&)
(12)
- P??“(&))
where D = d/d&. REMARK
3.3. For the polynomial
f(51) = 2
f(!f~)h(tl)
+
f(tl) up to 2m - 1 degrees, we thus have
PZ(51)
mil
where cj denotes appropriate Wi =
coefficients. The weights Wi can then be obtained as
P[L*(&)] , i = 1,.
Here each interpolation Li(ei) =
(13)
Cj(5J
j=O
i=l
. . , m.
(14
function Li(t,) should be the polynomial
up to 2m - 1 satisfying
Sij a
(1%
It is especially remarked that the above-mentioned conditions allow unlimited selection of Li(tl)* If we adopt the Lagrangian interpolation polynomials with respect to 5: (i = 1, . . . , m) as usual, then naturally (15) is equivalent to (8) with 1= 0. REMARK 3.4. The Jacobi polynomial Pz(&) is identical Hence JZ is called the Gauss-Legendre quadrature operator.
to the Legendre
polynomial.
4. Isoparametric mapping between n-simplex and n-cube with collapse Consider an n-simplex in the global Cartesian system (x1, . . . , x,). Then the volume given by V, = det[“$ is1
“2 xj-l*q]/n! j=1
.
V,
is
(16)
298
M. Okabe, Gauss-Jacobi
quadrature for n-simplex
Here xi denotes the Xj coordinate of vertex i with xb = 1, and ei is the (n + 1)st unit vector having only the ith component of unity. The volume coordinate @k related to appex k can be defined as ~k=det[“~~(SkiX,-,+(l_4i)X;,)eie,]/n!V., i=l ;=I
k=l,...,n+l
(17)
where x0 = 1. Next we consider an n-cube normalized to (-1,l) in the local parametric system (tl, . . . , 5”). The subcube in II - 1 dimensions of [I = -1 is collapsed into node 1, and the subcube in n - 2 dimensions of (I = 1 and & = - 1 into node 2. In the same way, we collapse the subcube in it - 1 dimensions of [j = 1 (j < 1) and 4 = - 1 into node 1 up to I = n - 1. The residual two nodes with & = -1 and & = 1 are numbered II and II + 1, respectively. The collapsing procedures in two and three dimensions are, for clarity, illustrated in Figs. 2 and 3, respectively. Then the volume coordinates ok of (17) can be rewritten in terms of ,$I to & such that w1=
1(1- 61))
,i,k = { ‘fi
+(l - &) ,
$(I+
j=l
wn+1
=
fi
$(I
+
k = 2,
. . . , II,
(18)
5;.) -
i=l
THEOREM
4.1. The n-simplex
1
in (x1, . . . , x,)
to the (-1, 1) n-cube in
2
Fig. 2. Normalized 2-cube where node 1D is collapsed to node 1.
can be normalized
Fig. 3. Normalized 3-cube with collapse of nodes ID, ID’ and 1D” into node 1, and node 2D into node 2.
299
M. Okabe, Gauss-Jacobiquadrature for n-simplex
(51,* - * 5”) through the parametric representation 9
n+l
Xi=
J$X$Oky j=l,**.,
(19)
?Z.
k=l
Here
wk
is given by (18).
THEOREM
4.2. The Jacobian inherent to (19) with (18) can be writtenas
PROOF. It is clear in (18) that the summation
of all volume coordinates
takes unity;
(21) Then (19) can be expressed as n+l
xi-x$=
C (x!-xf)~~,
n.
j= l,...,
(2%
k=2
Differentiation
of (22) with respect to & to & thus gives (20).
REMARK 4.3. In usual finite-element applications, polynomial trial functions in terms of & to & are used. To the linear element, for example, volume coordinates of (18) are directly applied. Then the mapping function space and the trial function space have the same interpolation basis (@k), and hence the parametric representation of (19) with (18) could appropriately be called the isoparametric mapping.
5. Gauss-Jacobi quadrature rules for n-simplex under the isoparametric mapping Over the n-simplex we define the analytical integral operator I by (23) and the approximate
integral operator J,,, by
Jmlf(W,.*.,O n+1)1=v,
5 f(d,
i=l
. . . , d+,)W .
(24)
Here the number of sampling points is denoted by AX We aim to integrate ~lynomi~s up to 2m - 1 in terms of x1 to x, exactly by J,,,. Then
300
M. Okabe, Gauss-Jacobi
quadrature for n-simplex
noting (18) and (19) J,,, should realize
Jrn[fi(tj)k]=I[fi((j)k]T j=l
(25)
k=0,...,2m-1.
j=l
We restrict here f to monomials up to 2m - 1 in terms of volume coordinates. seperated in variables such that f(W1, .
* * 3 %+I)
=
fi
f;.(tj>
Then f can be
-
(26)
j=l
THEOREM 5.1. For f of (26), let the Gauss-Jacobi quadrature operator J,,, be constructed as a product of one-dimensional operators JE (p = 0, . . . , n - 1) by n-1
Jm [f
(w, . . . , ontl>l = (n! K/2”)
n J!Z[f n-dSn-d
.
(27)
p=o
Then J,,, realizes (25). REMARK
5.2.
In the Gauss-Jacobi quadrature rules, coordinates of sampling points are given in a tensor product of zero points of P?([,+) through p from 0 to n - 1. The weights are also given in a product of those for J!. Thus the number M of sampling points to realize (25) can be expressed as M=m”.
(28)
REMARK
5.3.
REMARK
5.4. The Gauss-Jacobi
REMARK
5.5. For
In order to construct the quadrature rules in the form of (24) sampling points should be written in volume coordinates by the use of (18). quadrature [lo], and hence called the Gauss-Radau ones.
rules for 2-simplex are established
by Radau
f of (26), the analytical integral operator I by (23) can also be separated in
variables such that
I[f(@l, . . . , @,+I)]= (n !VJ2”),@I”-‘[f;.Gfj>l . Here I”-’ is given by (1). Notice that fi in (29) could be nonpolynomials. (19) we can derive the formula of the form
(29) Then noting (18) and
(30) Here lj denotes the real number satisfying lj > -1. For nonnegative integers lj, (30) is well known especially in two and three dimensions
[3].
M. Okabe, Gauss-Jacobi
301
quadrature for n-simplex
6. Semi-radial singularity mapping between n-simplex and n-cube with collapse With the collapsing procedure alternately expressed as 01
=
1 - {$(1+ &)}lIA)
ok
=
I.%1
+
51))~‘~
previously argued, let the volume coordinates
02
=
G(1+
SW
tj)}!i(l -
(l@l%l+
[k)
O,Cl = {4(l+ 51)}l” { fi $(I +
[j)}
4(1-
7
52)
ok of (17) be
,
k = 3,. . . , n,
(31)
-
j=2
Here A is a positive constant below unity (0 < h < 1). THEOREM
6.1. The n-simplex in (x1, . . . , x,)
(51, * - * , &) through THEOREM
the parametric
representation
can be normalized to the (- 1,l) of (19) but with ok of (31).
n-cube in
6.2. The Jacobian inherent to (19) with (31) can be written us
(JI = {n! Vd2”AXd(l+
t~)}~~-l fi {f(l +
&j)}“-j
.
(32)
j=2
THEOREM r=
6.3. Let r denote the Euclidean 5
1 j=l
(Xj-
radius in (x1, . . . , x,) centered at node 1 by
112 ,
Xi)”
and let p be the nondimensional
(33) radius by
pA = $(1+ 61) . Then the nondimensional is identical to O(r). PROOF.
Substituting
(34) radius p is proportional to the Euclidean
radius r in the sense that O(p)
(34) into (31), (19) and (33) yield “2(x!-xj)M:
(35)
k=2
where M;
=
f(l -
52))
k=3,...,n,
M:={~~f(l+6)};(1-~k),
M ,*+I= fi 6(1+ j-2
5j)
*
(36)
302
M. Okabe, Gauss-Jacobi
quadrature for n-simplex
Evidently the MZ functions (k > 1) of (36) have no singularities and hence O(r) is identical to O(p). We then approximate
over the domain of interest,
a quantity # over n-simplex by fb” in the trial function form of
n+l
chh=
c &Nk.
(37)
k=l
Here & denotes the #-value at node k (& = +t), and Nk is the shape function taking the value unity at node k and zero at other nodes. We adopt the linear interpolation basis for the normalized n-cube but with collapse in the form of (18) such that N, = $(I - !5),
Nk = {
‘fi$(I+ tj)}$(l5)
7
k=2,...,n,
j=l
N n+l
=
fi
t(l+
(j)
(38)
.
j=1
THEOREM 4.4. Under the parametric representation of (19) with (31) the first derivatives of the trial function 4” of (37) with (38) hold O(r *-I) singularities within the vicinity of node 1. PROOF, Noting that the nondimensional radius p takes the value zero at node 1 and unity at all other nodes, it is easy to verify that our trial function of (37) with (38) can reproduce the ph term such that n+l
p* = c Nk.
(39)
k=2
The first derivatives of C;bhcan then be expected to have O(pAi’) singularities. In evaluating the singularities, it suffices to examine
since the shape functions Nk of (39) satisfy n+l
2 Nk = 1.
(41)
k=l
For simplicity but without any loss of generality, we place node 1 at the origin and node k + 1 on the &-axis distant by Xk from the origin (k = 1,. . . , n). Then the Jacobian transformation relations can explicitly be written as
303
M. Okabe, Gauss-Jacobi quadrature for n-simplex
and d[j/aXi
=
P_‘( (I-
[j)
i
Sil
-
(1 +
tj)Si,
j-l}/X
l=j
j=3,***,?2,
(jfi$(I+ !$I)} 7
1=2
i--l )..‘,
n.
(43)
It is clear in (42) and (43) that a&/L&, has the p”-’ term and @$/8Xi(j > 1) has the p-l term. However, aNJag (j > 1, k > 1) takes a value of strictly zero at l1 = -1, and hence it has a factor of $(l + &), i.e. ph. It is further noted that Nk (k > 2) has the term of nf1; $(l+ 6,). Since aNk/atj = 0 (j > k), we can conclude that aNk/aXi is a product of the p”-’ term and another polynomial. The proportionality of p and r in Theorem 6.2 then guarantees that the first derivatives of 4” hold O(rA-‘) singularities as expected. REMARK 6.5. The parametric representation of (19) with (31) is thus termed the semiradial singularity mapping. In two dimensions, systematic de~vation of (31) is given in 1141. REMARK 6.6. Under the semi-radial singularity mapping, any higher order polynomial interpolations can be applied to the trial function approximations, which naturally correspond to the Taylor expansion with respect to p*.
7. Gauss-Jacobi quadrature rules for n4mplex In finite-element formulations tegration may appear in the form
Then noting Theorem for the integrand by
under the semi-radial singularity mapping
with semi-radial
singularity
mapping
elements,
the in-
6.4, it suffices to consider the analytical integral operator I of (23) only
(45) Noting further Theorem
I[f(%
6.2, our analytical integral operator I for f of (45) can be written as
* - * W,+l)] = {IZ!V~2”A}11+‘“-2”A[fi(51)1fi In-jlfi(Si)l * 7
We now construct the approximate THEOREM
j-2
(46)
integral operator J, so as to satisfy (25) for I of (46).
7.1. For I of (46) with (45), let the Gauss-Jacobi
quadrature operator J,,, be
304
M. Okabe, Gauss-Jacobi quadrature
for n-simplex
constructed in the product form as
(47) Then J, realizes (25). REPARK
7.2. In the Gauss-Jacobi quadratures under the semi-radial the number of sampling points to realize (25) are also given by (28).
singularity mapping,
REMARK
7.3. In two dimensions, .!, by (47) under the semi-radial singularity mapping can be identified with that by (27) under the isoparametric mapping (except for A). The Gauss-Radau rule is thus always adequate for arbitrary A, which is of great significance in practical applications. REMARK
7.4. The semi-radial
singularity mapping is designed so that it is identical to the mapping at the limit of A = 1. Thus naturally substitution of A = 1 into (47) gives
isoparametric (27). REMARK
7.5. Let & in (31) and (38) be replaced
a1 + &>y= 4(1+
by [I through (48)
5,) .
Then (37) represents the non~lynomi~ trial function under the traditional polynomial mapping in terms of & and & to &. Here the Gauss-Jacobi quadrature operator .T, of (47) should be used by transforming the abscissa 5: into gl also through (48). For details of this alternative approach in fracture mechanics, see [17-211. REPARK
7.6. In the preceding
formulations,
we tacitly assume n > 1. In the case of n = 1,
(47) is simply written as Jm[f(w, @*)I = (w2~)J!?“[f1(51)1
*
(49)
Thus the semi-radial singularity Gauss-Jacobi quadratures for the l-simplex are valid only for h > $. Otherwise, the energy in the form of (44) could never be bounded.
8. Concluding
remarks
We have thus established the Gauss-Jacobi quadrature rules for the n-simplex as natural generalization of the Gauss-Radau ones for the 2-simplex. The Gauss-Jacobi quadrature operators for the n-simplex are given in a product of one-dimensional operators, derivation of which depends strongly on collapse of normalized n-cube. Evidently the Gauss-Jacobi quadrature is not optimal. The same order accuracy may be kept with less sampling points as is realized in two dimensions [lo, 111. Furthermore, it is
M. Okabe, Gauss-Jacobi quadrature
for n-simplex
305
aesthetically unpleasing since the sampling points are asymmetric with respect to vertices in the global Cartesian system. Some unpreferable bias may be produced due to asymmetry in approximate integrations. In finite-element applications, however, exact results can occasionally be obtained even through approximate quadrature operators as are demonstrated in this paper. The GaussJacobi quadratures could then be highly effective. It is further emphasized in pursueing much more desirable formulas that the number M of sampling points by (28) could be an upper bound to realize (25). The classical problem to find the optimal number, which should be at least lower than M of (28), has never been completely solved even for the Z-simplex [9]. Appendix A. Coeflkients of Jacobi polynomials We prove here that the polynomial of (5) with coefficients & of (4) is identified with the Jacobi polynomial. We express the Jacobi polynomial Pz(&) in the form of (5) such that
P”,B(&)=
k$o ak{i(l +&)I” .
The coefficient of the highest term in P?(&)
a,
= 2m*&m/2m
The orthogonal
condition
(A.1) is known as
.
64.2)
(9) further gives
W(1+ 51)/2YP?(5,)] = 0 , j = 0, * * . , m - 1 *
(A.3)
Substituting (A.l) into (A.3) and practising analytical integrations by the use of (6), we have (4). Appendix B. Zero points of Jacobi ~lynomi~s Unique existence of Gauss-Jacobi quadrature operators is guaranteed by the fact that the Jacobi ~lynomial ~~(~~) has ~-independent zero points in f-l, 13, to which we present a brief proof. The orthogonal condition (9) gives
~~~(~*~P~(~1)] = 0,
j = 0,. . . , m - 1.
W)
Let & (Jo = 1, . . . , n) denote n independent points in [-1, 1) where P”,B(&) changes its sign. Then evidently (n”,,l (& - ~~)~~~(~~) has no sign change in [-1,1). Suppose IZ< m, then (B.1) yields
03.2)
M. Okabe,
306
Gauss-Jacobi
quadrature
for n-simplex
Since the weight {i(l + &)}@is positive in [-1, l), (B.2) gives a contradiction that P”,p(,$,) should be thoroughly zero. Thus naturally II = m, and hence the Jacobi polynomial I??([,) has mindependent zero points in [-1, 1). Noting further P?(l)
= 1)
(B-3)
all zero points 5’1”are in E-1, 11.
Acknowledgment
The author is grateful to Prof. Norio Takenaka and suggestions.
of Nihon University for his useful comments
References [1] J.H. Argyris, K.E. Buch, I. Fried, G. Mareczek and D.W. Scharpf, Some new elements for matrix displacement methods, Proc. 2nd Conf. Matrix Methods in Structural Mechanics, Air Force Inst. Tech., Wright Patterson Air Force Base, OH, 1968. [2] J.H. Argyris and I. Fried, The LUMINA element for the matrix displacement method, Aero. J. 72 (1968) 514-517. [3] O.C. Zienkiewicz, The Finite Element Method (McGraw-Hill, London, 3rd ed., 1977). [4] M. Okabe, Y. Yamada and I. Nishiguchi, Reconsideration of rectangular Lagrange families with hierarchyranking bases, Comput. Meths. Appl. Mech. Engrg. 23 (1980) 369-390. [S] M. Okabe, Generalized Lagrange family for the cube with special reference to infinite-domain interpolations, Comput. Meths. Appl. Mech. Engrg. 38 (1983) 153-168. [6] J.H. Argyris, Triangular elements with linearly varying strain for the matrix displacement method, J. Roy. Aero. Sot. 69 (1965) 711-713. [7] J.H. Argyris, I. Fried and D.W. Scharpf, The TET20 and TEA8 elements for the matrix displacement method, Aero. J. 72 (1968) 618-625. [8] M. Okabe and N. Kikuchi, Some general Lagrange interpolations over simplex finite elements with reference to derivative singularities, Comput. Meths. Appl. Mech. Engrg. 28 (1981) l-25. [9] G. Strang and G.J. Fix, An Analysis of the Finite Element Method (Prentice-Hall, Englewood Cliffs, NJ, 1973). [lo] Radau, J. de Math. (1880) 283. [ll] G.R. Cowper, Gaussian quadrature formulas for triangles, Internat. J. Numer. Meths. Engrg. 7 (1973) 405-408. [12] T. Moan, Experiences with orthogonal polynomials and “best” numerical integration formulas on a triangle: with particular reference to finite element applications, Z. Angew. Math. Mech. 54 (1974) 501-508. [13] W. Magnus and F. Oberhettinger, Formeln und St&e fur die Speziellen Funktionen der Mathematischen Physik (Springer-Verlag, Berlin, 1948). [14] M. Okabe, Fundamental theory of the semi-radial singularity mapping with applications to fracture mechanics. Comput. Meths. Appl. Mech. Engrg. 26 (1981) 53-73. [15] M. Okabe, Quadrilateral element with semi-radial singularity mapping, Internat. J. Numer. Meths. Engrg. 17 (1981) 104.5-1067. [16] M. Okabe and N. Kikuchi, Semi-radial singularity mapping theory for line singularities in fracture mechanics, Comput. Meths. Appl. Mech. Engrg. 36 (1983) 257-276.
M. Okabe, Gauss-Jacobi [17] W.S. Blackbum,
[18] [19] [20] [21]
quadrature for n-simplex
307
Calculation of stress intensity factors, in: J.R. Whiteman, ed., The Mathematics of Finite Elements and Applications (Academic Press, London, 1973) 327-336. J.E. Akin, The generation of elements with singularities, Intemat. J. Numer. Meths. Engrg. 10 (1976) 12491259. J.E. Akin, Elements for problems with line singularities, in: J.R. Whiteman, ed., The Mathematics of Finite Elements and Applications, Vol. 3 (Academic Press, London, 1978) 65-75. M. Stern, Families of consistent conforming elements with singular derivative fields, Internat. J. Numer. Meths. Engrg. 14 (1979) 409-421. T.J.R. Hughes and J.E. Akin, Techniques for developing ‘special’ finite element shape functions with particular reference to singularities, Internat. J. Numer. Meths. Engrg. 15 (1980) 733-751.