Finite difference approximations for the three ... - Springer Link
Journal of Applied Mathematics and Physics (ZAMP) Vol. 33, March 1982
0044-2275/82/002266-16 $ 4.70/0 9 Birkh~iuser Verlag Basel, 1982
Finite difference approximations for the three-dimensional Laplacian in irregular grids By E. E. Okon, Engineering Analysis Unit, Faculty of Engineering, University of Lagos, Lagos, Nigeria
1. Introduction
A useful finite-difference expression for the two-dimensional Laplacian Operator at a point which is surrounded by six coplanar nodes was given sometime ago by Winslow [1]. Recently another numerical approximation for the operator was derived for the general case in which both the number and locations of the surrounding points are arbitrary [2, 3]. In the present paper the techniques which were used to obtain these two approximate expressions are extended to three-dimensions to give expressions of the Laplacian at a point in space in terms of the values which a scalar point function assumes at a number of surrounding points. These expressions can be used to obtain approximate numerical values of the Laplacian at nodes generated in threedimensional regions by automatic techniques such as that derived by Cook [4]. The procedure for obtaining the first numerical expression, which is the three-dimensional analogue of a two-dimensional technique presented by Okon [2], is based on the representation of the Laplacian as a ratio of a surface integral and the volume of the region bounded by the surface, which region contains the point at which the Laplacian is to be evaluated. The alternative expression relies, like Winslow's result for two dimensions [1], on the assumption that the scalar field is linearly dependent on certain of its nodal values. 2. A finite difference expression for the Laplacian
2.1 The divergence expression It is well known that if Po is any point lying within a region R enclosed by a piecewise - smooth surface S, then for a twice continuously-differentiable scalar field ~0: S~Vco'ndS V~ ~O[Vo= L i m ~ (1) V~O
V
'
Vol. 33, 1982
Finite difference approximations
267
where n is the unit vector along the outward normal to the smooth portions of S, and V is the volume of R. If n points Pi (r0, i = 1, 2 . . . . , n surround the point Po taken as the origin, at which VZ~0is to be evaluated, then at the point P~ (-~ r 0, on the plane perpendicular to P0 P~, the vector V~0 which occurs in (1) can be expressed approximately as: V~lr=~r~ - (~ - ~~176 n~
(2)
!r l
where r~ = P0 P~, n being the unit vector along P0 P~, and ~9o, ~0~are, respectively, the values of fp at P0 and P~. The points P~ can be considered as the vertices of a polyhedron surrounding Po. The planes which bisect the n rays P0 Pi perpendicularly will then intersect, in general, to give an auxiliary polyhedron the interior and bounding surface of which can be taken, respectively in (1) as the domain with volume V, and the surface S. By joining P0 to each of the points P~, the primary polyhedron which has these points as its vertices can be considered as a collection of M non-overlapping tetrahedra of which the j-th tetrahedron shown in Figure 1 is a typical example. In this figure the three planes which pass through the midpoints P]j, P~j, P~j of the three rays Po P1./, Po P2~ and Po Ps.~, and are such that, for example, the ray Po P2j is normal to the plane through P~j, meet in pairs along the lines R~ Q12~, R~ Q23~ and R~ Qalj. Here Qlzj is the circumcentre of the triangle Po P~j, P2j, etc., and R~ is the centre of the spherical surface on which Ply, P2~ and P3j lie.
go
p~js
Po
Componentof the auxiliary ~r~ polyhedron cor responding ~to the jth tetrahedron. \
Q3j~ . . . . . .
Figure 1 ff2j" '2j Thej-th tetrahedron associated with the first approximation for the Lapladan at P
268
E.E. Okon
ZAMP
2.2 The derivation of the generalised finite difference expression Assuming that the value of V~0 at Pfii given by (2) is the same for all the points on the plane polygonal segment ~z~j Dr,~t ,,~ ~g23i p t 2j of the bounding surface of the auxiliary polyhedron (Figure 1), then the contribution of the node P2j to the surface integral in (1) can be written as: ~I v ~ . n d S = r2i" {R~ x (r12~- ~2si)} (~02i- ~00) &, 2[r2~ ]2 '
(3)
in which S2j is the polygon Q ~ R~ Q23j P~i, and R i and rx2~ are the position vectors of R~, Q~i, etc. Also, the contribution of the pyramid Po Q12~ R5 Q2si P ~ to the volume V of the component of the auxiliary polyhedron corresponding to Po P~ P2j Psi is 1 V2j = ~ r2i" {Rj x (r12j - r23~)} 9 (4) By using (3) and (4) we find that the contribution of the tetrahedron Po Plj P2~ Psi to the expressions. ~ V e . n ds and V in (1) can be written as: 8
1 ff v ~ . n dS = ,2 r~ r 2 L'~22J,2sj ",j" {Ri x ('s~ - ',2~)} ( ~ - ~0) s~ 2 ~j 2i s~ + rs2i r~i r2j" {Ri x (ra~i - r23j)} (q72i- - ~0o) +
V~= ~
r~i rz2jrsi" {R~ x
[rzi" {Ri x (rs~i -
r12~)}
+ r s j " { R i X (r23~ --
(r23i - ralj)} ((~3i - ~o)]; +
r2~" {Ri x (r~2~
(5)
r2si)}
rsai)}]
where ~i is the value of ~0at Pi, etc.
(6)
In order to use the relations (5) and (6) to construct a finite difference expression for the Laplacian, it is necessary to know the n u m b e r of non-overlapping tetrahedra with a c o m m o n vertex which can be constructed with n points surrounding a given point such as Po. Starting with any tetrahedron, it is clear that any additional vertex contributes additionally 2 faces and 3 edges. It is then immediately apparent that the polyhedron having as its vertices the n points surrounding P, has 2 n - 4 faces and 3 n - 6 edges. The required number M of tetrahedra is, therefore, 2 n - 4. By considering all the M tetrahedra associated with Po we infer from (5) and (6), after some rearrangement, that the Laplacian at Po can be expressed in the form: 2~-4 6 ~--'~ Ji V2rp] p0 _ ~=1 (7) n--4
ZK~
Vol. 33, 1982
Finite difference approximations
269
where: Jj = ((p~r (Oo) (~- Rj) + ((fl2j-
(00) (~' " Rj) + (093j- (90) (~" 9Rj);
2 ( 6 ' R3 + r2j ( ~ " R3; xj = r~j (6" R3 + r~j Rj =
~=
~' -
r~j (rzj • r3j) + r22~(r3j • rlj) + r ~ (r:tj • r2j) .
2 rlj- (r2j • r3~)
( r ~ - *3j" r~j) (~3~ • ~j) § (,~j - r~" r2~) (r~ • r2~) ; 2 (r3~ x rlj) 2 2 (rlj x r2j) z 2 (rlj x r2j) 2 (r92j- r2j" r3j) (rzj • r3j) 2 (rzs x r3j) z
2 (r2j • r3j) 2
4
;
( r 2 j - r3j" rlj) (r3j • rlj) 2 (r3j x r~) z
3. An alternative finite difference approximation
3.I The analytical procedure In the j-th tetrahedron shown in Figure 2, the points P~j, P~j and P~j are, as in Figure 1, the midpoints of the rays Po P~, Po P2j and P0 P3j. In place of the circumcentres of the triangular faces P Plj Pz~, etc., we introduce their
zo-Po
/
~Componcnt of the auxiliary I ~ polyhedron cot responding
Figure 2 Thej-th tetrahedron associated with the alternative approximation for the Laplacian at P
270
E.E. Okon
ZAMP
centroids Q12j, Q23j and (23~j. The three planes determined by the triads {P~i, Q12~, Q31j} {P~j, Qz3j, Q12j} and {P~r Qaaj, 023i} meet, in general, at the point T}. Taking Po as the origin as before, with the position vector of Plj as raj, etc., the point T~y has the position Ti = -/1 ( r l j + r2r + r3j). By considering all the M tetrahedra having Po as a c o m m o n vertex, we obtain an auxiliary polyhedron with vertices P0, Pi~, Ql~i, P~i, 023~, P3i, 031i and T}. When applied to the auxiliary polyhedron generated by tetrahedra such as that in Figure 2, (1) gives an approximate expression for t h e Laplacian at P0 in the form: 2n--4 3
Z Z 2n--4
in which, for example, S~ is the quadrilateral segment Pie Qr~r T~ Q3xr of the surface of the auxiliary polyhedron and V~ is the volume of t h e j - t h component of the auxiliary polyhedron defined earlier. The vector areas Air of the S~i, i = 1, 2, 3, are given by: 24 A a~= ei + r2~ x rsi ,
24 A 2~= ~q x r3~ x rxi,
24 A 3r = ~i + rai x r2i,
where 9 j -- rli x r2i + r2~ x rsj + rsi • raj, and Vi = -~1 rl~- (r2i • rs~). We now introduce the volume coordinates r r/, ~,/2 of a general point Q (~, t/, ~,/2) interior to the tetrahedron Po Pli P2j P3j, which are defined by: = vol Q P2j Psi Po/V,
r/= vol Q Psi Psi P1/v,
= VO1 Q Pli P2j Po/If,
~ + t1 + ~ + 12 = 1,
where vol Q Pzi Psi Po is the volume of the tetrahedron Q P2i Psi Po and V is the volume of Po Pxi P2j P3i. It is evident that the position vector r of Q is given by
r= ~ r~ + tl rzi + ~ rsj
(9)
and the coordinates of the vertices of the tetrahedron are Ply(I,0,0,0),
P2j (0,1, 0,0), Ps~ (0,0, 1,0) and Po (0,0, 0,1). Under the assumption that the scalar field in the tetrahedron in Figure 2 is linearly dependent on the values ~o, ~lj, (a~j, (Psi, which it assumes at the nodes Po, Ply, P2~ and Psi, we find that its value (0 at Q is defined by: (o - q~o= (~1~ - qo) ~ + ( ~ 2 i - ~o) 1/+ (~Osi- ~o) r
(10)
Vol. 33, 1982
Finite difference approximations
271
From (9) we infer immediately that rl = (rs~ x
4 = (r2s x r3h " r/3j, = ( n j x r2s) 9 r / 3 s ,
ns) " r/As,
As = (r2S x rsj) " rlj,
(11)
so that V~ = (r2s x r3j)lds, etc.
(12)
The relations (10), (11) and (12) show that within the j-th tetrahedron the gradient of the scalar field assumes a constant value which, when applied in (8), gives an approximate expression for the Laplacian at P0 in the form: 2n--4
Vzq~lpo -
4 Z J} j=l 2n--4
(13)
ZK;
j=1
in which 2} = ~j- {(~o~ - ~oo) (r2j x r3j) + (e2j - ~oo) (r3j x rl;)
+ (~o3s- ~oo) ( n s x ~2s)}/K;, K; = (rlj • r2~) 9 rs~,
and, as before, 9 j = rlj x r2j + r2j x r3j + r3j • r~j.
3.2 Variational considerations The Poisson equation V2f0 = g,
(14)
where 9 = g (x, y, z) is a given scalar point function, is known to be the Euler Lagrange equation corresponding to the functional [5] J (~0) =SSS {(v~0)2 + 2 9 ~0} d12.
(15)
We now consider the first variation of this functional by taking as 12 the three-dimensional region which has been subdivided into tetrahedra such as that in Figure 2, and over which ~0 is defined as a linear function of its nodal values. We apply the Ritz procedure [5] by defining N trial functions fi (x, y, z) i - 1, 2 . . . . , N corresponding to the total number N of nodes in 12. These functions are such that for all the nodes Pk, k = 1 , 2 , . . . , N,
A (Pk) = 1, =0,
i= k
i#:k.
272
E.E. Okon
ZAMP
Taking the node Po in Figure 2 as a typical example, the restriction to this tetrahedron of the trial function corresponding to this node is thus the function # (x, y, z) = 1 - ~ (x, y, z) - r/(x, y, z) - ( (x, y, z), the functions ~, ~, being given by (11). The complete trial function associated with this node then has as its support the union of the M (= 2 n - 4) tetrahedra having Po as their common vertex. In the functional J (~0) defined by (15) we write N
(o = ~ ~kfk (X, y, Z),
(16)
k~l
i n which ~ok is the value assumed by (o at Pk. We arrange for J (~o) to be staOJ tionary by evaluating ~ at each node Pk and equating it to zero. N o w ~J = 2 ~ ! V(o-
(V~a) +
dO
,
(17)
where Ok is the support offk. The tetrahedron in Figure 2 can be considered as thej-th subdomain Okj of Ok, and for this subdomain (16) gives:
By considering all the subdomains O k j , j = 1, 2 , . . . , M, we find that (17) becomes 8J 2~-4 -----2 Z fS~{-V~~ do(lS) O~Pk
jffil ~k~
Equations (10), (11) and (12) s h o w t h a t V~0. (v~ + V ~ + V ~ ) = ~ " {(~0~- ~0o)(r2~ x r3j ) + (q)2j - ~00) (r3j x rlj) + ((P35 - ~00) ( r l j x
r2j)}/A~
where, following the notation in (13), we have written ~ok as ~Oo.The first term on the right-hand side of (18) can, therefore, be written as 1 , - ~ {V(o" (V~ + W / + VC)} d O = - -~- Jj. ~kj
(19)
Also, on writing g - V2CalPo using (14), we have:
g
dO - @ K; w lPo.
(20)
~kj ~" ~j Equations (18), (19) and (20), together wit the extremum condition ~ = 0, therefore give the relation: 2n--4
-
Z J}+--'l W(~ ~=1 4
which is identical to (13).
2n--4
Z K}-0, j=l
Vol. 33, 1982
Finite difference approximations
273
4. Applications to some configurations 4.1 The rectangular parallelepiped
Figure 3 shows a rectangular parallelepiped with edges 2 a, 2 b, 2c, these edges being parallel to the unit vectors i, j, k. The centroid O of the parallelepiped is the origin of the Cartesian system. Each face of the parallelepiped is associated with eight tetrahedra having O and the centroid of the face as common vertices. In the expression for the Laplacian given by (7), we note that for the parallelepiped, n = 26, so that the number of tetrahedra is 48. It can be verified that for the typical tetrahedral subdomain shown in Figure 4, the contributions to the J and K expressions in (7) are:
j: (~ol- ~Oo)b c 4a
,
K:
+
abc.
The seven other tetrahedra associated with A B C D contribute to J and K expressions which are identical to those already given, so that for this face: ~J
= 2 b c (~o, -
~oo)la,
~'~ K
(21)
= 2 a b c.
For the face B F G C we have, similarly: ZJ=2ca(~o~-~Oo)/b,
ZK=Zabc.
On considering all the faces in turn we have in (7): 48
Y'~Js = [b c (~o, + ~o3- 2 ~oo)la + c a (~o2 + ~o, - 2 q~o)lb + a b (~o5 + ~o8- 2 ~Oo)/C] Y=I 48
~Kj=12abc. H
D
G
I
",
I -
# I
/ /
'--~-
Figure 3 A, A subdivision of a rectangluar parallelepiped
-
-"K.
I
\
I 1
\ \1
1~ /~J /
-
/
l/
",~',j7--~-,~--~, /"..4'1.,~r
-
'
0044-2275/82/002266-16 $ 4.70/0 9 Birkh~iuser Verlag Basel, 1982
Finite difference approximations for the three-dimensional Laplacian in irregular grids By E. E. Okon, Engineering Analysis Unit, Faculty of Engineering, University of Lagos, Lagos, Nigeria
1. Introduction
A useful finite-difference expression for the two-dimensional Laplacian Operator at a point which is surrounded by six coplanar nodes was given sometime ago by Winslow [1]. Recently another numerical approximation for the operator was derived for the general case in which both the number and locations of the surrounding points are arbitrary [2, 3]. In the present paper the techniques which were used to obtain these two approximate expressions are extended to three-dimensions to give expressions of the Laplacian at a point in space in terms of the values which a scalar point function assumes at a number of surrounding points. These expressions can be used to obtain approximate numerical values of the Laplacian at nodes generated in threedimensional regions by automatic techniques such as that derived by Cook [4]. The procedure for obtaining the first numerical expression, which is the three-dimensional analogue of a two-dimensional technique presented by Okon [2], is based on the representation of the Laplacian as a ratio of a surface integral and the volume of the region bounded by the surface, which region contains the point at which the Laplacian is to be evaluated. The alternative expression relies, like Winslow's result for two dimensions [1], on the assumption that the scalar field is linearly dependent on certain of its nodal values. 2. A finite difference expression for the Laplacian
2.1 The divergence expression It is well known that if Po is any point lying within a region R enclosed by a piecewise - smooth surface S, then for a twice continuously-differentiable scalar field ~0: S~Vco'ndS V~ ~O[Vo= L i m ~ (1) V~O
V
'
Vol. 33, 1982
Finite difference approximations
267
where n is the unit vector along the outward normal to the smooth portions of S, and V is the volume of R. If n points Pi (r0, i = 1, 2 . . . . , n surround the point Po taken as the origin, at which VZ~0is to be evaluated, then at the point P~ (-~ r 0, on the plane perpendicular to P0 P~, the vector V~0 which occurs in (1) can be expressed approximately as: V~lr=~r~ - (~ - ~~176 n~
(2)
!r l
where r~ = P0 P~, n being the unit vector along P0 P~, and ~9o, ~0~are, respectively, the values of fp at P0 and P~. The points P~ can be considered as the vertices of a polyhedron surrounding Po. The planes which bisect the n rays P0 Pi perpendicularly will then intersect, in general, to give an auxiliary polyhedron the interior and bounding surface of which can be taken, respectively in (1) as the domain with volume V, and the surface S. By joining P0 to each of the points P~, the primary polyhedron which has these points as its vertices can be considered as a collection of M non-overlapping tetrahedra of which the j-th tetrahedron shown in Figure 1 is a typical example. In this figure the three planes which pass through the midpoints P]j, P~j, P~j of the three rays Po P1./, Po P2~ and Po Ps.~, and are such that, for example, the ray Po P2j is normal to the plane through P~j, meet in pairs along the lines R~ Q12~, R~ Q23~ and R~ Qalj. Here Qlzj is the circumcentre of the triangle Po P~j, P2j, etc., and R~ is the centre of the spherical surface on which Ply, P2~ and P3j lie.
go
p~js
Po
Componentof the auxiliary ~r~ polyhedron cor responding ~to the jth tetrahedron. \
Q3j~ . . . . . .
Figure 1 ff2j" '2j Thej-th tetrahedron associated with the first approximation for the Lapladan at P
268
E.E. Okon
ZAMP
2.2 The derivation of the generalised finite difference expression Assuming that the value of V~0 at Pfii given by (2) is the same for all the points on the plane polygonal segment ~z~j Dr,~t ,,~ ~g23i p t 2j of the bounding surface of the auxiliary polyhedron (Figure 1), then the contribution of the node P2j to the surface integral in (1) can be written as: ~I v ~ . n d S = r2i" {R~ x (r12~- ~2si)} (~02i- ~00) &, 2[r2~ ]2 '
(3)
in which S2j is the polygon Q ~ R~ Q23j P~i, and R i and rx2~ are the position vectors of R~, Q~i, etc. Also, the contribution of the pyramid Po Q12~ R5 Q2si P ~ to the volume V of the component of the auxiliary polyhedron corresponding to Po P~ P2j Psi is 1 V2j = ~ r2i" {Rj x (r12j - r23~)} 9 (4) By using (3) and (4) we find that the contribution of the tetrahedron Po Plj P2~ Psi to the expressions. ~ V e . n ds and V in (1) can be written as: 8
1 ff v ~ . n dS = ,2 r~ r 2 L'~22J,2sj ",j" {Ri x ('s~ - ',2~)} ( ~ - ~0) s~ 2 ~j 2i s~ + rs2i r~i r2j" {Ri x (ra~i - r23j)} (q72i- - ~0o) +
V~= ~
r~i rz2jrsi" {R~ x
[rzi" {Ri x (rs~i -
r12~)}
+ r s j " { R i X (r23~ --
(r23i - ralj)} ((~3i - ~o)]; +
r2~" {Ri x (r~2~
(5)
r2si)}
rsai)}]
where ~i is the value of ~0at Pi, etc.
(6)
In order to use the relations (5) and (6) to construct a finite difference expression for the Laplacian, it is necessary to know the n u m b e r of non-overlapping tetrahedra with a c o m m o n vertex which can be constructed with n points surrounding a given point such as Po. Starting with any tetrahedron, it is clear that any additional vertex contributes additionally 2 faces and 3 edges. It is then immediately apparent that the polyhedron having as its vertices the n points surrounding P, has 2 n - 4 faces and 3 n - 6 edges. The required number M of tetrahedra is, therefore, 2 n - 4. By considering all the M tetrahedra associated with Po we infer from (5) and (6), after some rearrangement, that the Laplacian at Po can be expressed in the form: 2~-4 6 ~--'~ Ji V2rp] p0 _ ~=1 (7) n--4
ZK~
Vol. 33, 1982
Finite difference approximations
269
where: Jj = ((p~r (Oo) (~- Rj) + ((fl2j-
(00) (~' " Rj) + (093j- (90) (~" 9Rj);
2 ( 6 ' R3 + r2j ( ~ " R3; xj = r~j (6" R3 + r~j Rj =
~=
~' -
r~j (rzj • r3j) + r22~(r3j • rlj) + r ~ (r:tj • r2j) .
2 rlj- (r2j • r3~)
( r ~ - *3j" r~j) (~3~ • ~j) § (,~j - r~" r2~) (r~ • r2~) ; 2 (r3~ x rlj) 2 2 (rlj x r2j) z 2 (rlj x r2j) 2 (r92j- r2j" r3j) (rzj • r3j) 2 (rzs x r3j) z
2 (r2j • r3j) 2
4
;
( r 2 j - r3j" rlj) (r3j • rlj) 2 (r3j x r~) z
3. An alternative finite difference approximation
3.I The analytical procedure In the j-th tetrahedron shown in Figure 2, the points P~j, P~j and P~j are, as in Figure 1, the midpoints of the rays Po P~, Po P2j and P0 P3j. In place of the circumcentres of the triangular faces P Plj Pz~, etc., we introduce their
zo-Po
/
~Componcnt of the auxiliary I ~ polyhedron cot responding
Figure 2 Thej-th tetrahedron associated with the alternative approximation for the Laplacian at P
270
E.E. Okon
ZAMP
centroids Q12j, Q23j and (23~j. The three planes determined by the triads {P~i, Q12~, Q31j} {P~j, Qz3j, Q12j} and {P~r Qaaj, 023i} meet, in general, at the point T}. Taking Po as the origin as before, with the position vector of Plj as raj, etc., the point T~y has the position Ti = -/1 ( r l j + r2r + r3j). By considering all the M tetrahedra having Po as a c o m m o n vertex, we obtain an auxiliary polyhedron with vertices P0, Pi~, Ql~i, P~i, 023~, P3i, 031i and T}. When applied to the auxiliary polyhedron generated by tetrahedra such as that in Figure 2, (1) gives an approximate expression for t h e Laplacian at P0 in the form: 2n--4 3
Z Z 2n--4
in which, for example, S~ is the quadrilateral segment Pie Qr~r T~ Q3xr of the surface of the auxiliary polyhedron and V~ is the volume of t h e j - t h component of the auxiliary polyhedron defined earlier. The vector areas Air of the S~i, i = 1, 2, 3, are given by: 24 A a~= ei + r2~ x rsi ,
24 A 2~= ~q x r3~ x rxi,
24 A 3r = ~i + rai x r2i,
where 9 j -- rli x r2i + r2~ x rsj + rsi • raj, and Vi = -~1 rl~- (r2i • rs~). We now introduce the volume coordinates r r/, ~,/2 of a general point Q (~, t/, ~,/2) interior to the tetrahedron Po Pli P2j P3j, which are defined by: = vol Q P2j Psi Po/V,
r/= vol Q Psi Psi P1/v,
= VO1 Q Pli P2j Po/If,
~ + t1 + ~ + 12 = 1,
where vol Q Pzi Psi Po is the volume of the tetrahedron Q P2i Psi Po and V is the volume of Po Pxi P2j P3i. It is evident that the position vector r of Q is given by
r= ~ r~ + tl rzi + ~ rsj
(9)
and the coordinates of the vertices of the tetrahedron are Ply(I,0,0,0),
P2j (0,1, 0,0), Ps~ (0,0, 1,0) and Po (0,0, 0,1). Under the assumption that the scalar field in the tetrahedron in Figure 2 is linearly dependent on the values ~o, ~lj, (a~j, (Psi, which it assumes at the nodes Po, Ply, P2~ and Psi, we find that its value (0 at Q is defined by: (o - q~o= (~1~ - qo) ~ + ( ~ 2 i - ~o) 1/+ (~Osi- ~o) r
(10)
Vol. 33, 1982
Finite difference approximations
271
From (9) we infer immediately that rl = (rs~ x
4 = (r2s x r3h " r/3j, = ( n j x r2s) 9 r / 3 s ,
ns) " r/As,
As = (r2S x rsj) " rlj,
(11)
so that V~ = (r2s x r3j)lds, etc.
(12)
The relations (10), (11) and (12) show that within the j-th tetrahedron the gradient of the scalar field assumes a constant value which, when applied in (8), gives an approximate expression for the Laplacian at P0 in the form: 2n--4
Vzq~lpo -
4 Z J} j=l 2n--4
(13)
ZK;
j=1
in which 2} = ~j- {(~o~ - ~oo) (r2j x r3j) + (e2j - ~oo) (r3j x rl;)
+ (~o3s- ~oo) ( n s x ~2s)}/K;, K; = (rlj • r2~) 9 rs~,
and, as before, 9 j = rlj x r2j + r2j x r3j + r3j • r~j.
3.2 Variational considerations The Poisson equation V2f0 = g,
(14)
where 9 = g (x, y, z) is a given scalar point function, is known to be the Euler Lagrange equation corresponding to the functional [5] J (~0) =SSS {(v~0)2 + 2 9 ~0} d12.
(15)
We now consider the first variation of this functional by taking as 12 the three-dimensional region which has been subdivided into tetrahedra such as that in Figure 2, and over which ~0 is defined as a linear function of its nodal values. We apply the Ritz procedure [5] by defining N trial functions fi (x, y, z) i - 1, 2 . . . . , N corresponding to the total number N of nodes in 12. These functions are such that for all the nodes Pk, k = 1 , 2 , . . . , N,
A (Pk) = 1, =0,
i= k
i#:k.
272
E.E. Okon
ZAMP
Taking the node Po in Figure 2 as a typical example, the restriction to this tetrahedron of the trial function corresponding to this node is thus the function # (x, y, z) = 1 - ~ (x, y, z) - r/(x, y, z) - ( (x, y, z), the functions ~, ~, being given by (11). The complete trial function associated with this node then has as its support the union of the M (= 2 n - 4) tetrahedra having Po as their common vertex. In the functional J (~0) defined by (15) we write N
(o = ~ ~kfk (X, y, Z),
(16)
k~l
i n which ~ok is the value assumed by (o at Pk. We arrange for J (~o) to be staOJ tionary by evaluating ~ at each node Pk and equating it to zero. N o w ~J = 2 ~ ! V(o-
(V~a) +
dO
,
(17)
where Ok is the support offk. The tetrahedron in Figure 2 can be considered as thej-th subdomain Okj of Ok, and for this subdomain (16) gives:
By considering all the subdomains O k j , j = 1, 2 , . . . , M, we find that (17) becomes 8J 2~-4 -----2 Z fS~{-V~~ do(lS) O~Pk
jffil ~k~
Equations (10), (11) and (12) s h o w t h a t V~0. (v~ + V ~ + V ~ ) = ~ " {(~0~- ~0o)(r2~ x r3j ) + (q)2j - ~00) (r3j x rlj) + ((P35 - ~00) ( r l j x
r2j)}/A~
where, following the notation in (13), we have written ~ok as ~Oo.The first term on the right-hand side of (18) can, therefore, be written as 1 , - ~ {V(o" (V~ + W / + VC)} d O = - -~- Jj. ~kj
(19)
Also, on writing g - V2CalPo using (14), we have:
g
dO - @ K; w lPo.
(20)
~kj ~" ~j Equations (18), (19) and (20), together wit the extremum condition ~ = 0, therefore give the relation: 2n--4
-
Z J}+--'l W(~ ~=1 4
which is identical to (13).
2n--4
Z K}-0, j=l
Vol. 33, 1982
Finite difference approximations
273
4. Applications to some configurations 4.1 The rectangular parallelepiped
Figure 3 shows a rectangular parallelepiped with edges 2 a, 2 b, 2c, these edges being parallel to the unit vectors i, j, k. The centroid O of the parallelepiped is the origin of the Cartesian system. Each face of the parallelepiped is associated with eight tetrahedra having O and the centroid of the face as common vertices. In the expression for the Laplacian given by (7), we note that for the parallelepiped, n = 26, so that the number of tetrahedra is 48. It can be verified that for the typical tetrahedral subdomain shown in Figure 4, the contributions to the J and K expressions in (7) are:
j: (~ol- ~Oo)b c 4a
,
K:
+
abc.
The seven other tetrahedra associated with A B C D contribute to J and K expressions which are identical to those already given, so that for this face: ~J
= 2 b c (~o, -
~oo)la,
~'~ K
(21)
= 2 a b c.
For the face B F G C we have, similarly: ZJ=2ca(~o~-~Oo)/b,
ZK=Zabc.
On considering all the faces in turn we have in (7): 48
Y'~Js = [b c (~o, + ~o3- 2 ~oo)la + c a (~o2 + ~o, - 2 q~o)lb + a b (~o5 + ~o8- 2 ~Oo)/C] Y=I 48
~Kj=12abc. H
D
G
I
",
I -
# I
/ /
'--~-
Figure 3 A, A subdivision of a rectangluar parallelepiped
-
-"K.
I
\
I 1
\ \1
1~ /~J /
-
/
l/
",~',j7--~-,~--~, /"..4'1.,~r
-
'
