- Email: [email protected]
Interpolation in Two and Three Dimensions
Chapter 6
Interpolation in two and three dimensions Multidimensional interpolation is commonly encountered in numerical methods such as the Finite Element Method (FEM) and the Finite Volume Method (FVM) used for solving partial differential equations. It is a general practice in numerical methods to discretize a two (three) dimensional domain into large number of small areas (volumes) known as elements in FEM and volumes in FVM. These methods assume a functional form for variable within each sub domain based on the nodal values and solve for the variables at the nodes. Development of multidimensional interpolation has also applications in computer graphics where surfaces are represented by Bezier curves and NURBS (Non-uniform rational B-splines). The present chapter intends to introduce the readers to multi-dimensional interpolation, with simple examples, to be made use of in later chapters.
255
256
6.1 6.1.1
Chapter.6 Interpolation in two and three dimensions
Interpolation over a rectangle Linear interpolation
The simplest interpolation technique is to use a linear variation along the two directions, say x and y for a function f (x, y). An engineering application will be in using table of properties of a pure substance such as water (referred to as Steam Tables). A part of such a table would look like this:
x1 x2
y1 f1 f2
y2 f3 f4
where the columns correspond to constant y’s, rows correspond to constant x’s and the entries correspond to the function f (x, y). We assume that the function varies linearly with both x and y. This means that, for a fixed y, the function varies linearly with x (along a column in the table). Suppose we desire to have the interpolated value of the tabulated function at (x i , yi ) within the range x1 < x i < x2 , y1 < yi < y2 defining a rectangular region (Figure 6.1) in the x, y plane.
f4 f1
fi
f3 y = y2 f2 y = yi y = y1
x=
x1
x=
xi
x=
x2
We hold y fixed at y1 and obtain the interpolated value of the function at x i , y1 by linear interpolation as (x i − x1 ) (x i − x2 ) f1 + f 2 = w1x f 1 + w2x f 2 f (x i , y1 ) = (x1 − x2 ) (x2 − x1 ) based on L 1 , Lagrange polynomial of degree 1. The weights have been identified with suitable subscripts. Similarly we hold y fixed at y2 and obtain the interpolated value of the function at x i , y2 by linear interpolation as f (x i , y2 ) =
(x i − x2 ) (x i − x1 ) f3 + f 4 = w1x f 3 + w2x f 4 (x1 − x2 ) (x2 − x1 )
Now we hold x fixed at x i and obtain the interpolated value of function at x i , yi as f (x i , yi )
= =
(yi − y2 ) (yi − y1 ) f (x i , y1 ) + f (x i , y2 ) (y1 − y2 ) (y2 − y1 ) w1y f (x i , y1 ) + w2y f (x i , y2 )
(6.1)
6.1. Interpolation over a rectangle
257
where again the weights have been identified by suitable subscripts. Specifically, x in the subscript means the weight is a function of x while y in the subscript means the weight is a function of y. Substituting Equations 6.1 and 6.1 in Equations 6.1 we get f (x i , yi )
=
w1x w1y f 1 + w2x w1y f 2 + w1x w2y f 3 + w2x w2y f 4
=
w1 f 1 + w2 f 2 + w3 f 3 + w4 f 4
(6.2)
where the final weights w1 - w4 are functions of both x and y. For example, we have w1 = w1x w1y = 1
(x i − x2 ) (yi − y2 ) (x1 − x2 ) (y1 − y2 )
(6.3)
1
w1x
w2x
0
0
3
4
y = y2
0
w2y
1 Figure 6.2: Graphical representa-
tion of weight functions
x = x2
x = x1
(x i , yi )
y = y1
1
2
0
w1y
1
The reader may similarly write down the other weights also. Thus the interpolated value is a weighted sum of the function at the 4 nodes. We see that the weights are products of Lagrange weights in one dimension. We also see that the weight is unity at (x1 , y1 ) and vanishes for x i = x2 or yi = y2 which means that this weight vanishes at the three nodes given by (x2 , y1 ), (x1 , y2 ) and (x2 , y2 ). The reader may verify, in a similar fashion, the properties of the other three weights also. Figure 6.2 illustrates the construction of weight or shape functions over a rectangular domain.
Example 6.1 Following data has been taken from a table of property values of a substance, after suppressing information on the units. x↓ y→ 15 25
500 1286 1285
600 1334.2 1333.6
Obtain the value of the function at x = 20, y = 530 using linear interpolation.
Solution : The weights are calculated with x i = 20, yi = 530 as given below. w1 =
(20 − 25)(530 − 600) = 0.35 (15 − 25)(500 − 600
w2 =
(20 − 15)(530 − 600) = 0.35 (25 − 20)(500 − 600
258
Chapter.6 Interpolation in two and three dimensions
w3 =
(20 − 25)(530 − 500) = 0.15 (15 − 25)(600 − 500
w4 =
(20 − 15)(530 − 500) = 0.15 (25 − 15)(600 − 500
The interpolated value of the function is then obtained as f (20, 530) = 0.35 × 1286 + 0.35 × 1285 + 0.15 × 1334.2 + 0.15 × 1333.6 = 1300.02
Linear interpolation; alternate way of representation Consider the weight w1 given in Equation 6.3. We may rewrite it in full as x i yi − x2 yi − y2 x i + x2 y2 w1 = = a 10 + a 11 x i + a 12 yi + a 13 x i yi x1 y1 − x2 y1 − y2 x1 + x2 y2 where x2 y2 x1 y1 − x2 y1 − y2 x1 + x2 y2 x2 a 12 = − x1 y1 − x2 y1 − y2 x1 + x2 y2 a 10 =
y2 x1 y1 − x2 y1 − y2 x1 + x2 y2 1 a 13 = x1 y1 − x2 y1 − y2 x1 + x2 y2 a 11 = −
(6.4)
It is easily shown all the other weights have the same form but with different coefficients (the a’s). It is thus seen that the interpolating function may be written in the alternate form f (x i , yi ) = c 0 + c 1 x i + c 2 yi + c 3 x i yi where the c’s are constants determined in following equation set to obtain the c’s. 1 x1 y1 1 x2 y1 1 x1 y2 1 x2 y2
(6.5)
terms of the four function values. In fact, we have the x1 y1 x2 y1 x1 y2 x2 y2
c0 c1 c 2 c3
f1 f2 = f 3 f4
(6.6)
These equations may easily be solved to get the vector c of coefficients in terms of the function values at the corners of the rectangle.
6.1.2
Local coordinate system for a rectangular element
A simple linear transformation may be used to transform the rectangular domain in x, y to a square of sides =2 in ξ, η such that −1 < ξ, η < 1 with the origin at the center of the square (see Figure 6.3). The required transformation may easily shown to be x1 + x2 x2 − x1 y1 + y2 y2 − y1 x= + ξ; y = + η (6.7) 2 2 2 2 Instead of interpolating in (x, y) we may interpolate in (ξ, η).
3(x1 , y2 )
30 (−1, 1)
4(x2 , y2 )
η
40 (1, 1) (ξ, η)
(x i , yi )
1(x1 , y1 )
ξ
=⇒
2(x2 , y1 )
10 (−1, −1)
Figure 6.3: Coordinate transformation
20 (1, −1)
6.1. Interpolation over a rectangle
259
Such a transformation is useful in engineering applications such as the finite element method. The given data may now be recast as follows: ξ = −1 ξ=1
η = −1 f1 f2
η=1 f3 f4
Equation set 6.6 is recast as 1 1 1 1
−1 1 −1 1
−1 −1 1 1
1 −1 −1 1
c0 c1 c2 c3
f1 f2 = f 3 f4
(6.8)
These equations may be solved easily (eg. Gauss elimination method) to get − f1 + f2 − f3 + f4 f1 + f2 + f3 + f4 c1 = 4 4 − f1 − f2 + f3 + f4 f1 − f2 − f3 + f4 c2 = c3 = (6.9) 4 4 We would now like to write the interpolating function in terms of weighted sum as given earlier. 1 − ξ − η + ξη We easily see that w1 (ξ, η) will be given by which may be recast in the form 4 (1 − ξ) (1 − η) w1 (ξ, η) = (6.10) 2 2 It is deliberately written in the form of a product of two Lagrange polynomials of first degree. We may show similarly that the following hold: (1 + ξ) (1 − η) (1 − ξ) (1 + η) w2 (ξ, η) = w3 (ξ, η) = 2 2 2 2 (6.11) (1 + ξ) (1 + η) w4 (ξ, η) = 2 2 Of course we may revert to (x, y) once the interpolation has been worked out in (ξ, η).
c0 =
Coordinate transformation
The advantage of coordinate transformations is that operations such as differentiation, integration and transforms can be applied on a local coordinate system rather than the actual coordinates. The mapping of any arbitrary quadrilateral to the local coordinate system can be achieved in a number of ways such as affine, projective transformations. The discussion on these methods are beyond the scope of the book and the interested reader should look into books on computational graphics for the same. Isoparametric transformations are most commonly employed in numerical methods such as FEM and FVM. The transformation maps line to line such that both variables ( f ) and the real coordinates (x, y) are evaluated using the same shape functions. For an arbitrary quadrilateral considered below, the coordinate transformation can be simply written as ξ 4 30 (−1, 1) l 2 40 (1, 1) l2 3 (x, y) l1 (η, ξ) l3 l1 l3 =⇒ 2 η l4 1
10 (−1, −1) l 4
20 (1, −1)
260
Chapter.6 Interpolation in two and three dimensions
x = w1 x1 + w2 x2 + w3 x3 + w4 x4 ;
y = w1 y1 + w2 y2 + w3 y3 + w4 y4
where w’s are the weight functions (evaluated in (ξ, η) coordinates) we are familiar with. If we closely observe, the coordinate transformation equation is of the form x = a 0 + a 1 ξ + a 2 η + a 3 ξη
y = b 0 + b 1 ξ + b 2 η + b 3 ξη
(6.12)
The two dimensional domain need not be represented by straight lines but could be any four arbitrary curves in which case coordinate transformation becomes murkier.
Example 6.2 Redo example 6.1 by transforming the coordinates and using the formulation given above
Solution : The transformations are easily obtained as follows. With x2 = 25 and x1 = 15, we have x1 + x2 15 + 25 x2 − x1 25 − 15 x − 20 = = 20; = = 5. Hence x = 5ξ + 20 or ξ = . Similarly we 2 2 2 2 5 y − 550 may write the second transformation as y = 50η + 550 or η = . Thus the interpolation 50 20 − 20 530 − 550 point corresponds to ξ = = 0 and η = = −0.4. The weights are then 2 50 calculated using Equations 6.10 and 6.11. w1
=
w2
=
w3
=
w4
=
(1 − 0)[1 − (−0.4)] 4 (1 + 0)[1 − (−0.4)] 4 (1 − 0)[1 + (−0.4)] 4 (1 + 0)[1 + (−0.4)] 4
= 0.35 = 0.35 = 0.15 = 0.15
The interpolated value is the same as f (0, 0.4) = 1300.02 obtained earlier in Example 6.1. Program 6.1: Program to calculate weights of a linear quadrilateral element
1 2 3 4 5 6 7 8 9
6.1.3
function w = weightquadlin (x , y ) % Input % Output
x,y : w :
l o c a l c o o r d i n a t e s of the point weights at point x , y n = size ( x ) ; % no o f p o i n t s w = zeros (n ,4) ; % i n i t i a l i z e weight w (: ,1) = (1 - x ) .*(1 - y ) /4; % c a l c u l a t e w e i g h t w (: ,2) = (1+ x ) .*(1 - y ) /4; % c a l c u l a t e w e i g h t w (: ,3) = (1 - x ) .*(1+ y ) /4; % c a l c u l a t e w e i g h t w (: ,4) = (1+ x ) .*(1+ y ) /4; % c a l c u l a t e w e i g h t
(s)
at at at at
( −1 , −1) (1 , −1) ( −1 ,1) (1 ,1)
Interpolating polynomials as products of ‘lines’
Now we are ready for another interpretation of the interpolating function we have obtained above. Consider Equation 6.10 as an example. It contains two factors, the first one containing 1−ξ the equation of a line l 1 (ξ) = which is a line lying in a plane that is normal to the ξ, η plane 2
6.1. Interpolation over a rectangle
261
that contains the ξ axis. This line passes through the points l 1 (−1) = 1 and l 1 (1) = 0. Similarly 1−η the second factor represents the equation of a line given by l 2 (η) = which is a line lying in a 2 plane that is normal to the ξ, η plane that contains the η axis. This line passes through the points l 2 (−1) = 1 and l 2 (1) = 0. Thus the weight is the product of equations of lines l 1 and l 2 . Similar interpretation is possible for the other weights also. What is important to notice is that the weight is a product of two lines. These lines represent Lagrange polynomials of first degree (see Figure 6.4).
wi l 1 or l 2
1
l 3 or l 4 Figure 6.4: Four `lines' used in linear
h = d/2
interpolation over a standard square
d
0
−1
ξ or η
−1
Using the ‘lines’ defined above we may recast the weights as w1 = l 1 l 2 ; w2 = l 2 l 3 ; w3 = l 1 l 4 and w4 = l 3 l 4
(6.13)
As indicated in the figure the height of the ‘line’ is equal to half the distance from the farthest point. Hence we may use the side of the square as a proxy for the corresponding line by simply halving d itself (i.e. by dividing the distance by the total distance, which is 2). This is shown schematically in Figure 6.5. The nodes are numbered as indicated and the lines are identified with the sides of the square domain. The weights are shown at the corresponding nodes as product of ‘lines’. This way of visualizing the weights in terms of lines is particularly useful in dealing with higher order interpolation and also interpolation over a triangular domain.
1−η 2
l1 =
l3 =
1+ξ 2
(ξ, η)
w1 = l 1 l 2 l = 1 + η 4 2
6.1.4
w4 = l 3 l 4
1−ξ 2
w3 = l 1 l 4 l 2 =
Figure 6.5:
Rectangular element
in local coordinates showing the lines
w2 = l 2 l 3
Lagrange Quadratic rectangular element
Figure 6.6 shows a 9 noded quadratic rectangular element also known as Lagrange quadratic element. The interpolating function for the element will be of the form f (ξ, η) = a 0 + a 1 ξ + a 2 η + a 3 ξη + a 4 ξ2 + a 5 η2 + a 6 ξη2 + a 7 ξ2 η + a 8 ξ2 η2
(6.14)
262
Chapter.6 Interpolation in two and three dimensions η
(−1, 1) 7
l2
(0, 1) 8
l3
(ξ, η)
(−1, 0) 4
l5
5 (0, 0)
(1, 1) 9 l1 6 (1, 0)
Figure 6.6:
ξ
Lagrange quadratic
rectangular element
l6 1 (−1, −1)
l4
2 (0, −1)
3 (1, −1)
The constants a i have to be determined such that f (ξ, η) is the same as the given data at all the nine nodes. Alternatively, the interpolation function can be estimated using Lagrange approach such that 9 X f (ξ, η) = w i (ξ, η) f i i =1
We can define two independent quadratic Lagrange polynomial functions along η and ξ directions. The total weight functions, following the procedure used in the case of linear interpolation, is the product of the two quadratic polynomial functions. The weight functions are of the form w(ξ, η) = (a 0 + a 1 ξ + a 2 ξ2 )(b 0 + b 1 η + b 2 η2 ) = F(ξ)G(η) Let us decompose the weight function and determine the weights of two separate functions F(ξ) and G(η). Then the weight functions at the nine points would be w1 = F−1 G −1 w4 = F−1 G 0 w7 = F−1 G 1
w 2 = F 0 G −1 w5 = F 0 G 0 w8 = F 0 G 1
w3 = F1 G −1 w6 = F 1 G 0 w9 = F 1 G 1
Based on Equation 5.17, the Lagrange weights L 2 (ξ) and L 2 (η) are written in vector form as FT (ξ)
=
GT (η)
=
ξ(ξ − 1) ξ(ξ + 1) ξ(ξ − 1) ξ(ξ + 1) = − (ξ − 1)(ξ + 1) (1 − ξ2 ) 2 2 2 2 η(η − 1) η(η + 1) η(η − 1) η(η + 1) 2 − (η − 1)(η + 1) = (1 − η ) 2 2 2 2
The nodal weights may then be written down as products GFT and are given as
=
w1 w2 w3 w w5 w6 = GFT 4 w7 w8 w9 2 ξη(ξ − 1)(η − 1) (1 − ξ )η(η − 1) 4 2 2 ξ(ξ − 1)(1 − η ) (1 − ξ2 )(1 − η2 ) 2 2 ξη(ξ − 1)(η + 1) (1 − ξ )η(η + 1) 4 2
ξη(ξ + 1)(η − 1) 4 2 ξ(ξ + 1)(1 − η ) 2 ξη(ξ + 1)(η + 1) 4
(6.15)
6.1. Interpolation over a rectangle
263
Alternate approach: product of lines The intersection points of the lines l 1 to l 6 are the nodes of the square domain. The nodes can be represented in terms of these lines. The equation of these lines are l 1 =⇒ 1 − ξ = 0;
l 2 =⇒ 1 − η = 0;
l 4 =⇒ 1 + η = 0;
l 5 =⇒ η = 0;
l 3 =⇒ 1 + ξ = 0; l 6 =⇒ ξ = 0
Node 1: The lines passing through node 1 are l 3 and l 4 . At node 1, all the weights except w1 are zero. In order to achieve this, all the weight functions except w1 should have either l 3 or l 4 as factors. Similarly, node 1 would be a factor of remaining lines l 1 , l 2 , l 5 and l 6 . Therefore the weight function for node 1 can be written down as w1 =
ξη(ξ − 1)(η − 1) l1 l2 l5 l6 = l 1 (−1, −1)l 2 (−1, −1)l 5 (−1, −1)l 6 (−1, −1) 4
which is the same as that obtained earlier. Similar approach can be adopted to derive weight functions for the other nodes.
Example 6.3 Obtain value of the tabulated function at x = 0.14, y = 0.07 by quadratic interpolation. x↓ y→ 0 0.1 0.2
0 1.000000 1.010050 1.040811
0.1 1.010050 1.020201 1.051271
0.2 1.040811 1.051271 1.083287
Use a nine noded element to solve the problem.
Solution : We transform the domain to the standard square domain by the transformation ξ = 10(x − 0.1), η = 10(y − 0.1). The elements are numbered as shown in Figure 6.6. In the present example, the function is required at x = 0.14 which corresponds to ξ = 10(0.14 − 0.1) = 0.4 and y = 0.07 which corresponds to η = 10(0.07 − 0.1) = −0.3. Substituting these in the matrix of weights (Equation 6.15) we get −0.0234 0.1638 0.0546 −0.1092 0.7644 0.2548 w= 0.0126 −0.0882 −0.0294 The interpolated value of the function is then obtained as sum of product of weights and the function values. We thus obtain f (ξ = 0.4, η = −0.3) = f (x = 0.14, y = 0.07) =
9 X
w i f i = 1.024826
i =1
The reader may verify that linear interpolation using only the corner nodes would give a 2 2 value for this as 1.043259. The function used for generating the data was f (x, y) = e(x + y ) which has a value of f (0.14, 0.07) = 1.024803. The quadratic interpolated value has an error of 0.000024 while the linear interpolated value has an error of 0.018457 ≈ 0.019, with respect to the exact value. There is a considerable improvement in using quadratic interpolation over linear interpolation in the present case.
264
Chapter.6 Interpolation in two and three dimensions η
6(−1, 1)
7(0, 1)
l2
8(1, 1) (ξ, η)
l3
l6
l5
l1 5(1, 0) ξ
(−1, 0)4
Figure 6.7:
Eight noded rectan-
gular element
l7
l4
1(−1, −1)
6.1.5
l8
2(0, −1)
3(1, −1)
Quadratic eight noded rectangular element
The above element is also known as ‘serendipity’ quadratic element.1 The interpolation function in the domain is of the form f (ξ, η) = a 0 + a 1 ξ + a 2 η + a 3 ξη + a 4 ξ2 + a 5 η2 + a 6 ξη2 + a 7 ξ2 η
(6.16)
The present domain does not have a point at the origin. Therefore the present interpolation formula does not support the term involving ξ2 η2 . Unlike the previous example, the weight function of the present interpolation formula cannot be simplified as product of two univariate functions in ξ and η. Let us adopt the Lagrange interpolation approach based on product of lines. The nodes of the domain are the points of intersection of lines l 1 to l 8 . Hence these lines can be used to represent the nodes. Equations of lines indicated in Figure 6.7 are given as a table below. Equations of lines: Line l1 l2 l3 l4
Formula 1−ξ 1−η 1+ξ 1+η
Line l5 l6 l7 l8
Formula ξ+η−1 −ξ + η − 1 −ξ − η − 1 ξ−η−1
Calculation of nodal weights as product of lines: For node 1 the procedure is as follows: l 1 and l 2 must be factors of w1 as these will be zero at nodes 3,5,6,7 and 8. The function should be equal to zero at nodes 2 and 4 as well. This means l 7 is also a factor of the weight function. Therefore the weight function can be written as w1 =
l1 l2 l7 (1 − ξ)(1 − η)(−ξ − η − 1) = l 1 (−1, −1)l 2 (−1, −1)l 7 (−1, −1) 4
Similarly we treat other nodes to arrive at weight functions as given below. 1
For a serendipity element function values are specified only at points on the boundary.
6.1. Interpolation over a rectangle
Node Number 1 2 3 4 5 6 7 8
265
Weight formula in terms of product of lines (l 1 × l 2 × l 7 )/4 (l 1 × l 2 × l 3 )/2 (l 2 × l 3 × l 8 )/4 (l 1 × l 2 × l 4 )/2 (l 2 × l 3 × l 4 )/2 (l 1 × l 4 × l 6 )/4 (l 1 × l 3 × l 4 )/2 (l 3 × l 4 × l 5 )/2
Weight expressions (1 − ξ)(1 − η)(−ξ − η − 1)/4 (1 − ξ)(1 − η)(1 + ξ)/2 (1 − η)(1 + ξ)(ξ − η − 1)/4 (1 − ξ)(1 − η)(1 + η)/2 (1 − η)(1 + ξ)(1 + η)/2 (1 − ξ)(1 + η)(−ξ + η − 1)/4 (1 − ξ)(1 + ξ)(1 + η)/2 (1 + ξ)(1 + η)(ξ + η − 1)4
Note: One may be tempted to use lines η = 0 and ξ = 0 instead of l 5 to l 8 . However, such an interpolation polynomial would not be consistent with Equation 6.16.
Example 6.4 A certain function is available at eight points as given in the following table. Obtain value of the function at x = 0.14, y = 0.07 by quadratic interpolation. x↓ y→ 0 0.1 0.2
0 1.000000 1.010050 1.040811
0.1 1.010050 ··· 1.051271
0.2 1.040811 1.051271 1.083287
Solution : Note that this set of data is the same as that used in Example 6.3 obtained by removing the data at the origin (removing data at point 5 in Figure 6.6). Nodes are numbered as shown in Figure 6.7. In the present example interpolated value is desired at ξ = 0.4, η = −0.3. We obtain all the weights using the weight functions and tabulate the results. Node No. 1 2 3 4 5 6 7 8 Sum =
f (x, y) 1.000000 1.010050 1.040811 1.010050 1.051271 1.040811 1.051271 1.083287 ···
W ei ght -0.2145 0.546 -0.1365 0.273 0.637 -0.1785 0.294 -0.2205 1
W ei ght × f (x, y) -0.214500 0.551487 -0.142071 0.275744 0.669660 -0.185785 0.309074 -0.238865 1.024744
Bold entry in the table is the required value of the function at the interpolation point.
266
6.2
Chapter.6 Interpolation in two and three dimensions
Interpolation over a triangle
Interpolation over a triangle is an important exercise since most numerical methods for solving partial differential equations use triangular elements. Consider a triangular element shown in Figure 6.8.
3 (x3 , y3 )
l2
(x, y)
l1 Figure 6.8: Linear triangular element
1 (x1 , y1 )
l3
2 (x2 , y2 )
l 1 , l 2 and l 3 are the equations of the lines forming the edges of the triangle. As there are three points, the interpolating function is of the form2 f (x, y) = a 0 + a 1 x + a 2 y
(6.17)
The constants a 0 , a 1 and a 2 can be determined by solving the system of equations at the three nodes of the triangle. Alternatively, we can apply the Lagrange interpolation scheme for the 3 X triangle such that f (x, y) = wi f i . i =1
Consider node 1. The weight w1 should be equal to 1 at node 1 and should be equal to zero at other two nodes. l 2 and l 3 are equal to zero at node 1. Therefore the weight function w1 is proportional to l 1 and is of the form given by Equation 6.17. Meeting all the required conditions, w1 can be written as l 1 (x, y) w1 = l 1 (x1 , y1 ) In general, the weight function at node i can be written as wi =
l i (x, y) l i (x i , yi )
(6.18)
Local coordinate system for a triangular element l i (x, y) = η i . η i is nothing but the ratio of normal distance of point (x, y) from l i (x i , yi ) line l i to the normal distance of node i from the line. When the point lies on l i , η i = 0 and when the point is at node i, η i = 1. Hence, the following representation also serves as an alternative coordinate system for the triangular element with η i as the basis. Any point in (x, y) coordinates can be represented in terms of (η 1 , η 2 , η 3 ) as indicated in the following figure. Let the fraction
2
Number of nodes less than four would mean that there can be only three terms in the interpolating function. The x y term would be missing.
6.2. Interpolation over a triangle
267
3 (x3 , y3 )
3 (0, 0, 1)
(η 1 , η 2 , η 3 )
η2
(x, y)
η1 η3
1 (x1 , y1 )
2 (x2 , y2 )
1 (1, 0, 0)
2 (0, 1, 0)
Figure 6.9: Transformation of coordinates
Hence the equation of the three lines l 1 , l 2 and l 3 in local coordinate system would be η 1 = 0, η 2 = 0 and η 3 = 0. For a linear triangular element, w i = η i The area of a triangle is proportional to product of height and base. Hence η i is also the ratio of the area of the triangles made by points (x, y) and (x i , yi ) with line l i . This can be understood from Figure 6.10.
3 (x3 , y3 )
l2
l1 Figure
(x, y)
6.10:
Calculation
of
weight functions for a triangular element
1 (x1 , y1 )
2 (x2 , y2 )
l3
η 3 is the ratio of the area of the shaded triangle to the area of the triangle itself. The weight n X function for a linear triangular element can be generalized as w i = η i = A i /A. Since A i = A, i =1
sum of weights at any point is equal to 1 i.e. given by the following determinant ¯ ¯ 1 1 ¯¯ A3 = ¯ 1 2¯ 1
x1 x2 x
y1 y2 y
3 X
w i (x, y) = 1. The area of the shaded triangle is
i =1
¯ ¯ ¯ 1 ¯ = { x(y1 − y2 ) + (x2 − x1 )y + (x1 y2 − x2 y1 )} ¯ 2 ¯
When (x, y) lies on the line l 3 , A 3 = 0 and when (x, y) = (x3 , y3 ), A 3 = A, where A is the area of the triangular element given by ¯ ¯ ¯ 1 x1 y1 ¯ ¯ 1 ¯¯ A = ¯ 1 x2 y2 ¯¯ 2¯ 1 x3 y3 ¯
268
Chapter.6 Interpolation in two and three dimensions
η i is of the form (a i + b i x + c i y) where a i ,b i and c i are constants and are evaluated as follows
(x2 y3 − y2 x3 ) 2A (x3 y1 − y3 x1 ) a2 = 2A (x1 y2 − y1 x2 ) a3 = 2A a1 =
(y2 − y3 ) 2A (y3 − y1 ) b2 = 2A (y1 − y2 ) b3 = 2A b1 =
(x3 − x2 ) 2A (x1 − x3 ) c1 = 2A (x2 − x1 ) c1 = 2A c1 =
(6.19)
Interpretation using local coordinates: Any local coordinate can be represented by other two local coordinates i.e. η 3 = 1 − η 1 − η 2 . The right angled triangle on η 1 ,η 2 plane in Figure 6.11 represents the map of the actual triangle. In a sense, the transformation exercise can be looked upon as mapping triangle onto a right angled isosceles triangle. MATLAB program has been given
η2
2(0, 1) η2 = 1 − η1
Figure
6.11:
Representation of local
coordinates of a triangle
η1
3(0, 0)
1(1, 0)
below to calculate the coefficients for a triangular element Program 6.2: Calculate coe cients for local coordinates of a triangular element
1 function [a ,b ,c , A ] = trlocalcoeff ( p1 , p2 , p3 ) 2 % I n p u t p1 , p2 , p 3 : c o o r d i n a t e s o f t h e v e r t i c e s 3 4 5 6 7 8 9 10 11 12 13 14 15 16
p1 ( 1 ) = x1 ; p1 ( 2 ) % = y1 % p2 ( 1 ) = x2 ; p2 ( 2 ) = y2 ; p3 ( 1 ) = x3 ; p3 ( 2 ) = y3 % Output a,b, c : c o e f f i c i e n t s η = a + bx + c y % A : area of t r i a n g l e a (1) = p2 (1) * p3 (2) - p3 (1) * p2 (2) ; % c a l c u l a t e a
a (2) = p3 (1) * p1 (2) - p1 (1) * p3 (2) ; a (3) = p1 (1) * p2 (2) - p2 (1) * p1 (2) ; A = ( a (1) + a (2) + a (3) ) /2; a = a /(2* A ) ; b (1) = ( p2 (2) - p3 (2) ) /(2* A ) ; b (2) = ( p3 (2) - p1 (2) ) /(2* A ) ; b (3) = ( p1 (2) - p2 (2) ) /(2* A ) ; c (1) = ( p3 (1) - p2 (1) ) /(2* A ) ; c (2) = ( p1 (1) - p3 (1) ) /(2* A ) ; c (3) = ( p2 (1) - p1 (1) ) /(2* A ) ;
% A = 2* area o f % normailze a % calculate b
% calculate
triangle
c
The following Matlab function can be used to convert from actual coordinates to local coordinates. Program 6.3: Convert from
(x, y)
to local coordinates
1 function eta = trianglelocalcord (a ,b ,c ,x , y ) 2 % Input a , b , c : c o e f f i c i e n t s η = a + bx + c y 3 % x , y : coordinates of points
(η 1 , η 2 , η 3 )
of a triangle
6.2. Interpolation over a triangle
4 5 6 7 8 9
% Output
eta
:
269
local
coordinates
of
points
x,y
n = size (x ,1) ; eta = zeros (n ,3) ; for i =1:3 eta (: , i ) = a ( i ) + b( i ) * x + c ( i ) * y ; end
Quadratic triangular element Now we derive expressions for weight functions for a quadratic triangular element.
3
(x, y)
5
4
Figure 6.12:
Quadratic triangular ele-
ment
6
1
2
There are six points. Three nodes occupy the vertices of the triangle and three points are located at the midpoints of the three edges of the triangular element. The interpolation function would be of the form3 f (x, y) = a 0 + a 1 x + a 2 y + a 3 x y + a 4 x2 + a 5 y2 (6.20) l 1 , l 2 and l 3 are the equations of the lines forming the edges of the triangle. l 4 , l 5 and l 6 are the
3
l6
6
5
l2
l3
6.13:
Lines
representing
quadratic triangular element
l1 l5
l4 1
Figure
4
2
equation of the lines joining the midpoints of the edges of the triangle. Let us represent the nodes in terms of the triangular coordinate system introduced earlier. 1 =⇒ (1, 0, 0) 4 =⇒ (0.5, 0.5, 0)
2 =⇒ (0, 1, 0) 5 =⇒ (0, 0.5, 0.5)
3 =⇒ (0, 0, 1) 6 =⇒ (0.5, 0, 0.5)
(6.21)
From Equation 6.21, the equation of the lines can be written down as l 1 =⇒ η 1 = 0 3
l 2 =⇒ η 2 = 0
l 3 =⇒ η 3 = 0
With only six points available the quadratic is constrained to the form shown here. Recall that nine points gave the more general expression (Equation 6.14)
270
Chapter.6 Interpolation in two and three dimensions l 4 =⇒ η 1 − 0.5 = 0
l 5 =⇒ η 2 − 0.5 = 0
l 6 =⇒ η 3 − 0.5 = 0
Now, we derive the weight function for all the nodes. Consider node 1: The weight w1 should be equal to 1 at node 1 and should be equal to zero at the other nodes. Lines l 2 and l 3 are zero where as l 1 = 1 at node 1. Therefore, w1 should be proportional to l 1 . This means w1 = 0 at nodes 2, 3 and 5. As w1 should be zero at nodes 4 and 6, w1 should be proportional to l 4 also. Hence, meeting all the required conditions w1 can be written as l 1 (η 1 , η 2 , η 3 )l 4 (η 1 , η 2 , η 3 ) = 2η 1 (η 1 − 0.5) w1 = l 1 (1, 0, 0)l 4 (1, 0, 0) Consider Node 4: As node 4 is a point on l 3 , η 3 = 0. However, lines l 1 and l 2 are not equal to zero at node 4. This means w4 is proportional to l 1 l 2 . Therefore w4 = 0 at other nodes except node 4. Meeting all the required conditions w4 can be written as w4 =
l 1 (η 1 , η 2 , η 3 )l 2 (η 1 , η 2 , η 3 ) = 4η 1 η 2 l 1 (0.5, 0.5, 0)l 2 (0.5, 0.5, 0)
The weights for all other points may be similarly determined and are summarized below w1 = 2η 1 (η 1 − 0.5); w4 = 4η 1 η 2 ;
w2 = 2η 2 (η 2 − 0.5); w5 = 4η 2 η 3 ;
w3 = 2η 3 (η 3 − 0.5); w6 = 4η 3 η 1
Example 6.5 Function value is available at six points as shown in the table below. Point No. 1 4 2 5 3 6
x 0.2 0.25 0.3 0.325 0.35 0.275
y 0.6 0.35 0.1 0.4 0.7 0.65
f (x, y) 2.2255 1.8221 1.4918 2.0647 2.8577 2.5219
Obtain the value of the function at x = 0.28, y = 0.5 using linear and quadratic interpolation.
Solution : Figure 6.14 shows the node numbering system used in the analysis. Linear interpolation Nodes at the apex of triangle are used for linear interpolation. We evaluate all the 1 weights using the determinants presented in the text. Factor is dropped in all the area 2 calculations since the weights are normalized with the area of the triangle. The area of the triangle is proportional to A given by ¯ ¯ ¯ 1 0.2 0.6 ¯¯ ¯ 0.1 ¯¯ = 0.085 A = ¯¯ 1 0.3 ¯ 1 0.35 0.7 ¯ Weights are now calculated as given below. Weight are given by the following. ¯ ¯ ¯ 1 0.3 0.1 ¯¯ 1 ¯¯ 1 0.35 0.7 ¯¯ = 0.3765 η 1 = w1 = 0.085 ¯¯ 1 0.28 0.5 ¯
6.2. Interpolation over a triangle
271
0.8 3
0.7 0.6
6 1
Figure
(0.28,0.5)
0.5
y
of
for
Numbering linear
and
quadratic interpolation over a triangle showing the inter-
5
0.4
6.14:
nodes
4
polation point
0.3 0.2 2
0.1 0 0.15
0.2
0.25
x
0.3
0.35
¯ ¯ 1 1 ¯¯ 1 η 2 = w2 = 0.085 ¯¯ 1 ¯ ¯ 1 1 ¯¯ 1 η 3 = w3 = 0.085 ¯¯ 1
0.4 0.35 0.2 0.28
0.7 0.6 0.5
¯ ¯ ¯ ¯ = 0.2706 ¯ ¯
0.2 0.3 0.28
0.6 0.1 0.5
¯ ¯ ¯ ¯ = 0.3529 ¯ ¯
Note that the weights add up to 1. Using function values given in the table, the interpolated value of the function at x = 0.28, y = 0.5 is obtained as f (0.28, 0.5) = 0.3765 × 2.2255 + 0.2706 × 1.4918 + 0.3529 × 2.8577 = 2.2501 Quadratic interpolation For quadratic interpolation, the data at the midpoints of the sides of the triangle element are also used. The weights already calculated for linear interpolation may be used to obtain the solution in the present case. We interpret the weights calculated there in terms of local coordinate system introduced in the text. Hence the coordinates of the given point x = 0.28, y = 0.5 may be identified in terms of the local coordinates as η 1 = 0.3765, η 2 = 0.2706, and η 3 = 0.3529. The following program has been written to perform quadratic interpolation
p1 = [0.2 0.6]; % v e r t i c e s of the t r i a n g l e ( x , y ) p2 = [0.3 0.1]; p3 = [0.35 0.7]; f = [2.2255 1.4918 2.8577 1.8221 2.0647 2.5219]; % f ( x , y ) x = 0.28; y = 0.5; [a ,b , c ] = trlocalcoeff ( p1 , p2 , p3 ) ; % calculate coefficients eta = trianglelocalcord (a ,b ,c ,x , y ) ; % c a l c u l a t e η w (1) = 2* eta (1) *( eta (1) -0.5) ; % calculate weights w (2) = 2* eta (2) *( eta (2) -0.5) ; w (3) = 2* eta (3) *( eta (3) -0.5) ; w (4) = 4* eta (1) * eta (2) ; w (5) = 4* eta (2) * eta (3) ; w (6) = 4* eta (3) * eta (1) ; flin = sum ( eta * f (1:3) ) ; % linear interpolation fquad = sum ( w .* f ) ; % quadratic interpolation
272
Chapter.6 Interpolation in two and three dimensions
The weights for all the points for quadratic interpolation are tabulated below. η
Point No. 1 4 2 5 3 6
0.3765 ··· 0.2706 ··· 0.3529 ···
Weight formula = 2η 1 (η 1 − 0.5) = 4η 1 η 2 = 2η 2 (η 2 − 0.5) = 4η 2 η 3 = 2η 3 (η 3 − 0.5) = 4η 3 η 1
Weight -0.0930 -0.1242 -0.1038 0.4075 0.3820 0.5315
Using the weights and the corresponding function values, we get Point No. 1 4 2 5 3 6 Sum:
w -0.0930 -0.1242 -0.1038 0.4075 0.3820 0.5315 1
f 2.2255 1.4918 2.8577 1.8221 2.0647 2.5219 ···
w× f -0.2070 -0.1852 -0.2966 0.7425 0.7887 1.3403 2.1827
Bold entry in the last column is the desired value of the function at x = 0.28, y = 0.5. This is an improved value as compared to that obtained by linear interpolation. The function is known to be f (x, y) = e(x+ y) and has an exact value of 2.1815 at x = 0.28, y = 0.5. The error for linear interpolation is 0.0686 while it is only 0.0012 in the case of quadratic. There is a huge improvement! We have discussed interpolation polynomials for two dimensional geometries commonly used in FEM and FVM. However, there are many more interpolation techniques like Shepard interpolation, Wachspress interpolation of irregular polynomials, Bsplines and NURBS. Interested readers are advised to read the cited literature.
6.3 Interpolation in three dimensions The concepts from two dimensional interpolation can be easily extended to three and higher dimensions. The most common three dimensional shapes are tetrahedral ( extension of triangle) and hexahedral (extension of quadrilateral). Interpolation formulas can be represented as product of ‘planes’ in place of lines used in the two dimensional cases.
6.3.1
Hexahedral element
Hexahedral element is a solid volume with six faces making it analogous to a quadrilateral in two dimensions. The method of deriving the weight functions is the same, but with product of lines in two dimensions becoming product of planes in the case of three dimensions. Consider the cube in Figure 6.154 . 4
The cube is the local coordinate representation of an arbitrary hexahedral.
6.3. Interpolation in three dimensions 3(−1, 1, −1)
273 4(1, 1, −1)
η 8(1, 1, 1)
7(−1, 1, 1)
ξ
Figure 6.15: Hexahedral element
ζ 1(−1, −1, −1)
5(−1, −1, 1)
2(1, −1, −1)
6(1, −1, 1)
The weight functions for all the points have been given below (1 − ξ)(1 − η)(1 − ζ) 8 (1 − ξ)(1 + η)(1 − ζ) w3 = 8 (1 − ξ)(1 − η)(1 + ζ) w5 = 8 (1 − ξ)(1 + η)(1 + ζ) w7 = 8 w1 =
6.3.2
(1 + ξ)(1 − η)(1 − ζ) 8 (1 + ξ)(1 + η)(1 − ζ) w4 = 8 (1 + ξ)(1 − η)(1 + ζ) w6 = 8 (1 + ξ)(1 + η)(1 + ζ) w8 = 8 w2 =
(6.22)
Tetrahedral element η2
1(x1 , y1 , z1 )
(0, 1, 0) 2(x2 , y2 , z2 ) (x, y, z)
=⇒
η1
(1, 0, 0) 3(x3 , y3 , z3 )
4(x4 , y4 , z4 )
η3
(0, 0, 1)
Figure 6.16: Tetrahedral element and its transformation to local coordinates
A tetrahedral element is a volume with four faces and is analogous to a triangle in two dimensions. The derivation of weight functions for the volume element is similar to the one for triangles. Planes forming the volume are analogous to the lines forming the triangle. Similar to the case of a triangle, the tetrahedron can be transformed to local coordinates (η 1 , η 2 , η 3 , η 4 ). The transformation maps the actual tetrahedral unto right angled tetrahedral as shown in Figure 6.16 (similar to transformation of triangular coordinates) η 1 is the ratio of the distance of the point and distance of node 1 from plane opposite to the node i.e. plane containing nodes 2, 3 and 4. It is also the ratio of the volume of the tetrahedron made by the point with the plane containing nodes 2, 3 and 4 (shaded volume in Figure 6.17) to the volume of the tetrahedral element. It can be easily inferred that for a linear tetrahedral element w i = η i .
274
Chapter.6 Interpolation in two and three dimensions 1(1, 0, 0, 0) Figure
2(0, 1, 0, 0)
to
local
6.17:
Transformation
coordinates:
ratio
of
volume of the shaded tetrahedral to volume of tetrahedral element
3(0, 0, 1, 0)
4(0, 0, 0, 1)
Concluding remarks In this chapter, we have introduced, interpolation in two and three dimensions. These are very important since they are applied in all programs written for solving field problems in engineering. In fact the developments in interpolation techniques have paralleled the developments in FVM, FEM and other commercial software.
Interpolation in two and three dimensions Multidimensional interpolation is commonly encountered in numerical methods such as the Finite Element Method (FEM) and the Finite Volume Method (FVM) used for solving partial differential equations. It is a general practice in numerical methods to discretize a two (three) dimensional domain into large number of small areas (volumes) known as elements in FEM and volumes in FVM. These methods assume a functional form for variable within each sub domain based on the nodal values and solve for the variables at the nodes. Development of multidimensional interpolation has also applications in computer graphics where surfaces are represented by Bezier curves and NURBS (Non-uniform rational B-splines). The present chapter intends to introduce the readers to multi-dimensional interpolation, with simple examples, to be made use of in later chapters.
255
256
6.1 6.1.1
Chapter.6 Interpolation in two and three dimensions
Interpolation over a rectangle Linear interpolation
The simplest interpolation technique is to use a linear variation along the two directions, say x and y for a function f (x, y). An engineering application will be in using table of properties of a pure substance such as water (referred to as Steam Tables). A part of such a table would look like this:
x1 x2
y1 f1 f2
y2 f3 f4
where the columns correspond to constant y’s, rows correspond to constant x’s and the entries correspond to the function f (x, y). We assume that the function varies linearly with both x and y. This means that, for a fixed y, the function varies linearly with x (along a column in the table). Suppose we desire to have the interpolated value of the tabulated function at (x i , yi ) within the range x1 < x i < x2 , y1 < yi < y2 defining a rectangular region (Figure 6.1) in the x, y plane.
f4 f1
fi
f3 y = y2 f2 y = yi y = y1
x=
x1
x=
xi
x=
x2
We hold y fixed at y1 and obtain the interpolated value of the function at x i , y1 by linear interpolation as (x i − x1 ) (x i − x2 ) f1 + f 2 = w1x f 1 + w2x f 2 f (x i , y1 ) = (x1 − x2 ) (x2 − x1 ) based on L 1 , Lagrange polynomial of degree 1. The weights have been identified with suitable subscripts. Similarly we hold y fixed at y2 and obtain the interpolated value of the function at x i , y2 by linear interpolation as f (x i , y2 ) =
(x i − x2 ) (x i − x1 ) f3 + f 4 = w1x f 3 + w2x f 4 (x1 − x2 ) (x2 − x1 )
Now we hold x fixed at x i and obtain the interpolated value of function at x i , yi as f (x i , yi )
= =
(yi − y2 ) (yi − y1 ) f (x i , y1 ) + f (x i , y2 ) (y1 − y2 ) (y2 − y1 ) w1y f (x i , y1 ) + w2y f (x i , y2 )
(6.1)
6.1. Interpolation over a rectangle
257
where again the weights have been identified by suitable subscripts. Specifically, x in the subscript means the weight is a function of x while y in the subscript means the weight is a function of y. Substituting Equations 6.1 and 6.1 in Equations 6.1 we get f (x i , yi )
=
w1x w1y f 1 + w2x w1y f 2 + w1x w2y f 3 + w2x w2y f 4
=
w1 f 1 + w2 f 2 + w3 f 3 + w4 f 4
(6.2)
where the final weights w1 - w4 are functions of both x and y. For example, we have w1 = w1x w1y = 1
(x i − x2 ) (yi − y2 ) (x1 − x2 ) (y1 − y2 )
(6.3)
1
w1x
w2x
0
0
3
4
y = y2
0
w2y
1 Figure 6.2: Graphical representa-
tion of weight functions
x = x2
x = x1
(x i , yi )
y = y1
1
2
0
w1y
1
The reader may similarly write down the other weights also. Thus the interpolated value is a weighted sum of the function at the 4 nodes. We see that the weights are products of Lagrange weights in one dimension. We also see that the weight is unity at (x1 , y1 ) and vanishes for x i = x2 or yi = y2 which means that this weight vanishes at the three nodes given by (x2 , y1 ), (x1 , y2 ) and (x2 , y2 ). The reader may verify, in a similar fashion, the properties of the other three weights also. Figure 6.2 illustrates the construction of weight or shape functions over a rectangular domain.
Example 6.1 Following data has been taken from a table of property values of a substance, after suppressing information on the units. x↓ y→ 15 25
500 1286 1285
600 1334.2 1333.6
Obtain the value of the function at x = 20, y = 530 using linear interpolation.
Solution : The weights are calculated with x i = 20, yi = 530 as given below. w1 =
(20 − 25)(530 − 600) = 0.35 (15 − 25)(500 − 600
w2 =
(20 − 15)(530 − 600) = 0.35 (25 − 20)(500 − 600
258
Chapter.6 Interpolation in two and three dimensions
w3 =
(20 − 25)(530 − 500) = 0.15 (15 − 25)(600 − 500
w4 =
(20 − 15)(530 − 500) = 0.15 (25 − 15)(600 − 500
The interpolated value of the function is then obtained as f (20, 530) = 0.35 × 1286 + 0.35 × 1285 + 0.15 × 1334.2 + 0.15 × 1333.6 = 1300.02
Linear interpolation; alternate way of representation Consider the weight w1 given in Equation 6.3. We may rewrite it in full as x i yi − x2 yi − y2 x i + x2 y2 w1 = = a 10 + a 11 x i + a 12 yi + a 13 x i yi x1 y1 − x2 y1 − y2 x1 + x2 y2 where x2 y2 x1 y1 − x2 y1 − y2 x1 + x2 y2 x2 a 12 = − x1 y1 − x2 y1 − y2 x1 + x2 y2 a 10 =
y2 x1 y1 − x2 y1 − y2 x1 + x2 y2 1 a 13 = x1 y1 − x2 y1 − y2 x1 + x2 y2 a 11 = −
(6.4)
It is easily shown all the other weights have the same form but with different coefficients (the a’s). It is thus seen that the interpolating function may be written in the alternate form f (x i , yi ) = c 0 + c 1 x i + c 2 yi + c 3 x i yi where the c’s are constants determined in following equation set to obtain the c’s. 1 x1 y1 1 x2 y1 1 x1 y2 1 x2 y2
(6.5)
terms of the four function values. In fact, we have the x1 y1 x2 y1 x1 y2 x2 y2
c0 c1 c 2 c3
f1 f2 = f 3 f4
(6.6)
These equations may easily be solved to get the vector c of coefficients in terms of the function values at the corners of the rectangle.
6.1.2
Local coordinate system for a rectangular element
A simple linear transformation may be used to transform the rectangular domain in x, y to a square of sides =2 in ξ, η such that −1 < ξ, η < 1 with the origin at the center of the square (see Figure 6.3). The required transformation may easily shown to be x1 + x2 x2 − x1 y1 + y2 y2 − y1 x= + ξ; y = + η (6.7) 2 2 2 2 Instead of interpolating in (x, y) we may interpolate in (ξ, η).
3(x1 , y2 )
30 (−1, 1)
4(x2 , y2 )
η
40 (1, 1) (ξ, η)
(x i , yi )
1(x1 , y1 )
ξ
=⇒
2(x2 , y1 )
10 (−1, −1)
Figure 6.3: Coordinate transformation
20 (1, −1)
6.1. Interpolation over a rectangle
259
Such a transformation is useful in engineering applications such as the finite element method. The given data may now be recast as follows: ξ = −1 ξ=1
η = −1 f1 f2
η=1 f3 f4
Equation set 6.6 is recast as 1 1 1 1
−1 1 −1 1
−1 −1 1 1
1 −1 −1 1
c0 c1 c2 c3
f1 f2 = f 3 f4
(6.8)
These equations may be solved easily (eg. Gauss elimination method) to get − f1 + f2 − f3 + f4 f1 + f2 + f3 + f4 c1 = 4 4 − f1 − f2 + f3 + f4 f1 − f2 − f3 + f4 c2 = c3 = (6.9) 4 4 We would now like to write the interpolating function in terms of weighted sum as given earlier. 1 − ξ − η + ξη We easily see that w1 (ξ, η) will be given by which may be recast in the form 4 (1 − ξ) (1 − η) w1 (ξ, η) = (6.10) 2 2 It is deliberately written in the form of a product of two Lagrange polynomials of first degree. We may show similarly that the following hold: (1 + ξ) (1 − η) (1 − ξ) (1 + η) w2 (ξ, η) = w3 (ξ, η) = 2 2 2 2 (6.11) (1 + ξ) (1 + η) w4 (ξ, η) = 2 2 Of course we may revert to (x, y) once the interpolation has been worked out in (ξ, η).
c0 =
Coordinate transformation
The advantage of coordinate transformations is that operations such as differentiation, integration and transforms can be applied on a local coordinate system rather than the actual coordinates. The mapping of any arbitrary quadrilateral to the local coordinate system can be achieved in a number of ways such as affine, projective transformations. The discussion on these methods are beyond the scope of the book and the interested reader should look into books on computational graphics for the same. Isoparametric transformations are most commonly employed in numerical methods such as FEM and FVM. The transformation maps line to line such that both variables ( f ) and the real coordinates (x, y) are evaluated using the same shape functions. For an arbitrary quadrilateral considered below, the coordinate transformation can be simply written as ξ 4 30 (−1, 1) l 2 40 (1, 1) l2 3 (x, y) l1 (η, ξ) l3 l1 l3 =⇒ 2 η l4 1
10 (−1, −1) l 4
20 (1, −1)
260
Chapter.6 Interpolation in two and three dimensions
x = w1 x1 + w2 x2 + w3 x3 + w4 x4 ;
y = w1 y1 + w2 y2 + w3 y3 + w4 y4
where w’s are the weight functions (evaluated in (ξ, η) coordinates) we are familiar with. If we closely observe, the coordinate transformation equation is of the form x = a 0 + a 1 ξ + a 2 η + a 3 ξη
y = b 0 + b 1 ξ + b 2 η + b 3 ξη
(6.12)
The two dimensional domain need not be represented by straight lines but could be any four arbitrary curves in which case coordinate transformation becomes murkier.
Example 6.2 Redo example 6.1 by transforming the coordinates and using the formulation given above
Solution : The transformations are easily obtained as follows. With x2 = 25 and x1 = 15, we have x1 + x2 15 + 25 x2 − x1 25 − 15 x − 20 = = 20; = = 5. Hence x = 5ξ + 20 or ξ = . Similarly we 2 2 2 2 5 y − 550 may write the second transformation as y = 50η + 550 or η = . Thus the interpolation 50 20 − 20 530 − 550 point corresponds to ξ = = 0 and η = = −0.4. The weights are then 2 50 calculated using Equations 6.10 and 6.11. w1
=
w2
=
w3
=
w4
=
(1 − 0)[1 − (−0.4)] 4 (1 + 0)[1 − (−0.4)] 4 (1 − 0)[1 + (−0.4)] 4 (1 + 0)[1 + (−0.4)] 4
= 0.35 = 0.35 = 0.15 = 0.15
The interpolated value is the same as f (0, 0.4) = 1300.02 obtained earlier in Example 6.1. Program 6.1: Program to calculate weights of a linear quadrilateral element
1 2 3 4 5 6 7 8 9
6.1.3
function w = weightquadlin (x , y ) % Input % Output
x,y : w :
l o c a l c o o r d i n a t e s of the point weights at point x , y n = size ( x ) ; % no o f p o i n t s w = zeros (n ,4) ; % i n i t i a l i z e weight w (: ,1) = (1 - x ) .*(1 - y ) /4; % c a l c u l a t e w e i g h t w (: ,2) = (1+ x ) .*(1 - y ) /4; % c a l c u l a t e w e i g h t w (: ,3) = (1 - x ) .*(1+ y ) /4; % c a l c u l a t e w e i g h t w (: ,4) = (1+ x ) .*(1+ y ) /4; % c a l c u l a t e w e i g h t
(s)
at at at at
( −1 , −1) (1 , −1) ( −1 ,1) (1 ,1)
Interpolating polynomials as products of ‘lines’
Now we are ready for another interpretation of the interpolating function we have obtained above. Consider Equation 6.10 as an example. It contains two factors, the first one containing 1−ξ the equation of a line l 1 (ξ) = which is a line lying in a plane that is normal to the ξ, η plane 2
6.1. Interpolation over a rectangle
261
that contains the ξ axis. This line passes through the points l 1 (−1) = 1 and l 1 (1) = 0. Similarly 1−η the second factor represents the equation of a line given by l 2 (η) = which is a line lying in a 2 plane that is normal to the ξ, η plane that contains the η axis. This line passes through the points l 2 (−1) = 1 and l 2 (1) = 0. Thus the weight is the product of equations of lines l 1 and l 2 . Similar interpretation is possible for the other weights also. What is important to notice is that the weight is a product of two lines. These lines represent Lagrange polynomials of first degree (see Figure 6.4).
wi l 1 or l 2
1
l 3 or l 4 Figure 6.4: Four `lines' used in linear
h = d/2
interpolation over a standard square
d
0
−1
ξ or η
−1
Using the ‘lines’ defined above we may recast the weights as w1 = l 1 l 2 ; w2 = l 2 l 3 ; w3 = l 1 l 4 and w4 = l 3 l 4
(6.13)
As indicated in the figure the height of the ‘line’ is equal to half the distance from the farthest point. Hence we may use the side of the square as a proxy for the corresponding line by simply halving d itself (i.e. by dividing the distance by the total distance, which is 2). This is shown schematically in Figure 6.5. The nodes are numbered as indicated and the lines are identified with the sides of the square domain. The weights are shown at the corresponding nodes as product of ‘lines’. This way of visualizing the weights in terms of lines is particularly useful in dealing with higher order interpolation and also interpolation over a triangular domain.
1−η 2
l1 =
l3 =
1+ξ 2
(ξ, η)
w1 = l 1 l 2 l = 1 + η 4 2
6.1.4
w4 = l 3 l 4
1−ξ 2
w3 = l 1 l 4 l 2 =
Figure 6.5:
Rectangular element
in local coordinates showing the lines
w2 = l 2 l 3
Lagrange Quadratic rectangular element
Figure 6.6 shows a 9 noded quadratic rectangular element also known as Lagrange quadratic element. The interpolating function for the element will be of the form f (ξ, η) = a 0 + a 1 ξ + a 2 η + a 3 ξη + a 4 ξ2 + a 5 η2 + a 6 ξη2 + a 7 ξ2 η + a 8 ξ2 η2
(6.14)
262
Chapter.6 Interpolation in two and three dimensions η
(−1, 1) 7
l2
(0, 1) 8
l3
(ξ, η)
(−1, 0) 4
l5
5 (0, 0)
(1, 1) 9 l1 6 (1, 0)
Figure 6.6:
ξ
Lagrange quadratic
rectangular element
l6 1 (−1, −1)
l4
2 (0, −1)
3 (1, −1)
The constants a i have to be determined such that f (ξ, η) is the same as the given data at all the nine nodes. Alternatively, the interpolation function can be estimated using Lagrange approach such that 9 X f (ξ, η) = w i (ξ, η) f i i =1
We can define two independent quadratic Lagrange polynomial functions along η and ξ directions. The total weight functions, following the procedure used in the case of linear interpolation, is the product of the two quadratic polynomial functions. The weight functions are of the form w(ξ, η) = (a 0 + a 1 ξ + a 2 ξ2 )(b 0 + b 1 η + b 2 η2 ) = F(ξ)G(η) Let us decompose the weight function and determine the weights of two separate functions F(ξ) and G(η). Then the weight functions at the nine points would be w1 = F−1 G −1 w4 = F−1 G 0 w7 = F−1 G 1
w 2 = F 0 G −1 w5 = F 0 G 0 w8 = F 0 G 1
w3 = F1 G −1 w6 = F 1 G 0 w9 = F 1 G 1
Based on Equation 5.17, the Lagrange weights L 2 (ξ) and L 2 (η) are written in vector form as FT (ξ)
=
GT (η)
=
ξ(ξ − 1) ξ(ξ + 1) ξ(ξ − 1) ξ(ξ + 1) = − (ξ − 1)(ξ + 1) (1 − ξ2 ) 2 2 2 2 η(η − 1) η(η + 1) η(η − 1) η(η + 1) 2 − (η − 1)(η + 1) = (1 − η ) 2 2 2 2
The nodal weights may then be written down as products GFT and are given as
=
w1 w2 w3 w w5 w6 = GFT 4 w7 w8 w9 2 ξη(ξ − 1)(η − 1) (1 − ξ )η(η − 1) 4 2 2 ξ(ξ − 1)(1 − η ) (1 − ξ2 )(1 − η2 ) 2 2 ξη(ξ − 1)(η + 1) (1 − ξ )η(η + 1) 4 2
ξη(ξ + 1)(η − 1) 4 2 ξ(ξ + 1)(1 − η ) 2 ξη(ξ + 1)(η + 1) 4
(6.15)
6.1. Interpolation over a rectangle
263
Alternate approach: product of lines The intersection points of the lines l 1 to l 6 are the nodes of the square domain. The nodes can be represented in terms of these lines. The equation of these lines are l 1 =⇒ 1 − ξ = 0;
l 2 =⇒ 1 − η = 0;
l 4 =⇒ 1 + η = 0;
l 5 =⇒ η = 0;
l 3 =⇒ 1 + ξ = 0; l 6 =⇒ ξ = 0
Node 1: The lines passing through node 1 are l 3 and l 4 . At node 1, all the weights except w1 are zero. In order to achieve this, all the weight functions except w1 should have either l 3 or l 4 as factors. Similarly, node 1 would be a factor of remaining lines l 1 , l 2 , l 5 and l 6 . Therefore the weight function for node 1 can be written down as w1 =
ξη(ξ − 1)(η − 1) l1 l2 l5 l6 = l 1 (−1, −1)l 2 (−1, −1)l 5 (−1, −1)l 6 (−1, −1) 4
which is the same as that obtained earlier. Similar approach can be adopted to derive weight functions for the other nodes.
Example 6.3 Obtain value of the tabulated function at x = 0.14, y = 0.07 by quadratic interpolation. x↓ y→ 0 0.1 0.2
0 1.000000 1.010050 1.040811
0.1 1.010050 1.020201 1.051271
0.2 1.040811 1.051271 1.083287
Use a nine noded element to solve the problem.
Solution : We transform the domain to the standard square domain by the transformation ξ = 10(x − 0.1), η = 10(y − 0.1). The elements are numbered as shown in Figure 6.6. In the present example, the function is required at x = 0.14 which corresponds to ξ = 10(0.14 − 0.1) = 0.4 and y = 0.07 which corresponds to η = 10(0.07 − 0.1) = −0.3. Substituting these in the matrix of weights (Equation 6.15) we get −0.0234 0.1638 0.0546 −0.1092 0.7644 0.2548 w= 0.0126 −0.0882 −0.0294 The interpolated value of the function is then obtained as sum of product of weights and the function values. We thus obtain f (ξ = 0.4, η = −0.3) = f (x = 0.14, y = 0.07) =
9 X
w i f i = 1.024826
i =1
The reader may verify that linear interpolation using only the corner nodes would give a 2 2 value for this as 1.043259. The function used for generating the data was f (x, y) = e(x + y ) which has a value of f (0.14, 0.07) = 1.024803. The quadratic interpolated value has an error of 0.000024 while the linear interpolated value has an error of 0.018457 ≈ 0.019, with respect to the exact value. There is a considerable improvement in using quadratic interpolation over linear interpolation in the present case.
264
Chapter.6 Interpolation in two and three dimensions η
6(−1, 1)
7(0, 1)
l2
8(1, 1) (ξ, η)
l3
l6
l5
l1 5(1, 0) ξ
(−1, 0)4
Figure 6.7:
Eight noded rectan-
gular element
l7
l4
1(−1, −1)
6.1.5
l8
2(0, −1)
3(1, −1)
Quadratic eight noded rectangular element
The above element is also known as ‘serendipity’ quadratic element.1 The interpolation function in the domain is of the form f (ξ, η) = a 0 + a 1 ξ + a 2 η + a 3 ξη + a 4 ξ2 + a 5 η2 + a 6 ξη2 + a 7 ξ2 η
(6.16)
The present domain does not have a point at the origin. Therefore the present interpolation formula does not support the term involving ξ2 η2 . Unlike the previous example, the weight function of the present interpolation formula cannot be simplified as product of two univariate functions in ξ and η. Let us adopt the Lagrange interpolation approach based on product of lines. The nodes of the domain are the points of intersection of lines l 1 to l 8 . Hence these lines can be used to represent the nodes. Equations of lines indicated in Figure 6.7 are given as a table below. Equations of lines: Line l1 l2 l3 l4
Formula 1−ξ 1−η 1+ξ 1+η
Line l5 l6 l7 l8
Formula ξ+η−1 −ξ + η − 1 −ξ − η − 1 ξ−η−1
Calculation of nodal weights as product of lines: For node 1 the procedure is as follows: l 1 and l 2 must be factors of w1 as these will be zero at nodes 3,5,6,7 and 8. The function should be equal to zero at nodes 2 and 4 as well. This means l 7 is also a factor of the weight function. Therefore the weight function can be written as w1 =
l1 l2 l7 (1 − ξ)(1 − η)(−ξ − η − 1) = l 1 (−1, −1)l 2 (−1, −1)l 7 (−1, −1) 4
Similarly we treat other nodes to arrive at weight functions as given below. 1
For a serendipity element function values are specified only at points on the boundary.
6.1. Interpolation over a rectangle
Node Number 1 2 3 4 5 6 7 8
265
Weight formula in terms of product of lines (l 1 × l 2 × l 7 )/4 (l 1 × l 2 × l 3 )/2 (l 2 × l 3 × l 8 )/4 (l 1 × l 2 × l 4 )/2 (l 2 × l 3 × l 4 )/2 (l 1 × l 4 × l 6 )/4 (l 1 × l 3 × l 4 )/2 (l 3 × l 4 × l 5 )/2
Weight expressions (1 − ξ)(1 − η)(−ξ − η − 1)/4 (1 − ξ)(1 − η)(1 + ξ)/2 (1 − η)(1 + ξ)(ξ − η − 1)/4 (1 − ξ)(1 − η)(1 + η)/2 (1 − η)(1 + ξ)(1 + η)/2 (1 − ξ)(1 + η)(−ξ + η − 1)/4 (1 − ξ)(1 + ξ)(1 + η)/2 (1 + ξ)(1 + η)(ξ + η − 1)4
Note: One may be tempted to use lines η = 0 and ξ = 0 instead of l 5 to l 8 . However, such an interpolation polynomial would not be consistent with Equation 6.16.
Example 6.4 A certain function is available at eight points as given in the following table. Obtain value of the function at x = 0.14, y = 0.07 by quadratic interpolation. x↓ y→ 0 0.1 0.2
0 1.000000 1.010050 1.040811
0.1 1.010050 ··· 1.051271
0.2 1.040811 1.051271 1.083287
Solution : Note that this set of data is the same as that used in Example 6.3 obtained by removing the data at the origin (removing data at point 5 in Figure 6.6). Nodes are numbered as shown in Figure 6.7. In the present example interpolated value is desired at ξ = 0.4, η = −0.3. We obtain all the weights using the weight functions and tabulate the results. Node No. 1 2 3 4 5 6 7 8 Sum =
f (x, y) 1.000000 1.010050 1.040811 1.010050 1.051271 1.040811 1.051271 1.083287 ···
W ei ght -0.2145 0.546 -0.1365 0.273 0.637 -0.1785 0.294 -0.2205 1
W ei ght × f (x, y) -0.214500 0.551487 -0.142071 0.275744 0.669660 -0.185785 0.309074 -0.238865 1.024744
Bold entry in the table is the required value of the function at the interpolation point.
266
6.2
Chapter.6 Interpolation in two and three dimensions
Interpolation over a triangle
Interpolation over a triangle is an important exercise since most numerical methods for solving partial differential equations use triangular elements. Consider a triangular element shown in Figure 6.8.
3 (x3 , y3 )
l2
(x, y)
l1 Figure 6.8: Linear triangular element
1 (x1 , y1 )
l3
2 (x2 , y2 )
l 1 , l 2 and l 3 are the equations of the lines forming the edges of the triangle. As there are three points, the interpolating function is of the form2 f (x, y) = a 0 + a 1 x + a 2 y
(6.17)
The constants a 0 , a 1 and a 2 can be determined by solving the system of equations at the three nodes of the triangle. Alternatively, we can apply the Lagrange interpolation scheme for the 3 X triangle such that f (x, y) = wi f i . i =1
Consider node 1. The weight w1 should be equal to 1 at node 1 and should be equal to zero at other two nodes. l 2 and l 3 are equal to zero at node 1. Therefore the weight function w1 is proportional to l 1 and is of the form given by Equation 6.17. Meeting all the required conditions, w1 can be written as l 1 (x, y) w1 = l 1 (x1 , y1 ) In general, the weight function at node i can be written as wi =
l i (x, y) l i (x i , yi )
(6.18)
Local coordinate system for a triangular element l i (x, y) = η i . η i is nothing but the ratio of normal distance of point (x, y) from l i (x i , yi ) line l i to the normal distance of node i from the line. When the point lies on l i , η i = 0 and when the point is at node i, η i = 1. Hence, the following representation also serves as an alternative coordinate system for the triangular element with η i as the basis. Any point in (x, y) coordinates can be represented in terms of (η 1 , η 2 , η 3 ) as indicated in the following figure. Let the fraction
2
Number of nodes less than four would mean that there can be only three terms in the interpolating function. The x y term would be missing.
6.2. Interpolation over a triangle
267
3 (x3 , y3 )
3 (0, 0, 1)
(η 1 , η 2 , η 3 )
η2
(x, y)
η1 η3
1 (x1 , y1 )
2 (x2 , y2 )
1 (1, 0, 0)
2 (0, 1, 0)
Figure 6.9: Transformation of coordinates
Hence the equation of the three lines l 1 , l 2 and l 3 in local coordinate system would be η 1 = 0, η 2 = 0 and η 3 = 0. For a linear triangular element, w i = η i The area of a triangle is proportional to product of height and base. Hence η i is also the ratio of the area of the triangles made by points (x, y) and (x i , yi ) with line l i . This can be understood from Figure 6.10.
3 (x3 , y3 )
l2
l1 Figure
(x, y)
6.10:
Calculation
of
weight functions for a triangular element
1 (x1 , y1 )
2 (x2 , y2 )
l3
η 3 is the ratio of the area of the shaded triangle to the area of the triangle itself. The weight n X function for a linear triangular element can be generalized as w i = η i = A i /A. Since A i = A, i =1
sum of weights at any point is equal to 1 i.e. given by the following determinant ¯ ¯ 1 1 ¯¯ A3 = ¯ 1 2¯ 1
x1 x2 x
y1 y2 y
3 X
w i (x, y) = 1. The area of the shaded triangle is
i =1
¯ ¯ ¯ 1 ¯ = { x(y1 − y2 ) + (x2 − x1 )y + (x1 y2 − x2 y1 )} ¯ 2 ¯
When (x, y) lies on the line l 3 , A 3 = 0 and when (x, y) = (x3 , y3 ), A 3 = A, where A is the area of the triangular element given by ¯ ¯ ¯ 1 x1 y1 ¯ ¯ 1 ¯¯ A = ¯ 1 x2 y2 ¯¯ 2¯ 1 x3 y3 ¯
268
Chapter.6 Interpolation in two and three dimensions
η i is of the form (a i + b i x + c i y) where a i ,b i and c i are constants and are evaluated as follows
(x2 y3 − y2 x3 ) 2A (x3 y1 − y3 x1 ) a2 = 2A (x1 y2 − y1 x2 ) a3 = 2A a1 =
(y2 − y3 ) 2A (y3 − y1 ) b2 = 2A (y1 − y2 ) b3 = 2A b1 =
(x3 − x2 ) 2A (x1 − x3 ) c1 = 2A (x2 − x1 ) c1 = 2A c1 =
(6.19)
Interpretation using local coordinates: Any local coordinate can be represented by other two local coordinates i.e. η 3 = 1 − η 1 − η 2 . The right angled triangle on η 1 ,η 2 plane in Figure 6.11 represents the map of the actual triangle. In a sense, the transformation exercise can be looked upon as mapping triangle onto a right angled isosceles triangle. MATLAB program has been given
η2
2(0, 1) η2 = 1 − η1
Figure
6.11:
Representation of local
coordinates of a triangle
η1
3(0, 0)
1(1, 0)
below to calculate the coefficients for a triangular element Program 6.2: Calculate coe cients for local coordinates of a triangular element
1 function [a ,b ,c , A ] = trlocalcoeff ( p1 , p2 , p3 ) 2 % I n p u t p1 , p2 , p 3 : c o o r d i n a t e s o f t h e v e r t i c e s 3 4 5 6 7 8 9 10 11 12 13 14 15 16
p1 ( 1 ) = x1 ; p1 ( 2 ) % = y1 % p2 ( 1 ) = x2 ; p2 ( 2 ) = y2 ; p3 ( 1 ) = x3 ; p3 ( 2 ) = y3 % Output a,b, c : c o e f f i c i e n t s η = a + bx + c y % A : area of t r i a n g l e a (1) = p2 (1) * p3 (2) - p3 (1) * p2 (2) ; % c a l c u l a t e a
a (2) = p3 (1) * p1 (2) - p1 (1) * p3 (2) ; a (3) = p1 (1) * p2 (2) - p2 (1) * p1 (2) ; A = ( a (1) + a (2) + a (3) ) /2; a = a /(2* A ) ; b (1) = ( p2 (2) - p3 (2) ) /(2* A ) ; b (2) = ( p3 (2) - p1 (2) ) /(2* A ) ; b (3) = ( p1 (2) - p2 (2) ) /(2* A ) ; c (1) = ( p3 (1) - p2 (1) ) /(2* A ) ; c (2) = ( p1 (1) - p3 (1) ) /(2* A ) ; c (3) = ( p2 (1) - p1 (1) ) /(2* A ) ;
% A = 2* area o f % normailze a % calculate b
% calculate
triangle
c
The following Matlab function can be used to convert from actual coordinates to local coordinates. Program 6.3: Convert from
(x, y)
to local coordinates
1 function eta = trianglelocalcord (a ,b ,c ,x , y ) 2 % Input a , b , c : c o e f f i c i e n t s η = a + bx + c y 3 % x , y : coordinates of points
(η 1 , η 2 , η 3 )
of a triangle
6.2. Interpolation over a triangle
4 5 6 7 8 9
% Output
eta
:
269
local
coordinates
of
points
x,y
n = size (x ,1) ; eta = zeros (n ,3) ; for i =1:3 eta (: , i ) = a ( i ) + b( i ) * x + c ( i ) * y ; end
Quadratic triangular element Now we derive expressions for weight functions for a quadratic triangular element.
3
(x, y)
5
4
Figure 6.12:
Quadratic triangular ele-
ment
6
1
2
There are six points. Three nodes occupy the vertices of the triangle and three points are located at the midpoints of the three edges of the triangular element. The interpolation function would be of the form3 f (x, y) = a 0 + a 1 x + a 2 y + a 3 x y + a 4 x2 + a 5 y2 (6.20) l 1 , l 2 and l 3 are the equations of the lines forming the edges of the triangle. l 4 , l 5 and l 6 are the
3
l6
6
5
l2
l3
6.13:
Lines
representing
quadratic triangular element
l1 l5
l4 1
Figure
4
2
equation of the lines joining the midpoints of the edges of the triangle. Let us represent the nodes in terms of the triangular coordinate system introduced earlier. 1 =⇒ (1, 0, 0) 4 =⇒ (0.5, 0.5, 0)
2 =⇒ (0, 1, 0) 5 =⇒ (0, 0.5, 0.5)
3 =⇒ (0, 0, 1) 6 =⇒ (0.5, 0, 0.5)
(6.21)
From Equation 6.21, the equation of the lines can be written down as l 1 =⇒ η 1 = 0 3
l 2 =⇒ η 2 = 0
l 3 =⇒ η 3 = 0
With only six points available the quadratic is constrained to the form shown here. Recall that nine points gave the more general expression (Equation 6.14)
270
Chapter.6 Interpolation in two and three dimensions l 4 =⇒ η 1 − 0.5 = 0
l 5 =⇒ η 2 − 0.5 = 0
l 6 =⇒ η 3 − 0.5 = 0
Now, we derive the weight function for all the nodes. Consider node 1: The weight w1 should be equal to 1 at node 1 and should be equal to zero at the other nodes. Lines l 2 and l 3 are zero where as l 1 = 1 at node 1. Therefore, w1 should be proportional to l 1 . This means w1 = 0 at nodes 2, 3 and 5. As w1 should be zero at nodes 4 and 6, w1 should be proportional to l 4 also. Hence, meeting all the required conditions w1 can be written as l 1 (η 1 , η 2 , η 3 )l 4 (η 1 , η 2 , η 3 ) = 2η 1 (η 1 − 0.5) w1 = l 1 (1, 0, 0)l 4 (1, 0, 0) Consider Node 4: As node 4 is a point on l 3 , η 3 = 0. However, lines l 1 and l 2 are not equal to zero at node 4. This means w4 is proportional to l 1 l 2 . Therefore w4 = 0 at other nodes except node 4. Meeting all the required conditions w4 can be written as w4 =
l 1 (η 1 , η 2 , η 3 )l 2 (η 1 , η 2 , η 3 ) = 4η 1 η 2 l 1 (0.5, 0.5, 0)l 2 (0.5, 0.5, 0)
The weights for all other points may be similarly determined and are summarized below w1 = 2η 1 (η 1 − 0.5); w4 = 4η 1 η 2 ;
w2 = 2η 2 (η 2 − 0.5); w5 = 4η 2 η 3 ;
w3 = 2η 3 (η 3 − 0.5); w6 = 4η 3 η 1
Example 6.5 Function value is available at six points as shown in the table below. Point No. 1 4 2 5 3 6
x 0.2 0.25 0.3 0.325 0.35 0.275
y 0.6 0.35 0.1 0.4 0.7 0.65
f (x, y) 2.2255 1.8221 1.4918 2.0647 2.8577 2.5219
Obtain the value of the function at x = 0.28, y = 0.5 using linear and quadratic interpolation.
Solution : Figure 6.14 shows the node numbering system used in the analysis. Linear interpolation Nodes at the apex of triangle are used for linear interpolation. We evaluate all the 1 weights using the determinants presented in the text. Factor is dropped in all the area 2 calculations since the weights are normalized with the area of the triangle. The area of the triangle is proportional to A given by ¯ ¯ ¯ 1 0.2 0.6 ¯¯ ¯ 0.1 ¯¯ = 0.085 A = ¯¯ 1 0.3 ¯ 1 0.35 0.7 ¯ Weights are now calculated as given below. Weight are given by the following. ¯ ¯ ¯ 1 0.3 0.1 ¯¯ 1 ¯¯ 1 0.35 0.7 ¯¯ = 0.3765 η 1 = w1 = 0.085 ¯¯ 1 0.28 0.5 ¯
6.2. Interpolation over a triangle
271
0.8 3
0.7 0.6
6 1
Figure
(0.28,0.5)
0.5
y
of
for
Numbering linear
and
quadratic interpolation over a triangle showing the inter-
5
0.4
6.14:
nodes
4
polation point
0.3 0.2 2
0.1 0 0.15
0.2
0.25
x
0.3
0.35
¯ ¯ 1 1 ¯¯ 1 η 2 = w2 = 0.085 ¯¯ 1 ¯ ¯ 1 1 ¯¯ 1 η 3 = w3 = 0.085 ¯¯ 1
0.4 0.35 0.2 0.28
0.7 0.6 0.5
¯ ¯ ¯ ¯ = 0.2706 ¯ ¯
0.2 0.3 0.28
0.6 0.1 0.5
¯ ¯ ¯ ¯ = 0.3529 ¯ ¯
Note that the weights add up to 1. Using function values given in the table, the interpolated value of the function at x = 0.28, y = 0.5 is obtained as f (0.28, 0.5) = 0.3765 × 2.2255 + 0.2706 × 1.4918 + 0.3529 × 2.8577 = 2.2501 Quadratic interpolation For quadratic interpolation, the data at the midpoints of the sides of the triangle element are also used. The weights already calculated for linear interpolation may be used to obtain the solution in the present case. We interpret the weights calculated there in terms of local coordinate system introduced in the text. Hence the coordinates of the given point x = 0.28, y = 0.5 may be identified in terms of the local coordinates as η 1 = 0.3765, η 2 = 0.2706, and η 3 = 0.3529. The following program has been written to perform quadratic interpolation
p1 = [0.2 0.6]; % v e r t i c e s of the t r i a n g l e ( x , y ) p2 = [0.3 0.1]; p3 = [0.35 0.7]; f = [2.2255 1.4918 2.8577 1.8221 2.0647 2.5219]; % f ( x , y ) x = 0.28; y = 0.5; [a ,b , c ] = trlocalcoeff ( p1 , p2 , p3 ) ; % calculate coefficients eta = trianglelocalcord (a ,b ,c ,x , y ) ; % c a l c u l a t e η w (1) = 2* eta (1) *( eta (1) -0.5) ; % calculate weights w (2) = 2* eta (2) *( eta (2) -0.5) ; w (3) = 2* eta (3) *( eta (3) -0.5) ; w (4) = 4* eta (1) * eta (2) ; w (5) = 4* eta (2) * eta (3) ; w (6) = 4* eta (3) * eta (1) ; flin = sum ( eta * f (1:3) ) ; % linear interpolation fquad = sum ( w .* f ) ; % quadratic interpolation
272
Chapter.6 Interpolation in two and three dimensions
The weights for all the points for quadratic interpolation are tabulated below. η
Point No. 1 4 2 5 3 6
0.3765 ··· 0.2706 ··· 0.3529 ···
Weight formula = 2η 1 (η 1 − 0.5) = 4η 1 η 2 = 2η 2 (η 2 − 0.5) = 4η 2 η 3 = 2η 3 (η 3 − 0.5) = 4η 3 η 1
Weight -0.0930 -0.1242 -0.1038 0.4075 0.3820 0.5315
Using the weights and the corresponding function values, we get Point No. 1 4 2 5 3 6 Sum:
w -0.0930 -0.1242 -0.1038 0.4075 0.3820 0.5315 1
f 2.2255 1.4918 2.8577 1.8221 2.0647 2.5219 ···
w× f -0.2070 -0.1852 -0.2966 0.7425 0.7887 1.3403 2.1827
Bold entry in the last column is the desired value of the function at x = 0.28, y = 0.5. This is an improved value as compared to that obtained by linear interpolation. The function is known to be f (x, y) = e(x+ y) and has an exact value of 2.1815 at x = 0.28, y = 0.5. The error for linear interpolation is 0.0686 while it is only 0.0012 in the case of quadratic. There is a huge improvement! We have discussed interpolation polynomials for two dimensional geometries commonly used in FEM and FVM. However, there are many more interpolation techniques like Shepard interpolation, Wachspress interpolation of irregular polynomials, Bsplines and NURBS. Interested readers are advised to read the cited literature.
6.3 Interpolation in three dimensions The concepts from two dimensional interpolation can be easily extended to three and higher dimensions. The most common three dimensional shapes are tetrahedral ( extension of triangle) and hexahedral (extension of quadrilateral). Interpolation formulas can be represented as product of ‘planes’ in place of lines used in the two dimensional cases.
6.3.1
Hexahedral element
Hexahedral element is a solid volume with six faces making it analogous to a quadrilateral in two dimensions. The method of deriving the weight functions is the same, but with product of lines in two dimensions becoming product of planes in the case of three dimensions. Consider the cube in Figure 6.154 . 4
The cube is the local coordinate representation of an arbitrary hexahedral.
6.3. Interpolation in three dimensions 3(−1, 1, −1)
273 4(1, 1, −1)
η 8(1, 1, 1)
7(−1, 1, 1)
ξ
Figure 6.15: Hexahedral element
ζ 1(−1, −1, −1)
5(−1, −1, 1)
2(1, −1, −1)
6(1, −1, 1)
The weight functions for all the points have been given below (1 − ξ)(1 − η)(1 − ζ) 8 (1 − ξ)(1 + η)(1 − ζ) w3 = 8 (1 − ξ)(1 − η)(1 + ζ) w5 = 8 (1 − ξ)(1 + η)(1 + ζ) w7 = 8 w1 =
6.3.2
(1 + ξ)(1 − η)(1 − ζ) 8 (1 + ξ)(1 + η)(1 − ζ) w4 = 8 (1 + ξ)(1 − η)(1 + ζ) w6 = 8 (1 + ξ)(1 + η)(1 + ζ) w8 = 8 w2 =
(6.22)
Tetrahedral element η2
1(x1 , y1 , z1 )
(0, 1, 0) 2(x2 , y2 , z2 ) (x, y, z)
=⇒
η1
(1, 0, 0) 3(x3 , y3 , z3 )
4(x4 , y4 , z4 )
η3
(0, 0, 1)
Figure 6.16: Tetrahedral element and its transformation to local coordinates
A tetrahedral element is a volume with four faces and is analogous to a triangle in two dimensions. The derivation of weight functions for the volume element is similar to the one for triangles. Planes forming the volume are analogous to the lines forming the triangle. Similar to the case of a triangle, the tetrahedron can be transformed to local coordinates (η 1 , η 2 , η 3 , η 4 ). The transformation maps the actual tetrahedral unto right angled tetrahedral as shown in Figure 6.16 (similar to transformation of triangular coordinates) η 1 is the ratio of the distance of the point and distance of node 1 from plane opposite to the node i.e. plane containing nodes 2, 3 and 4. It is also the ratio of the volume of the tetrahedron made by the point with the plane containing nodes 2, 3 and 4 (shaded volume in Figure 6.17) to the volume of the tetrahedral element. It can be easily inferred that for a linear tetrahedral element w i = η i .
274
Chapter.6 Interpolation in two and three dimensions 1(1, 0, 0, 0) Figure
2(0, 1, 0, 0)
to
local
6.17:
Transformation
coordinates:
ratio
of
volume of the shaded tetrahedral to volume of tetrahedral element
3(0, 0, 1, 0)
4(0, 0, 0, 1)
Concluding remarks In this chapter, we have introduced, interpolation in two and three dimensions. These are very important since they are applied in all programs written for solving field problems in engineering. In fact the developments in interpolation techniques have paralleled the developments in FVM, FEM and other commercial software.