A study of solving two-dimensional potential problems using trial functions defined on finite elements and global trial functions defined on boundary elements

Advances in Engineering Software 15 (1992) 15-25

A study of solving two-dimensional potential problems using trial functions defined on finite elements and global trial functions defined on boundary elements Cindy Shrout California State University, Fullerton, California, USA and Williamson and Schmid, Tustin, California, USA

Sharon Gee Rockwell International, Seal Beach, California, USA

&

Dennis McGuigan Interstate Electronics Corporation, Anaheim, California, USA There are many methods that solve potential problems in two dimensions. Two such methods will be described. The first method uses trial functions defined on finite elements. Here, the given problem domain is discretized into finite elements. Then, evaluation points are determined on the boundary (in terms of the boundary conditions) and in the interior (in terms of the linear operator applied to the interior). The second method uses global trial functions defined on boundary elements. This method does not require the interior evaluation points, but rather depends on the evaluation points placed on the boundary. Test cases were run for both of these methods, and a discussion of the results and the methods are included in this paper.

L2(cl(f~), d#), which has the inner product (f,g)

INTRODUCTION

defined in the following manner: The Best Approximation Method can be used to solve modern-day problems involving linear operators, such as steady-state heat transfer and groundwater flow, among other problems. The method of solving two-dimensional potential problems (and other linear operator equations) by the Best Approximation Method requires an understanding of inner products and generalized Fourier series as applied to the following setting:

(f, g) = Ifg d# One way to define the necessary inner product for the development of a generalized Fourier series is to let # be one measure #1 on f~ and another measure #2 on F. One choice for a plane region would be for #~ to be the usual two-dimensional Lebesgue measure df~ on ~2 and for #2 to be the usual arc length measure dF on F. Then, an inner product is given by the following:

(f'g) = I~fg dFt + j fg

DEVELOPMENT Let f~ be a region in R m with boundary F and denote the closure of ~ by cl(f~). Consider the Hilbert space

Consider a boundary value problem consisting of an operator L defined on domain D(L) contained in L2(f~) and mapping it into L2(f~), and a boundary condition operator B defined on domain D(B) contained in L2(f~ ) and mapping it into L2(r). The domains of L and B

Advances in Engineering Software 0965-9978/92/$05.00 © 1992 Elsevier Science Publishers Ltd. 15

C. Shrout, S. Gee, D. McGuigan

16

have to be chosen so that for f i n D(L), Lfis in L2(f~), and for f i n D(B), Bfis in L2(F). With this restriction, we can construct a linear operator T mapping its domain D ( T ) = D(L)t_JD(B) into L2(cl(f~)) by the following:

Tf(x) =Lf(x) for x E a Tf(x) =Bf(x) for x on F. Then there is a single operator T on the Hilbert space L2(cl(f~)) which incorporates both the operator L and the boundary condition operator B, and which is linear if both L and B are linear. Given a linear operator relationship L(O) = h on f~ and auxiliary conditions O = O b on the boundary F. Then, by choosing a set of m linearly independent functions {~}m which spans an m-dimensional space S m, we have the inner product defined for elements of S m by (u,v) for u,v ~ S m in the following manner: (u,v)=Jr uvdF+I~

LuLvd~'l.

Clearly, this satisfies all the necessary conditions for inner products. From this, we can define the norm as follows:

I1 11

=

EJ~gid~i.

Generalized Fourier coefficients For a Hilbert space H, let V be a subspace that is spanned by n linearly independent vectors denoted {Yi}i~=l, Then, by using the Gram-Schmidt orthogonalization procedure, 1'3 a set of orthonormal vectors {~i}i~=1 can be found which also spans V. Let x be a given function that is to be approximated in the subspace V. Then, the best approximation of x in V, in the sense that the norm of the error (using the chosen inner product and associated norm) is minimized, can be written as a linear combination of the orthonormal vectors: x' = EcjO~,

where c~ = (x, ~3j).

Now, (x'Ok) = (EcjOj, Ok), which can be rewritten as the following:

(x',

The weighted inner product The choice of inner product above treats the boundary conditions and the interior approximation as being of equal importance (that is, the inner product is just the sum of the two). In most cases, it is desirable to weight the terms of the inner product. For example, due to dimensionality, the error for either the boundary or the interior is usually much larger in magnitude than the other. An obvious way to minimize the influence is to weight the inner product as the following:

(f,g)

.Io where0
+e.I/g dr, 1.

By choosing e > 0.5, the emphasis is given to matching the boundary conditions, while for e < 0-5, the interior approximation is emphasized (for e = 0-5, the interior and boundary approximations are weighted the same).

(u,u) 112.

The generalized Fourier series approach can now be used to obtain the 'best' approximation x' E S m using the inner product and this corresponding norm. Numerical quadrature is used to approximate this inner product, f~ is broken up into elements. The inner product is then given by:

Jfg

Since the set of vectors {Oj} is orthonormal, this simplifies to ck = (x',Ok). The coefficients {ok} thus determined are called the generalized Fourier coefficients 2.

:

A P P R O X I M A T I O N BY TRIAL F U N C T I O N S DEFINED ON FINITE ELEMENTS Consider the set of finite elements { Q i } , i = 1~2 , . . . , q, consisting of convex quadrilaterals whose union approximates the problem domain f~. In each finite element, an approximation of x (the function we are approximating over the entire domain) will be generated. Vector representation of the trial functions at a given node To simplify the bookkeeping, a geometry vector is initially stated. Its purpose is to specify the order in which boundary and interior evaluation points are represented in the associated function vectors. In this paper, the evaluation points on the boundary are listed first and then the interior points are listed. Let k represent the number of nodes for a specified problem. Then there are k function vectors {Yi}, each having dimension equal to the sum of the number of evaluation points on the boundary plus the number of interior points (in the entire domain). For each element, we will use eight nodes. One node will be placed at each vertex and one at each midpoint between adjacent vertices of the quadrilateral (see Fig. 1). Consider the kth element in our discretization of the domain. Also consider the set of eight basis functions { 1, x, y, x 2, y2, x2y, xy2, x2y2}. Let Nik(x, y) = ai + a2x +

A study o f solving two-dimensional potential problems

(I,3) (0,3) I

8

2

I

3

(2,3) 19

23

30

34

I:I

8

24

7

35

25

31

36

14

9

20

N~'(X'Y)t,

(3,3)

12

(0,2)

/

~ _I

(3,2)

(0,1)

I0

5

15

9

26

6

37

16

21

27

32

38

Y

..

25

/ / ~

4

17

II

a~//"

,,' /

ii

31

^

!

/

/ ~ 1

"/

t ~

(3,1)

Fig. 3. Trial function for node 27 (figure 1, rectangle 9). 6

17

4

8

22

7 (0,0)

,0)

28

5

39

29

33

4O

(2,0)

(3,0)

Fig. 1. a3y + a4 x2 + asy ~ + a6x2y + a7xy 2 + a8x2y 2 be the trial function for node i on element k. Now, set Ni~(x,y) to be 1 at node i, 0 at each of the other seven nodes in element k, and 0 outside of element k (see Figs 2 and 3). Let (xi, yi) be the coordinates of node i. Then the coefficients of N¢ can be determined by the matrix equation: xlylxlylxly

1

al

l

X2

Y2

X~

y~

2 2 Y22 X22 y2X2 y2X2

a2

l

x8

Y8

X~

y~

x ~ y 8 x 8 ysx8 2 2 Y82

8

Determining the generalized Fourier coefficients

1

Likewise, the coefficients of the basis functions can be determined for the other seven trial functions (N2~,..., Nsk) for this finite element. Once the coefficients {ai}, i = l, 2 , . . . , 8, of the trial functions for a node are found, then the coordinates of the evaluation points on the boundary of the element are plugged into the trial function. The result is then stored in the appropriate location in the function vector for that node (Yi), a s specified by the geometry vector.

Na(X,Y)

/

Then, the linear operator is applied to the trial functions associated with the node, and the coordinates of the interior evaluation points for the finite element are plugged in, with the result being stored into the function vector associated with the node being worked on, once again, as specified by the geometry vector. Likewise, this process is repeated for each node in the finite element being worked on. Once the process is completed for this finite element, each of the other finite elements are done in turn. Finally, a global geometry matrix is formed from the function vectors { Yi} just formed.

Y J

Fig. 2. Trial function for node 21 (figure 1, rectangle 9).

Let G = (si,.. •, Sm, q~,.. •, qn) be the designated geometry vector for a potential problem with m boundary evaluation points whose exact values are known, and n interior evaluation points. Let k denote the number of nodes obtained from discretizing the domain into p finite elements. Define a vector F of dimension rn + n, whose first m elements are the given values for the boundary evaluation points ordered according to the geometry vector G. The last n elements of the vector are the exact values of the linear equation L f = h applied to the interior evaluation points, again ordered according to G. In the case of potential problems, this exact value is known to be identically zero for all interior points, and so the last n elements of the exact solution vector F will always be zero in this paper. Let { Yi}, i -- 1 , 2 , . . . , k be the set of vectors obtained by applying the trial functions for each finite element according to the method detailed above. Then, for node i, Yi is the function vector for approximating the function x at that node. Using the Gram-Schmidt orthogonalization procedure, a set of k orthonormal vectors is constructed from the set of vectors derived above. Let {hi) denote this set of orthonormal vectors. Then the generalized Fourier coefficients are given by the following: "vi = (F, hi) , for i = 1 , 2 , . . . , k, and where F is the exact solution vector described above.

C. Shrout, S. Gee, D. McGuigan

18

This is the best approximation in the subspace R k, in that this approximation minimizes the L2 error for the given weighted inner product. Using a back-substitution procedure, the Fourier coefficients can be expressed in terms of the original basis: ~7-1h I :

~]'r~y 1

The trial functions used to construct these vectors were chosen to be identically zero at each node except the node corresponding to the function vector Yl. In addition, the trial functions were assigned the value zero outside of the subdomain covered by the element that it is defined in. The evaluation of E~-/*yiat node i is thus reduced to evaluating the ith term, as all other function values are zero at this node. Since the ith trial function was chosen to be one at node i, the approximation of the potential function at node i is just ~-/*. Error estimate for the finite element trial functions method

The solution of the least-squares Best Approximation Problem is unique and can be written as: !

*

x = E'~.y/, where ~ are the Fourier coefficients. As an error bound, we can use the Pythagorean Theorem to get the result known as Bessel's Inequality:

X2 >-- ~(gi, x) 2 In the problem at hand, we express Bessel's Inequality the following way:

(X, X) >_ ~(gi,

X) 2

But,

(x,x) = J'x2 dI" + J(Lx)2 d~, and

E(gi, x) = Ea; 2.

examination of the inner product used:

(f,g) = Ivfg dF + I~ LfLg dQ Since the linear operator is the Laplacian operator, choosing a basis satisfying Laplace's equations ensures that the second term is always zero. The problem is then reduced to one of minimizing the error on the boundary. Consider the set of boundary elements {Qi} i = l , . . . q , consisting of line segments whose union approximates the problem boundary. Define nodes at the endpoints, between adjacent line segments. In this method, global trial functions will be used in the process of approximating a function x over the problem domain. Let x' be the approximation obtained for x. This method will generate the value of x at any point interior to, on, or exterior to the problem boundary. After the nodes have been chosen, evaluation points must be selected along the problem boundary. Then, for each node, a branch cut, which is a ray that starts out at the node, and points away from the interior of ~, is defined. Also, a trial function is defined for each node. For this paper, we will consider the trial functions as functions of two dependent Variables: Ri/, the measure of distance from node i to evaluation point j and @i, the angle measured counterclockwise, whose initial side is the nodal branch cut and whose terminal side is a segment drawn from the evaluation point j to node i. Vector representation of the trial functions at a given node

The geometry vector for this method is simply an ordering of the evaluation points chosen along the boundary. The trial function for each node is represented by a vector whose length equals the number of evaluation points. The trial function is evaluated at each of the evaluation points. Then, the value obtained is inserted in the function vector in the position specified by the geometry vector. After the trial function vectors have been formed, they are stored in a global geometry matrix. Determining the generalized Fourier coefficients

So, I x 2 dE + I(Lx)2 dR >_ ~]a] 2.

A P P R O X I M A T I O N BY G L O B A L TRIAL F U N C T I O N S D E F I N E D ON B O U N D A R Y ELEMENTS This method uses global trial functions that are constructed from basis elements that satisfy Laplace's equations. The reason behind this is clear upon

The generalized Fourier coefficients are determined the same way as they were in the first metod, with two exceptions. First, n, the number of interior evaluation points, is equal to zero. The second difference is in the determination of the final approximation function. Let ~-~* be the Fourier coefficients obtained after the backsubstitution process. The approximation at any point exterior to, interior to, or on the boundary of ~ can be determined by summing up the products of the ~-*s and the basis functions for each of the nodes:

~j~ = ~T~B(Rij , Oij)

A study of solving two-dimensionalpotential problems w h e r e j varies from one to the number of result points, i varies from one to the number of nodes, and B(Rij, Oij) is the given basis function whose variables are based on a result point and a node.

Error estimate for real variable boundary element method Because we are still dealing with the least squares Best Approximation Problem, the error bounds can be found by using Bessel's Inequality like in the method using finite elements.

USER'S GUIDE The input data for the PASCAL programs describes the domain of the region to be applied on by the approximation by finite elements or by boundary elements and will define the division of the domain and the set-up of the associated geometry vector or boundary condition vector. There are three methods of inputting this data: l) create a data file (named TESTFILE.PAS) which the program calls; 2) input the data interactively while the program is running; 3) use the data from the last run, which is stored in a file called PREVIOUS.PAS.

following options: A) f(x,y) = 2xy B) f ( x , y ) = x ~ - y~ C) f(x,y) = 4x - 5y D) f(x,y) is unknown; only values at certain points are known. If this option is selected, then the user will input the F vector; 11) the nodes for each finite element and the side number, if any, that is on the boundary. When creating a data file, the file should consist of the data listed above. When the program is run, the outputs are written to external files as follows: 1) O U T F I L E . F I L - - contains the solution vector, the exact solution (if known), the Laplacian approximation surface, and the exact surface. 2) F O U T . F I L - - contains the contents of the F vector and the calculated values in the F matrix. 3) A M A T O U T . F I L - - contains the calculated values in the A matrix. 4) T R I A L O U T . F I L - - contains the values calculated to be the coefficients of the basis functions used in the trial functions for each node. The computer code for this program may be obtained from the authors. For the real variable boundary element program, the following input is needed:

The program will prompt the user to select the method desired for a program run. For the finite elements program, the following input is needed: 1) the number of elements that the region is divided into (NumberOfElements); 2) the number of element sides that are on the boundary (DimBoundary); 3) the number of nodes that are element vertices (NumVertices); 4) the number of nodes that are midpoints of an element side (NumMidPoints); 5) the number of evaluation points that should be used for each side of an element that is on the boundary between a vertex and the midpoint (HalfNumEval); 6) the node number (j) and the coordinates of each of the vertex nodes (NC[ j][1], NC[ j][2]); 7) the node number (j) and the node numbers of the adjacent vertex nodes of each midpoint node; 8) the weighting factor for the norm and inner product (e2); 9) whether the exact solution is known (m); 10) how elements of the F vector will be evaluated. For this input, the user will select from the

19

2) 3) 4) 5) 6)

7)

8)

the choice of basis functions (BasisFunct). For this input, the user will select from the following options: A) B(Rj, Oj) = R~(cos(Oj) ln(Rj) - O~ sin(Oj)) B) B(Rj, Oi) = R2(sin(Oj)ln(R~) - O~cos(Oj)) C) ln(R~); the number of nodes (NumberOfNodes); the number of evaluation points between any two adjacent nodes (NumEval); the node number (j), the coordinates of the nodes (NC[ Jill], NC[ j][2]), and the angle of the branch cut for each node (BranchCuts[j]); the number of points in the list of results (NumResult); the list of results numbering system and the coordinates of the point to be included in the list of results (j, ResultPoints[ j][1], ResultPoints[ j][2]); the boundary conditions (BoundaryFunct). For this input, the user will select from the options listed above for selecting the F vector input for the finite elements method; whether the exact solution is known (t). When creating a data file, the file should consist of the data listed above.

C. Shrout, S. Gee, D. McGuigan

20

W h e n the p r o g r a m is run, the outputs are written to external files as follows: 1) F _ O U T F I L E . F I L - - will contain the F matrix, the F vector, and the G a m m a vector. 2) O U T F I L E . F I L - - will contain the list o f results and the exact solution (if known). The c o m p u t e r code for this p r o g r a m is found in the Appendix

TEST CASES The following section documents the work done on two separate regions used as input spaces. Case 1

The d o m a i n o f this case is a square with sides along the x and y axes. F o r the finite elements program, it was divided up into nine elements, all squares o f uniform size. The F vector was set up so that the nodes on the left side o f the square contained the value o f 100, the nodes on the right side o f the square contained the value o f zero, and the remaining nodes contained values that are linear in nature. Figure 1 shows the discretized d o m a i n of this region for the finite element method.

F o r the real variable b o u n d a r y element program, the b o u n d a r y conditions were m a d e to be identical to that o f the finite elements p r o g r a m run. The test results for this run are shown in Table 1. The region in Case 1 is a slightly more complex one than the one shown in Fig. 5. We discretized it into more finite elements, and we also had one element entirely on the interior o f the region. The calculated solution vector for the finite elements p r o g r a m reveals that we are able to approximate the function values at the nodes to a high degree o f accuracy (see Table 1). F o r the real-variable b o u n d a r y element program, Case 1 is virtually identical to the type o f region in Fig. 5. W h e n the basis function is o f the form Rj(sin(Oj)ln(Rj)-Ojcos(~)j)) or o f the form Rj(cos(O2) ln(Rj)- O~sin(Oj)), then we get a fair approximation to the function values at the nodal points on the b o u n d a r y o f the surface. W h e n the basis function is o f the form ln(R~), the approximation isn't quite as good. One reason for this is that the basis function doesn't take into consideration the angle between the branch cut and the evaluation points. The results o f Case 1 illustrate a very important point: the selection o f the basis functions is very important in solving problems o f this type. F o r the finite elements program, the basis functions used were

Table 1. Output for region in Fig. 1 (Case 1)

Point coords (0,3) (0,2'5) (0,2) (0,1'5) (0,1) (0,0"5) (0,0) (0.5,0) (1,0) (1.5,0) (2,0) (2.5,0) (3,0) (3,0.5) (3,1) (3,1 '5) (3,2) (3,2.5) (3,3) (2.5,3) (2,3) (1-5,3) (1,3) (0-5,3)

F.E. approx,

B.E. approx. (at

B.E. approx. (b)

B.E. approx (c/

100

100.000 00

100 100 100 100 100 100

100'000 00 100'000 00 100.000 00 100.000 00 100'000 00 100.000 00

100.332 46 100.000 26 99.982 17 100.019 59 100-004 12 99"998 31 100.145 08 83.33422 66.658 84 50.023 37 33.347 83 16.662 15 0.271 80 0.000 09 -0"013 02 0.074 36 0.039 22 -0.008 57 0-459 18 16-666 13 33-310 31 50.070 58 66.695 50 82.32760

100-172 10 100-008 71 99"996 36 100-000 36 100.002 01 99'998 16 99.938 80 83.32420 66.675 08 50.001 08 33-319 32 16'66763 -0.529 55 -0.026 45 0.012 07 -0'000 02 -0.013 02 0.002 85 -0.296 25 16.65806 33.333 34 49-999 26 66-669 67 83-33338

99.963 46 100.005 78 100.007 47 98.211 58 100.007 47 100.005 78 99.963 46 83.311 11 66-680 14 53.920 66 33.282 51 16.75203 27.962 91 0"057 36 -0.044 82 9-629 75 -0-044 82 0.057 36 27.962 91 16.666 13 33.282 51 53-920 66 66.680 14 83.311 11

Exact

83.33333 66'666 67 50 33.333 33 16"66667 0 0 0 0 0 0 0 16-66667 33-333 33 50 66.666 67 83.33333

83.33333 66"666 66 50'000 00 33.333 34 16.66667 0-000 00 0-000 00 0-000 00 0"000 00 0.000 00

0'000 00 0.000 00

16.66667 33.333 34 50.000 00 66.666 66 83.33333

(a)basis function = Rj(cos(Oj)ln(Rj)- Oj sin(Oj)) (b)basis function Rj(sin(Oj)ln(Rj) - Oj cos(Oj)) (")basis function ln(Rj)

A study of solving two-dimensional potential problems

(2,3) 6

(0,25) (0,2)

(5,3)

9

14 5

18

II

15

19

l ~ I I ° 4

2

1751

(2,1 5)

5 ~ 1 2 (0,

(775,1)

20

3 13

21 (5,0)

16

(2,0)

Fig. 4.

fl(x,y) = 1, f 2 ( x , y ) = x ,

f 3 ( x , y ) = y , f a ( x , y ) = X 2, f6(x,y) = x2y, f7(x,y) xy 2, and fs(x,y) = x2y 2. Using these basis functions, we can approximate any function that is a linear combination of these functions exactly as our results show. Any function being approximated that is not a linear combination of these functions can still be approximated; however, the calculated values won't be equal to the exact values of the function (to be approximated) at the given nodes. fs(x,y)

y2,

=

Case 2 The domain of this case is an unusually shaped eightsided figure with two sides that lay partially on the x and y axes. F o r the finite elements program, it was divided up into five quadrilaterals of various sizes and shapes. The F vector was set up such that the nodes on the boundary were defined using the function f ( x , y ) - - x 2 _ y 2 . Figure 4 shows the discretized domain of the region for the finite element method. For the real variable boundary element program, the boundary conditions were made to be identical to that of the finite elements program run. The test results for this run are in Table 2 below.

Table 2. Output for region in Fig. 4 (Case 2) Point coords

Exact

F.E. approx.

B.E. approx. (a)

(2,3) (1,2'75) (0,2) (0,1'25) (0,0"5) (1,0"25) (2,0) (3-5,0) (5,0) (6"375,0'5) (7'75,1) (8"375,1'375) (9,1'75) (7,2'375) (5,3) (3"5,3) (2,3) (1,2"75)

-6'25 -5"062 5 -4 -1"5625 -0"25 0"937 5 4 12"25 25 40-640625 59"0625 68-25 77-937 5 43"359 375 16 3"25 -5 -6'562 5

-6"250 00 -5'062 50 -4"000 00 -1"56250 -0"250 00 0"937 50 4"000 00 12'250 00 25-000 00 40-39063 59"06250 68-25000 77"937 50 43"35938 16'00000 3"250000 -5'00000 -6'562 50

-6"209 35 -5'046 01 -3"976 90 -1'51466 -0" 122 84 0'905 64 3'959 27 12'201 26 24"942 60 40"359 62 59"06547 68"24285 77"767 03 43"249 66 15'88606 3"171 61 -5'06461 -6-603 49

(a)basis function = Rj(COS(Oj)ln(Rj)- Oysin(Oj))

21

In Case 2, we modified the shape of the region so that we would not have all of the sides parallel to the X-axis or the Y-axis. Examination of the values in Table 2 reveals that, with the finite elements program, we can approximate the function f ( x ) = x 2 _y2 exactly at all of the nodal points. This is because the function being approximated can be written as a linear combination of the basis functions. For the real-variable boundary element program, we get good approximations (to the function f(x) = x 2 - y2 along the boundary of the surface) for basis functions of the form Rj(sin(Oj)ln(RflOjCOS(Oj))and Rj(cos(O~)ln(R~) - Oj sin(Off).

ERROR ANALYSIS For this research, we used the Laplacian operator to approximate the values of the interior points. We set L ~ = h such that we have the homogeneous case of h = 0 (i.e. L ~ = 0). This simplifies Bessel's Inequality to: Jr ~2 d£ >_ ~ a ; 2 For

the

domain

in

Fig.

5

using

the

function

f(x,y) = 2xy for the boundary condition, the line integral on the left side of the inequality above evaluates to 42, while the summation on the right-hand side of the inequality was calculated to be 34.824 73 for the method involving finite elements. For the method involving boundary elements, the right-hand side was calculated to be 33.827 82. The error estimate is given by ~ : ( ~ , ~ ) - - ~ ( a ; ) 2, which for the two methods is 7.17527 and 8.172 18, respectively. There are three ways to improve the error for these methods. The first is by increasing the number of evaluation points which are used in calculating the approximation. The next way to improve the error is to increase the number of nodes, either on the boundary, in the domain, or both of these. The third way to improve the error is to use a different set of basis functions. Our efforts reveal that the selection of basis functions is very

(0,1)

(05,1)

(0,05)

I

3 (0,0)

(I,I)

(I5,1)

(1,0.5) 2

5 (05,0)

10 (1,0)

(15,0)

Fig. 5.

(2,1) 12 (2,051

13 (2,0)

22 (0,1)

C. Shrout, S. Gee, D. McGuigan ~

.

~

r

- -

--

(2,1)

The following listing is from MAIN.PAS, and contains the main routine and the functions and procedures in the program (for brevity, we did not include the procedure gramonly, which is responsible for performing the Gram-Schmidt orthogonalization; it is part of the source code used for the method involving finite elements): *************************************************************************

(o,o) "

. . . .

(2,o)

* * *

Sharon Gee Dennis McGuigan Cindy Shrout

* Program Units * * *

Fig. 6. Approximate boundary for Fig. 5. important to calculating a reasonable approximation.

-

* * * * * * * * * *

Main Program Procedure CheckForZero Function Distance Procedure InpufData Procedure ReadData

Procedure ReadLast

* *

Function Calculate Function FCalculate * * Procedure PrintOut * **************************************************************************

Approximate boundary

Program B o u n d a r y E l e m e n t s ;

The real-variable boundary element program estimates the values of ff at the result points asked for. As a measure of error, this procedure may be reversed to find the points that match the boundary conditions exactly. If enough of the points are found, then an approximate boundary can be drawn. As an example of this, the domain in Fig. 5 was used, as was the function f ( x , y ) = 2xy to determine the boundary values. The resulting approximate boundary is shown in Fig. 6. It shows that the approximate boundary matches the domain boundary very well, except where x gets close to, but is not equal to, zero.

CLOSING REMARKS

************************************************************************* THIS PROGRAM APPROXIMATES 2 DIMENSIONAL REAL VARIABLE BOUNDARY ELEMENTS ************************************************************************

const M a x R o d e s = 50; M a x E v a i P t s = 100: type C o o r d T y p e = a r r a y [1..2] of reai~ V e c t o r T y p e = a r r a y [ l . . M a x E v a I P t s ] o£ real; M a t r l x T y p e = a r r a y [ l . . M a x N o d e s 3 of V e c t o r T y p e ~ E v a i T y p e = a r r a y [ l . . M a x E v a l P t s ] of C o o r d T y p e ; S o l u T y p e = a r r a y [ l . . M a x N o d e s ] of real; vat F,B,W,V text; NC array[l.,MaxNodes] BranchCut s array[l..MaxNodes] OammaVecLor soluType; I, J, DimF, B o u n d a r y F u n c t integer; N u m E r a l , N u m b e r o fN o d e s Integer; ptr, N u m R e s u l t , B a s i s F u n c t integer; PX, PY, QX, QY, R J real; T, U char; Zero, exact, again boolean; FMat rix matrixtype; FVector, Boundar yValue vectortype; E v a l P o l n t s , Resul t s P o l n t s evaltype; {* B r l n g in p r o c e d u r e s and functions {$I a : G r a m o n l y . p a s } {$I a : F i n d E v a l P o i n t s . p a s }

For the purposes of this paper, we used the Laplacian operator as our linear operator and set the value of the operator at all interior points equal to zero; however, it should be noted that the method involving finite elements can be used to approximate functions using many other linear operators just as easily.

of C o o r d T y p e ; of i n t e g e r ;

from o t h e r f l i e s x)

(* I n c r e a s e s t a c k s i z e for s u b r o u t i n e s {$M 6 5 5 2 0 , 0 , 6 5 5 5 6 0 }

*)

******************************************************************** (* P R O C E D U R E C H E C K FOR Z E R O *) {* ................. *) (* T h i s p r o c e d u r e d e t e r m i n e s w h e t h e r a n u m b e r is so c l o s e to *) (* zero that if will he c o n s i d e r e d as such. *) ******************************************************************** Procedure CheckForZero(number:

real};

begin Z e r o : = false; If { n u m b e r < 0 . 0 0 0 0 0 1 ) Z e r o : : true~

and

{number

> -0.000001)

then

end~ ******************************************************************** (* (*

REFERENCES 1. Hromodka II, T.V., Yen, C.C. & Pinder, G.F. Approximation M e t h o d - An Introduction, Berlin, 1987. 2. Luenberger, D.G. Optimization By Vector Space John Wiley & Sons, New York, USA, 1969. 3. Johnson, L.W. & Riess, R.D. Numerical Addison-Wesley Publishing Company, Inc., Mass., 1982.

FUNCTION DISTANCE .................

*) *)

(* T h l s f u n c t i o n c a l c u l a t e s the d i s t a n c e b e t w e e n two p o l n t s *) (* in a plane. *) ********************************************************************

The Best Springer, Methods, Analysis, Reading,

APPENDIX This appendix contains the source code used to implement the method using global trial functions defined on boundary elements in approximating twodimensional potential problems.

Function Distance(Xl,Yi,X2,Y2

: real)

: real;

begin distance end;

:= s q r t ( s q r ( X 2

- Xl)

+ sqr(Y2

- YI))~

******************************************************************** (* PROCEDUHE INPUTDATA *) (* ................... *) (* T h i s p r o c e d u r e is u s e d to input d a t a w h e n the o p e r a t o r has *) (* s e l e c t e d i n t e r a c t i v e input. *} ******************************************************************** Procedure InputData; vat J, I : i n t e g e r ;

begin w r i t e l n ( ' C h o o s e the basls functlon d e s i r e d . ' } : wrlfeln; writeln(' i. R * c o s ( t h e t a ) * In(E) sin(theta)'); wrlteln~ writeln{' 2. E * s i n ( t h e t a } * In(R) cos{theta)'): wrlteln: w r l t e l n { ' 3. In(R)'): writeln; readln(BaslsPunct)~ wrlteln(B,BaslsFunct); writeln; w r l t e { ' E n t e r the n u m b e r of nodes: '); readln(NumberofNodes); writeln(B,NumbetofNodes); writeln~

-

R

*

theta

*

+

R

*

theta

*

A study of solving two-dimensionalpotential problems w r i t e l n ( ' E n t e r the n u m b e ~ of e v a l u a t i o n ' ) ; w r i t e ( ' p o i n t s b e t w e e n any 2 a d j a c e n t nodes: '); readln(NumEval); writeln(B,NumEval); for ] :: 1 to N u m b e r o f N o d e s do begin w r i t e ( ' E n t e r the n u m b e r c o r r e s p o n d i n g to the node: '); readln(J); writeln(B,J): w r i t e ( ' E n t e r the x c o o r d i n a t e of n o d e ' , J , ' : '); readln(NC[J][1]); writeln(B,NC[J][l]); w r i t e ( ' E n t e r the y c o o r d i n a t e of n o d e ' , J , ' : '); readln(NC[J][2]); writeln(B,NC[J][2]); end: wrlteln; w r i t e l n ( ' T h e f o l l o w i n g s e c t i o n is to e n t e r the angle') w r l t e l n ( ' o f the b r a n c h cut for e a c h n o d e . ' ) ; for I := 1 to N u m b e r o f R o d e s do begin w r i t e ( ' E n t e r the N o d e Number: ')~ readln(J); writeln(B,J)~ w r l t e l n ( ' E n t e r the b r a n c h cut a n g l e for n o d e ' , I ) readln(BranchCuts[I]); wrlteln(B,BranchCuts[I]); end: wrlteln; w r i t e l n ( ' T h e f o l l o w i n g s e c t i o n is to enter the p o i n t s ' ) ; w r i t e ! n ( ' f o r the list o£ r e s u l t s . ' ) ; w r l t e l n ( ' E n t e r the n u m b e r of p o i n t s in the list of r e s u l t s . ' ) ~ readln(NumResult); wrlt e l n ( B , N u m R e s u ~ t ) ; For [ :: i to N u m R e s u l t do begin w r l t e l n ( ' E n t e r the result p o l n t n u m b e r . ' ) ; readln(J); writeln(B,J); w r l t e l n ( ' E n t e r the x c o o r d i n a t e for this p o i n t . ' ) ; readln(ResultsPolnts[J][l]); wrlteln(B,ResultsPoints[J~[i]); w r l t e l n ( ' E n t e r the y c o o r d l n a t e fer thls p o l n t . ' ) ~ read[n(EesultsPolnts[J][2]); writeln(B,Resu]tsPoints[J~[2]); end; writeln('Choose the b o u n d a r y f u n c t i o n d e s i r e d . ' ) ; writeln; w r i t e l n ( ' i. 2 * X * Y'); writeln; w r i t e ! n ( ' 2. X * X - Y * Y'); wrlteln; w r i t e l n ( ' 3. 4 * X - 5 * Y'); writeln; w r i t e l n ( ' 4. U s e r input point by p o i n t ' ) ; writeln; read]n(BoundaryFunct); writeln(B,BoundaryFunct); If ( B o u n d a r y F u n c t = 4) t h e n begln F o r J :: 1 to (NumEral * N u m b e r O f N o d e s ) do begin w r i t e ( ' E n t e r the e v a l u a t i o n point n u m b e r : '): readln(J); writeln(B,J); w r i t e ( ' E n t e r the b o u n d a r y v a l u e at p o i n t ',J); readln(BoundaryValue[J]): writeln(B,BoundaryValue[J3)~ end~ end: w r i t e l n ( ' I s the exact s o l u t i o n k n o w n ? (Y/N)'); readln(T); writeln(S,T); If (T = 'Y') or (T = 'y') t h e n exact := true; end: ******************************************************************** (* P R O C E D U R E R E A D D A T A *) (* .................. *) (* T h i s p r o c e d u r e is u s e d to i n p u t d a t ~ w h e n the o p e r a t o r has *) (* d e c l d e d to u s e t h e d a t a s t o r e d in t e s t f i ! e . p a s . *) ******************************************************************** Procedure ReadData; vat I,J : i n t e g e r ; begln readln(F,BasisFunct); readln(F,NumberofNodes); readln(F,NumEval); for I :: i to N u m b e r o f N o d e s do begin readln(F,J); readln(F,NC[J][l]); readln(F,NC[J][2]); end; for I := i to N u m b e r o f N o d e s do begin readln(F,J); readln(F,BranchCuts[J]); end; readl n ( F , N u m R e s u l t ) ; for I := i to N u m R e s u l t do begin readln(F,J); readln(F,Resu~tsPoints[J][1])~ readln(F,ResultsPoints[J~[2~); end~ readln(F,BoundaryPunct); If ( B o u n d a r y F u n c t = 4) t h e n begin For J := 1 to (NumEval * N u m b e r O f N o d e s ) do begin readln{F,J); readln(F,BoundaryValue[J]); end; end; readln(F,T): (* Y if exact s o l u t i o n is k n o w n *) If (T : 'Y') or (T = 'y') t h e n e x a c t := true; end; ******************************************************************** (* PROCEDURE E E A D L A S T (* .................. (* T h i s p r o c e d u r e is u s e d to i n p u t d a t a w h e n the o p e r a t o r has (* d e c i d e d to use the d a t a s t o r e d in p r e v i o u s . p a s . ******************************************************************** Procedure ReadLast; vat

*) *) *) *)

I,J : integer; begln readln(B,BasisFunct)~ readln(B,NumberofNodes): read|n(B,NumEval); for [ :: i to N u m b e r o f N o d e s do begin readln(B,J); readln(B,NC[J][l]); readln(B,NC[J][2]); end; for i :- 1 to N u m b e r o f N o d e s do begln readln(B,J): readln(B,BranchCuts[J])~ end: readln(B,NumResult)~ for i := i to N u m R e s u l t do begln readln(B~J); readln(B,ResultsPolnts[J][l]) readln(B,ResultsPoints[J][23) end~ readln(B,BoundaryFunct); readln(B,T); (* Y if exact s o l u t i o n If (T = 'Y') or (T = 'y') t h e n exact :: true; end~

23

is k n o w n

*)

******************************************************************** (* F U N C T I O N C A L C U L A T E *) (* .................. ~> (* A n o d e and an e v a l u a t l o n p o i n t are s e n t to thls function. *) (* It first d e t e r m i n e s the a n g l e b e t w e e n the b r a n c h cut at the *) (* n o d e and the line s e g m e n t b e t w e e n the n o d e and the *) (* e v a l u a t i o n point. T h e n it c a l c u l a t e s the d ~ s t a n c e b e t w e e n *) (* the n o d e a n d the e v a l u a t i o n point. F i n a l l y it e v a l u a t e s *) (* the basls f u n c t i o n for the a n g l e a n d d i s t a n c e found. *) ******************************************************************** F u n c t l o n C a l c u l a t e ( X l , YI, X2, Y2: ~eal; K : integer) : real; var TX, TY, temp, angle~ angleA, a n g l e B : real; a n g l e f i g u r e d , zero_x, z e r o _ y : boolean; begin a n g l e f i g u r e d := false; angle := 0.0; TX :: X2 - Xl; TY Y2 YI; CheckForZero(TX); Z e r o _ x :: zero; CheckForZero(TY); Z e r o _ y := zero; If z e r o _ x a n d z e r o _ y t h e n begin a n g l e := 0.0; a n g l e f i g u r e d := true; end; If not a n g l e f i g u r e d t h e n begin If Z e r o _ x then begin If Y1 > Y2 t h e n begin If B r a n c h C u t s [ K ] : 0 t h e n a n g l e := (3"Pi)/2; If B r a n c h C u t s [ K ] : 90 t h e n a n g l e :: pi; If B r a n c h C u t s [ K ] = 180 t h e n a n g l e := p1/2; If B r a n c h C u t s [ K ] : 270 t h e n a n g l e :: 0.0; a n g l e f i g u r e d := true; end~ If Y1 < Y2 t h e n begin If B r a n c h C u t s [ K ] = 0 t h e n a n g l e := pi/2; If B r a n c h C u t s [ K ] = 90 t h e n angle := 0.0; If B r a n c h c u t s [ K ] = 180 t h e n a n g l e := (3"PI)/2; If B r a n c h C u t s [ K ] = 270 then a n g l e := pi; a n g l e f l g u r e d := true: end; end~ end; If not a n g l e f i g u r e d then begin If Z e r o _ y t h e n begin If X1 > X2 t h e n begin If B r a n c h C u t s [ K 3 = 0 t h e n a n g l e := Pi; If B r a n c h C u t s [ K ] : 90 t h e n a n g l e := pi/2; If B r a n c h C u t s [ K ] = 180 t h e n a n g l e := 0.0~ If B r a n c h C u t s [ K 3 = 270 t h e n a n g l e :: (3"Pi)/2 a n g l e f i g u r e d := true; end; If Xl < X2 t h e n begin If B r a n c h C u t s [ K ] = 0 t h e n a n g l e := 0.0; If B r a n c h C u t s [ K ] = 90 t h e n a n g l e := (3"Pi)/2; If B r a n c h C u t s [ K ] = 180 t h e n a n g l e := Pi~ If B r a n c h c u t s [ K ] = 270 t h e n a n g l e := Pi/2~ a n g l e f i g u r e d :: true; end; end; end; If not a n g l e f i g u r e d t h e n begin t e m p := TY/TX; a n g l e A :: A r c T a n ( t e m p ) ; If (TY > 0.0) a n d (TX > 0.0) t h e n begin a n g ] e B := angleA; If B r a n c h C u t s [ K ] = 0 t h e n a n g l e := angleB; If B r a n c h C u t s [ K ] = 90 t h e n a n g l e := (3"Pi)/2 + angleB; If B r a n c h C u t s [ K ] = 180 t h e n a n g l e :: Pi + angleR; If B r a n c h C u t s [ K ~ : 270 t h e n a n g l e :: (Pi/2) + angleR; end; If (TY > 0.0) and (TX < 0.0) t h e n begin a n g l e r := Pi + angleA; [f B r a n c h C u t s [ K ] = 0 then a n g l e :: angleB~ If B r a n c h C u t s [ K ] = 90 t h e n a n g l e := a n g l e r - (Pi/2): If B r a n c h C u t s [ K ] = 180 then a n g l e := Pi + angleB; [f B r a n c h C u t s [ K ] : 270 t h e n a n g l e := (Pi/2) + angleR; end; If (TY < 0.0) and (TX < 0.0) t h e n begin a n g l e B := Pi + ang]eA; If B r a n c h C u t s [ K ] = 0 then a n g l e := angleR; If B r a n c h C u t s [ K ] = 90 t h e n a n g l e := a n g l e B (Pi/2); If B r a n c h C u t s [ K ] = 180 t h e n a n g l e := a n g l e B - Pl;

C. Shrout, S. Gee, D. McGuigan

24

If B r a n c h C u t s [ K ] = 270 then a n g l e :- (Pi/2) * angi~B; end; If (TY < 0.0) and (TX > 0.0) t h e n begin a n g ] e B := 2*P~ * angleA; If S r a n c h C u t s [ K ] = 0 then a n g l e := angleB; If B r a n c h C u t s [ K ] = 90 then a n g l e := a n g l e B - (Pi/2); If B r a n c h C u t s [ K ] = 180 then a n g l e := a n g l e B - Pl; If B r a n c h C u t s [ K ] = 270 t h e n a n g l e := a n g l e B - (3"PI)/2; end; end; RJ := D l s t a n c e ( X l , Y i , X 2 , Y 2 ) ; CheckForZero(R3); If z e r o t h e n C a l c u l a t e := 0.0 else begin C a s e B a s i s F u n c t of i: C a l c u l a t e :: RJ * ( c o s ( a n g l e ) * ~n(RJ) angle * sin(angle)); 2: C a l c u l a t e := EJ * ( s l n ( a n g l e ) * in(RJ) + angle * cos(angle)); ~: C a l c u l a t e :- in(EJ); end~ (* c a s e *) end;

NC[I][2], ResultsPo~nts[J)[l~, R e s u l t s P o i n t s [ J S [ 2 ] , I)): end; wrlteln(W, end~

end;

FUNCTION FCALCULATR ...................

*) *)

(* T h i s f u n c t i o n e v a l u a t e s the b o u n d a r y c o n d i t i o n at an *) (* e v a l u a t i o n point. *) ******************************************************************** F u n c t i o n F C a l c u l a t e ( X, Y : real; J : i n t e g e r ) : begin C a s e B o u n d a r y F u n c t of l: P C a l c u [ a t e := 2 * X * Y; 2: F C a l c u l a t e := X * X - Y * Y; 3: F C a l c u l a t e := 4 * X - 5 * Y; 4: F C a l c u l a t e := B o u n d a r y V a l u e [ J ] ; end~ (* case *) end;

•* * * * * * * * * * * * * * * * * * * * * *

real~

PROCEDURE PRINTOUT ..................

*) *)

(* T h i s p r o c e d u r e w r l t e s the F M a t r i x , the F V e c t o r a n d the *) (* s o l u t i o n v e c t o r g e n e r a t e d by the m o d i f i e d G r a m - S c h m i d t *) (" p r o c e d u r e to F _ O u t f i l e . p a s . T h e r e s u l t s , and the exact *) (* s o l u t i o n if k n o w n are w r i t t e n to o u t f i l e . p a s . *) ******************************************************************** Procedure PrintOut~ vat i,j : i n t e g e r ; temp : real; begin w r i t e ( W , 'The B a s i s F u n c t l o n s are : '); C a s e B a s i s F u n c t of i : w r i t e l n ( W , ' R j * ( c o s ( t h e t a ) * in(Rj) - t h e t a * s i n ( t h e t a ) ) ' ) ;

Main Program Body

*****************************

readln(U): writeln; If (U : 'Y') or (U = 'y') t h e n begin Reset(B); Readlast; end else begin Reset(F); ReadData; end~ end: D i m F := N u m E r a l * N u m b e r o f N o d e s ;

2 : w r i t e l n ( w , ' R j * ( s i n ( t h e t a ) * in(Ej) + t h e t a * c o s ( t h e t a ) ) ' ) : 3 : wrlteln(w,'in(Rj)'): end; (* c a s e *~ writeln(W); w r i t e ( V , 'The B a s i s F u n c t i o n s are : '}; C a s e B a s i s F u n c t o£ 1 : w r i t e ] n ( V , ' R j * ( c o s ( t h e t a ) * in(Rs) ~ thet~ * s i n ( t h e t a ) ) ' ) ; 2 : w r i t e l n ( V , ' R j * ( s i n ( t h e t a ) * In(Ej) + t h e t a * c o s ( t h e t a ) ) ' ) ;

I

:=

i~

ptr := !; w h i l e (I <= R u m b e r O f B o d e s ) do begln If (i = N u m b e r O f N o d e s ) t h e n F i n d E v a l _ P o i n t s ( p t r , i, I, N u m E r a l , E v a l P o i n t s ) else F i n d E v a l _ P o i n t s ( p t r , i, i+l, NttmEval, E v a l P o i n t s ) ;

3 : wrlteln(V,'in(Rj)'); end~ (* c a s e *) writeln(V); w r i t e ( W , 'The B o u n d a r y F u n c t i o n is : '); C a s e B o u n d a r y P u n c t of 1 : w r i t e l n ( w , '2 * X * Y'); 2 : w r i t e l n ( w , 'X * X - Y * Y')~ 3 : w r i t e l n ( w , '4 * X ~ 5 * Y'); 4 : w r i t e l n ( W , '~ser d e f i n e d p o i n t by p o i n t . ' ) ; end; (* o a s e *) writeln(W); w r i t e ( v , 'The B o u n d a r y F u n c t i o n is : '); C a s e B o u n d a r y P u n c t of 1 : w r l t e l n ( v , '2 * X * Y'); 2 : w r i t e l n ( V , 'X * X - Y * Y'); 5 : w r i t e l n ( v , '4 * X - 5 * Y'); 4 : w r l t e l n ( v , 'User d e f i n e d p o i n t by p o i n t . ' ) ; end; (* c a s e *} writeln(V); w r l t e ] n ( V , 'The F M a t r i x is -• ' ) ; For i := 1 to R u m b e r O f N o d e s do begln writeln(v); w ~ i t e l n ( v , ' F',i,' : '): for j:= i to D i m F do be~in if ((j m o d 8) = 0 ) t h e n writeln(V,' ',FMatrix[l~[J]:10:5) elS~rite(V,' ',FMatrix[I][J]:10:5): end~ end; writeln(V)~ writeln(V); w r i t e l n ( V , ' T h e F V e c t o r is : '); for j := 1 to D I m F do begin If ((j m o d 8) = 0) then w r i t e l n ( V , ' ', F V e c t o r [ $ ] : 1 0 : 5 ) else w r i t e ( V , ' ', P V e c t o r [ J ] : 1 0 : 5 ) end; writeln(V)~ w r i t e l n ( V , ' T h e s o l u t l o n v e c t o r is : '); write[n(V)7 for i := I to N u m b e r O f N o d e s do w r i t e l n ( V , ' x ' , i , ' = ',Ga,m~aVector[i]:10:5); writeln(V); If not e x a c t t h e n begln w r i t e l n ( W , 'The r e s u l t s llst is : '); writeln(W); for J := 1 to N u m R e s u l t do beglD t e m p := 0.0; for I :: i to ~ u m b e r O f R o d e s do begln t e m p := t e m p * (GanlmaVector[I] * C a l c u l a t e ( N C [ I ~ [ l ~ ,

- ',temp:10:5);

begln a s s i g n (F, ' a : t e s t f i l e . p a s ' ) ; a s s i g n (B, ' a : p r e v i o u s . p a s ' ) ; a s s i g n (W, ' a : o u t f l l e . f i i ' ) ; a s s i g n (V, ' a : F _ O u t f i l e . f i ! ' ) ; A g a i n := true; w h i l e A g a i n do begin rewrite(W)~ rewrite(V)~ w r i t e ( ' D o y o u want to i n p u t d a t a I n t e r a c t i v e l y ? (Y/N) '): readln(T); e x a c t :: false; write[n: if (T = 'Y') or (T = 'y') then begln rewrite(B): InputData; close(B); end else begin w r i t e ( ' D o y o u w a n t to r e - r u n last i n t e r a c t i v e test? (Y/N)

******************************************************************** (* (*

',J,'

end else begin w r i t e ( w , 'The r e s u l t s list is : '): writeln(W, ' T h e e x a c t s o l u t i o n is : '); writeln(W); for J := 1 to N u m R e s u l t do begin t e m p := 0.0; for I := 1 to N u m b e r O f N o d e s do begin temp := temp + (uam~naVector[I] * C a l c u l a t e ( N O [ f i l l ] , NC[I][2], ResultsPoints[J][l], R e s u l t s P o i n t s [ J ] [ 2 ] , I)); end; w r i t e ( w , 'Result ',J,' = ' , t e m p : 1 0 : 5 ) ; write(W, ' '); writeln(W,FCalcu]ate(NesultsPoints[J][l], ResultsFoints[J][2], J):10:5); end; end~ w r i t e l n ( ' P r l n t out c o m p l e t e d . ' ) ; end;

******************************************************************** (* (*

'Result

I

::

I

•

i;

~nd;

for

I := 1 to N u m b e r o f N o d e s do begin for J :: i to D i m F do begin PX :: N C [ I ] [ I ] ; PY :: N C [ I ] [ 2 ] ; QX := E v a ] P o i n t s [ J ] [ 1 3 : QY := E v a l P o i n t s [ J ] [ 2 ] ; F M a t r i x [ I ] [ J ] := C a l c u l a t e ( P X , P Y , Q X , Q Y , I ) ;

end; (* for J *) end; (* for I *) For J :- i to D I m F do F V e c t o r [ J ] := F C a l c u l a t e ( E v a l P o i n t s [ J ] [ 1 ] , E v a I P o i n t s [ J ] [ 2 ] , J); g r a m o n l y ( F V e c t o r , F M a t r i x , DimF, 0, N u m b e r O f N o d e s , 1.0, Gan~,aVector); PrintOut; w r i t e ( ' D o you want to run a n o t h e r test c a s e ? (Y/N}'); readln(U); writeln; if (U = 'N') or (U = 'n') t h e n a g a i n := false; end; (* w h i l e *) close(W); close(v); end. (* End of P r o g r a m *)

T h e f o l l o w i n g l i s t i n g is f r o m f i n d e v a l p o i n t s . p a s , and c o n t a i n s the p r o c e d u r e s n e c e s s a r y to c a l c u l a t e the c o o r d i n a t e s of all of the e v a l u a t i o n p o i n t s u s e d in the real v a r l a b l e b o u n d a r y e l e m e n t s p r o g r a m : ******************************************************************** * * *

Sharon Gee Dennis McOulgan Cindy Shrout

* * Program Units * * * Variables * * * *

-

-

Procedure Procedure Procedure Procedure

FindEval_Foints Down Up Equal

Ptr : index array point

i n t e g e r v a r l a b l e c o n t a l n i n g the of the l o c a t i o n in the c o o r d i n a t e of w h e r e the n e x t e v a l u a t i o n is.

Start

: i n t e g e r v a r i a b l e c o n t a i n l n g the

A study of solving two-dimensional potential problems first n o d e n u m b e r on the s i d e w h e r e e v a l u a t i o n p o i n t is to be placed. Finish : integer last n o d e n u m b e r evaluation point

the

v a r i a b l e c o n t a i n i n g the on the s i d e w h e r e the is to be placed.

NumEral8 : integer variable containing the n u m b e r of e v a l u a t i o n p o l n t s to be p l a c e d b e t w e e n the n o d e s on a b o u n d a r y . A c t u a l _ P t s : a r r a y c o n t a i n i n g the coordi n a t e s of the e v a l u a t i o n points. ********************************************************************* Procedure

FindEval_Polnts

* * * * * * * * * * * *

(vat Ptr: integer; Start, Finlsh, N u m E r a l 8 : integer: va~ A c t u a l _ P t s : e v a l t y p e );

var dl, d2 alltheway x_space, y _ s p a c e temp done

: : : : :

real; integer; =eel; integer; boolean~

********************************************************************* (* P R O C E D U R E UP (* . . . . . . . . . . . . (* T h i s p r o c e d u r e c a l c u l a t e s a c o o r d i n a t e of the e v a l u a t i o n (* p o i n t in the c a s e w h e r e the c o o r d i n a t e of the first n o d e (* is less t h a n the c o s r d i n a t e of the s e c o n d node. ********************************************************************* procedure

Up

amount xyl, a space vat s p o t var index

*) *) *) *) *)

: integer; : integer; : real; evaltype; : i n t e g e r );

vat m : integer; begin for m := 1 to a m o u n t do begln spot[index,a] := N C [ x y l ~ a ] i n d e x := i n d e x + l; end: end; (~ p r o c e d u r e up *)

+ (m'space)/

********************************************************************* (* P R O C E D U R E D O W N <* . . . . . . . . . . . . . [* T h z s p r o c e d u r e c a l c u l a t e s a c o o r d l n a t e of the e v a l u a t i o n (* p o i n t in the c a s e w h e r e the c o o r d i n a t e of the first n o d e (* is g r e a t e r t h a n the c o o r d i n a t e of the s e c o n d node. ********************************************************************* procedure Down ( amount : integer; xyl, a : integer; space : ~eal; vat spot : evaltype; vat index : i n t e g e r ); vat : integer; m begin for m := 1 to a m o u n t do begln spot[index,a] :- N C [ x y i , a ] - [m'space); index ;- index + I; end; end;

(* p r o c e d u r e

*) *> *) "I *)

d o w n *)

********************************************************************* (* P R O C E D U R E E Q U A L (* ............ (* T h i s p r o c e d u r e c a l c u l a t e s a c o o r d l n a t e of the e v a l u a t i o n (* point in the c a s e w h e r e the c o o r d l n a t e of the flrst n o d e (* is egual to the c o o r d i n a t e of the s e c o n d node. ********************************************************************* p r o c e d u r e Equal ( a m o u n t : integer; xyl, a : integer; var spot : evaltype; var index : i n t e g e r ); vat m : integer; begin for m := 1 to a m o u n t do begin spot[index,a] := N C [ X y I , a ] ;

*) *) "I *) *)

end;

i n d e x := i n d e x end; (* p r o c e d u r e egual

+ i; *)

********************************************************************* (* P R O C E D U R E R I N D E V A L _ P O I N T S *) (* . . . . . . . . . . . . . . . . . . . . . . . . . *) (* T h i s p r o c e d u r e c h e c k s the o r i e n t a t i o n of the two n o d e s *) (* b e t w e e n w h i c h the e v a l u a t i o n p o i n t ( s ) is (are) to be placed, *) (* a n d calls the p r e c e d i n t t h r e e p r o c e d u r e s as needed. *) ********************************************************************* begin d o n e := false; a l l t h e w a y := N u m E r a l S ; if ( ( [ N C [ s t a r t , l ] ) >= ( N C [ f i n l s h , l ] ) ) and ( ( R C [ s t a r t , 2 ] 1 >= ( N C [ f i n i s h , 2 ] ) ) ) then begin X _ s p a c e := ( N C [ s t a r t , l ] ~ N C [ f i n i s h , l ] ) / ( A l l T h e W a y + i); y _ s p a c e := ( N C [ s t a r t , 2 ] - N c [ f i n i s h , 2 ] ) / ( A l l T h e W a y + i); t e m p := ptr; if ( ( N C [ s t a r t , l ] ) • ( N C [ f i n i s h ~ l ] ) ) t h e n D o w n ( a l l t h e w a y , start, i, x_space, a c t u a l _ p t $ , pt~ )~ if ( ( N C [ s t a r t , l ] ) = ( N C [ f i n i s h , l ] ) ) t h e n Egual ( a l l t h e w a y , start, i, a c t u a l _ p t s , ptr ); if ( ( N O [ s t a r t , 2 ] ) • ( N C [ f i n i s h , 2 ] ) ) t h e n D o w n ( a l l t h e w a y , start, 2, y_space, a c t u a l _ p t s , t e m p ); if ( ( N C [ s t a r t , 2 ] ) = ( N C [ f i n i s h , 2 ] ) ) t h e n Egual ( a l ] t h e w a y , start, 2, a c t u a l _ p t s , t e m p ); d o n e := true; end; if [not done) t h e n begin if ( ( ( N C [ s t a r t , l ] ) <= ( N C [ f i n i s h , l ] ) ) and ( ( N O [ s t a r t , 2 ] ) >= ( N C [ f i n l s h , 2 ] ) ) ) then begin x _ s p a c e := ( N C [ f i n i s h , l ] - N C [ s t a r t , l ] ) / ( A l l T h e W a y + i); y _ s p a c e := ( N C [ s t a r t , 2 ] - N C [ f i n i s h , 2 ] ) / ( A l l T h e W a y + 1): temp := ptr; if ( ( N C [ s t a r t , l ] ) < ( N C [ f i n i s h , l ] ) ) t h e n up ( a l l t h e w a y , start, i, x_space, a c t u a l _ p t s , ptr ); if [ ( N C [ s t a r t , l ] ) = ( N C [ f i n i s h , l ] ) ) t h e n Equal ( a l l t h e w a y , start, i, a c t u a l _ p t s , ptr ); if ( ( N C [ s t a r t , 2 ] ) > ( N C [ f i n i s h , 2 ] ) ) t h e n D o w n ( a l l t h e w a y , start, 2, y_space, a c t u a l _ p t s , t e m p ); if ( ( N C [ s t a r t , 2 ] ) = ( N C [ f i n i s h , 2 ] ) ) t h e n Equal ( a l l t h e w a y , start, 2, a c t u a l _ p t s , t e m p ); d o n e := true; end; end; If (not done) t h e n begin if ( ( ( N C [ s t a r t , l ] ) <= ( N C [ f i n i s h , l ] ) ) and ( ( N C [ s t a r t , 2 ] ) <= ( N C [ f i n i s h , 2 ] ) ) 1 then begin x _ s p a c e := ( N C [ f i n i s h , l ] - N C [ s t a r t , l ] ) / ( A l l T h e W a y + i); y _ s p a c e :: ( N C [ f l n i s h , 2 ] - N C [ s t a r t , 2 ] ) / ( A l i T h e W a y + i); temp := ptr; if ( ( N C [ s t a r t , l ] ) < ( N C [ f i n i s h , l ] ) ) then Up ( a l l t h e w a y , start, i, x_space, a c t u a | _ p t s , ptr I; if ( ( R C [ s t a r t , l ] ) = ( N C [ f i n i s h , l ] ) ) t h e n Equal ( a l l t h e w a y , start, i, a c t u a l _ p t s , p t r ); if ( ( N C [ s t a r t , 2 ] ) < ( N C [ f i n i s h , 2 ] ) ) t h e n Up ( a l l t h e w a y , start, 2, y_space, a c t u a l _ p t s , t e m p )~ if ( ( N C [ s t a r t , 2 ] ) = ( N C [ f i n i s h , 2 ] ) ) t h e n Equal ( a l l t h e w a y , start, 2, a c t u a l _ p t s , t e m p ); d o n e := true; end; end; if (not donel t h e n begin if ( ( ( N C [ s t a r t , l ] ) >= ( N C [ f i n i s h ~ l ] ) ) and ( ( N C [ s t a r t , 2 ] ) <= ( N C [ f i n i s h , 2 ] ) ) ) t h e n begln x _ s p a c e := ( N C [ s t a r t , l ] - N C [ f i n i s h , l ] ) / ( A l ] T h e W a y + I); y _ s p a c e := ( N C [ f i n i s h , 2 ] - N C [ s t a r t , 2 ] ) / ( A l l T h e W a y + i); t e m p := ptr; if ( ( N C [ s t a r £ , l ] ) > ( N C [ f i n i s h , l ] ) ) t h e n D o w n ( a l l t h e w a y , start, i, x_space, a c t u a l _ p t s , ptr ); if ( ( N O [ s t a r t , l ] ) = ( R C [ f i n i s h , l ] ) ) t h e n Equal ( a l l t h e w a y , start, I, a c t u a l _ p t s , ptr ); if ( ( N C [ s t a r t , 2 ] ) < ( N C [ f i n i s h , 2 ] ) ) t h e n Up ( a l l t h e w a y , start, 2, y_space, a c t u a l _ p t s , t e m p ); if ( ( N C [ s t a r t , 2 ] ) = ( N C [ f i n i s h , 2 ] ) ) t h e n Egual ( a l l t h e w a y , start, 2, a c t u a l _ p t s , t e m p ); d o n e := true; end: end; end; (* p r o c e d u r e find eval p o i n t s *)

25