next up previous contents
Next: Algorithm Outline and Discussion Up: Numerical Solution Previous: Numerical Solution

System Matrix Assembling and the Skyline Storage Scheme


Since the nodes of the dynamic fish model are not fully connected by viscoelastic units, the effective stiffness matrix tex2html_wrap_inline2566 and, consequently (since tex2html_wrap_inline2514 is diagonal), the system matrix tex2html_wrap_inline2578 , are sparse matrices. Sparsity can be exploited to increase the computational efficiency of the factorization algorithms [Press et al.1986]. We have employed the Skyline Storage Scheme for storing tex2html_wrap_inline2578 , which is popular in numerical methods for finite element analysis [Bathe and Wilson1976]. This storage scheme is most efficient for storing sparse matrices that are symmetric and banded. We will now explain our strategy in more detail.

The tex2html_wrap_inline2528 system matrix tex2html_wrap_inline2578 and tex2html_wrap_inline2540 effective load matrix tex2html_wrap_inline2614 are computed according to Eq. (gif) and Eq. (gif), respectively. To assemble tex2html_wrap_inline2578 , we first assemble the instantaneous, internal effective viscoelastic stiffnesses tex2html_wrap_inline2618 of the system (see Eq. (gif)). Denoting the entries of tex2html_wrap_inline2578 as tex2html_wrap_inline2622 and the rows of tex2html_wrap_inline2614 as tex2html_wrap_inline2626 , the following procedure assembles tex2html_wrap_inline2578 and tex2html_wrap_inline2614 :

  1. Initialize tex2html_wrap_inline2632 and tex2html_wrap_inline2634 .
  2. For each viscoelastic unit tex2html_wrap_inline2344 , assemble the tex2html_wrap_inline2638 term of Eq. (gif) into tex2html_wrap_inline2578 :


    and assemble the tex2html_wrap_inline2642 term of Eq. (gif) into tex2html_wrap_inline2644 :


  3. Assemble the remaining terms of tex2html_wrap_inline2578 and tex2html_wrap_inline2614 . For tex2html_wrap_inline2650 :


Fig. gif shows the structure of the tex2html_wrap_inline2528 matrix tex2html_wrap_inline2578 which depends solely on the connectivity of the n=23 nodes in the biomechanical fish model. Note that the matrix is symmetric, sparse and banded. Bandedness means that all off-diagonal entries beyond the bandwidth of the matrix are zero. Because tex2html_wrap_inline2578 is symmetric, we can state the bandedness condition as


where tex2html_wrap_inline2662 is the bandwidth of tex2html_wrap_inline2578 and tex2html_wrap_inline2666 is called the half-bandwidth of tex2html_wrap_inline2578 . For example, if tex2html_wrap_inline2670 , then tex2html_wrap_inline2578 is a diagonal matrix. For our dynamic fish model, tex2html_wrap_inline2674 (see Fig. gif).


The skyline storage scheme exploits the fact that the many zeros outside the skyline of tex2html_wrap_inline2578 are unnecessary in the factorization solution of Eq. gif. Hence, an effective way of storing the symmetric tex2html_wrap_inline2578 would be to store only the upper triangular half of the matrix below the skyline, including the diagonal.gif The skyline storage scheme stores tex2html_wrap_inline2578 into a one-dimensional array tex2html_wrap_inline2688 . In addition, it uses an auxiliary array to properly address the entries of tex2html_wrap_inline2578 stored in tex2html_wrap_inline2688 .

We explain the matrix storage and addressing scheme via a small example system matrix tex2html_wrap_inline2694 of order 4 shown in Fig. gif. Let tex2html_wrap_inline2696 be the row number of the first nonzero entry in column i of tex2html_wrap_inline2694 . The variables tex2html_wrap_inline2696 , i=0, 1, 2, ..., define the skyline of the matrix and tex2html_wrap_inline2706 are the column heights. Furthermore, the half-bandwidth of the matrix, tex2html_wrap_inline2666 , equals tex2html_wrap_inline2710 . With the column heights defined, we can now store all entries of tex2html_wrap_inline2694 below its skyline consecutively in tex2html_wrap_inline2688 . Fig. gif shows how the entries of tex2html_wrap_inline2694 are stored in tex2html_wrap_inline2688 . In addition to tex2html_wrap_inline2688 , we also define an array MAX tex2html_wrap_inline2688 , which stores the indexes of the diagonal entries of tex2html_wrap_inline2694 , tex2html_wrap_inline2726 , in tex2html_wrap_inline2688 . Note that MAX tex2html_wrap_inline2688 (i) is equal to the sum of column heights up to the (i-1)th column plus i. Hence, the number of stored entries in the ith column of tex2html_wrap_inline2694 is equal to MAX tex2html_wrap_inline2688 (i+1) - MAX tex2html_wrap_inline2688 (i) and the entry indexes are MAX tex2html_wrap_inline2688 (i), MAX tex2html_wrap_inline2688 (i)+1, MAX tex2html_wrap_inline2688 (i)+2, ..., MAX tex2html_wrap_inline2688 (i+1) - 1. Therefore, storing tex2html_wrap_inline2694 in the form of tex2html_wrap_inline2688 and MAX tex2html_wrap_inline2688 allows the entries of tex2html_wrap_inline2694 to be easily addressed.


A factorization solution of Eq. gif using skyline storage is substantially more efficient compared to a Gaussian elimination solution which has tex2html_wrap_inline2780 complexity in terms of the number of operations required (an operation is defined as a multiplication and an addition), where n is the order of tex2html_wrap_inline2578 . For factorization with skyline storage, the number of operations required is tex2html_wrap_inline2786 in the worst case of constant column heights (i.e. tex2html_wrap_inline2788 for all i). Clearly, the larger the difference between n and tex2html_wrap_inline2666 , the greater the advantage. In simulating the dynamics of the artificial fish, we have n=23 and tex2html_wrap_inline2674 which represents a significant improvement in computational efficiency. Note that the column heights tex2html_wrap_inline2800 of a system matrix tex2html_wrap_inline2578 and, hence, the effectiveness of the storage scheme, depends on how the nodes in the discrete dynamic model are numbered. For example, using our current numbering scheme, tex2html_wrap_inline2674 (Fig. gif), however if we number node 22 as node 1 (referring to Fig. gif), tex2html_wrap_inline2806 . Thus in general, one should number the nodes in such a way that all connections are as local as possible. We can frequently determine a reasonable nodal point numbering by inspection. However, this numbering may not be particularly easy to generate in certain cases [Bathe and Wilson1976].gif

next up previous contents
Next: Algorithm Outline and Discussion Up: Numerical Solution Previous: Numerical Solution
Xiaoyuan TuJanuary 1996