Algorithm for the discretised BIE

In this section we shall describe our algorithm to compute the integrals in the discretised variational equation in (3.100) efficiently. The method presented here will be used to solve the discretised integral equation for the crack opening displacement in conjunction with certain iterative solvers for matrix equations. Before we start, we have to define some words precisely. Boundary elements refer to small disjoint patches that cover the crack S. Each boundary element has local nodes on its boundary or in its interior. Local shape functions are defined in the interior of each element, and each of these local shape functions is associated with a local node where it takes the value of 1, while at other local nodes it is equal to 0. Some local nodes of different boundary elements may share a global position called a global node. The union of all the local shape functions that takes the value of 1 at a global node is called a basis function, or a global shape function associated with that global node. The combination of such global shape function and the union of the relevant boundary elements are called the global boundary element associated with the global node. We take Gauss points on each global boundary element to compute the outer integral on the right-hand side of (3.100).

We now move on to the description of the algorithm. We assume that the function $\phi_i$ is known on the crack, and is discretised as

$$\sum_J N^J(\mbox{\boldmath$ x $ }) \phi_i^J$$ (3.111)

where $N^J(\mbox{\boldmath$\space x $ })$ is the basis function associated with the global node x^J and $\phi_i^J$ is the nodal value of $\phi_i$ (See Fig.3.27). The global nodes are taken only in the interior of S so that the regularity condition in (3.43) is satisfied. Our objective is to compute a vector whose (3I+a-3)th components are equal to the right-hand side of (3.100) with $\psi_a$ replaced by N^I.

  

Figure 3.27: Piecewise linear shape function N^J and global node x^J

  

Figure 3.28: Global boundary element belongs to a cell or not

Steps 1-2. Same as the steps 1-2 in the algorithm described in chapter 2. In FM-GBIEM a global boundary element is considered to belong to a cell if the associated global node is in a cell (See Fig.3.28).

Step 3. Computation of the multipole moments (Upward)
First compute the multipole moments associated with leaves using (3.53) and (3.54). Here the multipole moments associated with a cell C mean the integrals in (3.53) and (3.54) with O taken as the centroid of C and $\phi_i$ understood as

$$\sum_{x^I\in C} N^I(\mbox{\boldmath$ x $ }) \phi_i^I.$$

Now consider a non-leaf cell C of level l. We compute the multipole moments associated with C by adding all the multipole moments of C's children after shifting the origin from the centroids of C's children to that of C via (3.57) and (3.58). We repeat this procedure tracing the tree structure of cells obtained in step 2 upward (decreasing l) until we reach level 2 cells.

Step 4. Computation of the coefficients of the local expansion (Downward)

We first compute the coefficients of the local expansion of cells of level 2 according to the definition using (3.60) and (3.61). After that, the coefficients of the local expansion of a level l cell C are computed recursively in 2 steps. Namely, we first add together the contributions of the forms in (3.60) and (3.61) from cells of the level l which are well-separated from C. This is followed by the addition of the coefficients of the local expansion of the parent of C with the origin shifted from the centroid of the parent to that of C via (3.64) or (3.65).

Step 5. Evaluation of the integral in (3.100)
We now compute the right-hand side of (3.100) with the test function $\psi_i$ replaced by a basis function N^I: this integral is denoted by a^I. Let C be a cell of the level l to which x^I belongs, and let C' be a cell of level l adjacent to C. We then compute the following sum if either C or C' is a leaf:

$$\sum_{J: x^J\in C'}\int_{S} N^I(\mbox{\boldmath$ x $ }) \ \... ...ath$ y $ }) N^J(\mbox{\boldmath$ y $ }) dS_{y} dS_x \phi_d^J.$$ (3.112)

This computation is carried out directly with the help of the regularised expression on the right-hand side of (3.104). When C is a leaf, we also compute the right-hand side of (3.106) using the coefficients of the local expansion associated with C for $L^{1,2}_{\cdot,N,M}$ and N^I for $\psi_i$. The sum of all the terms thus computed as one increases l from 2 to its maximum value is equal to the desired quantity, i.e., a^I.

In the actual computation the 4th and 5th steps are carried out at the same time in a single do loop with respect to l. It is seen that the computational cost for the above algorithm is of the order of N if one truncates the infinite sums in (3.60), etc. at a fixed number P. Notice, however, that the size of the boundary elements has to be fine enough so that the conditions of validity of various expansions such as the one stated after (3.105) are not violated.