sci.math.research

SUNA23, Substructuring
======================

This article is essentially an elaboration of some theory which was developed
in the following (IMHO exellent) monograph:

       F.J. Peters, "Sparse matrices and Substructures",
       Mathematical Centre Tracts 119, Amsterdam 1980.

(Disclaimer: in fact, below is all I could manage to comprehend of it. Shame;-)

Suppose that part of a (Finite Element) problem can be viewed upon in such a
way that there exist a distinction between so-called "internal" and "external"
unknowns. Unknowns are called internal if their values are uniquely determined
by the values of external unknowns. Therefore internal variables can always be
expressed into external variables. Hence they can be _eliminated_, leaving us
hopefully with a smaller problem, containing only external unknowns. In Finite
Element contexts, people often talk about super-elements (: bottom-up approach)
or sub-structures (: top-down approach). Whether substructuring is feasable can
be judged best by examining the problem at hand more closely. We will meet an
example in the planned follow-up article.

Back to math now. Mark the internal unknowns with (i) and the external unknowns
with (e). The equations for an element can then be decomposed in matrix form as
follows:

     | Qii   Qie | | Wi |   | Fi |
     |           | |    | = |    |
     | Qei   Qee | | We |   | Fe |

Here: Q = element matrix-blocks, W = vector of unknowns, F = load vector.

We are working towards elimination of the internal nodes:

                                           -1
  Qii . Wi + Qie . We = Fi   -->   Wi = Qii  . (Fi - Qie . We)

  Qei . Wi + Qee . We = Fe   -->

           -1                -1
  Qei . Qii  . Fi - Qei . Qii  . Qie . We + Qee . We = Fe

                               -1
Define:    Qr = Qee - Qei . Qii   . Qie
                                              --->   Qr . We = Fr
                              -1
           Fr = Fe - Qei . Qii  . Fi          Qr = "reduced" element matrix
                                              Fr = "reduced" load vector

In the sequel, it will be assumed that all element matrices are symmetric and
positive definite, which is the case with least squares finite element methods.

Algorithm:  1. Invert  Qii  in place  -->  Qii ;
---------   2. Calculate  Qei . Qii   -->  Qei ;
            3. Calculate  Qee - Qei . Qie   -->  Qee : reduced matrix ;
            4. Calculate  Fe = Fe - Qei.Fi  -->  Fe  : reduced vector .

Then the global system of equations will be assembled from reduced elements,
and hopefully it can be solved. Having done that, the internal unknowns can be
found by backsubstitution:
                                   -1
                           Wi = Qii  . (Fi - Qie . We)

Below is the basic (-: yes!) "subroutine" that implements the whole procedure.

1720 REM   SUBSTRUCTURING
1730 REM   --------------
1740 REM   M = NUMBER OF INTERNAL NODES   | II IE |
1750 REM   N = TOTAL NUMBER OF NODES      | EI EE |
1760 REM   E = SYMMETRIC POSITIVE DEFINITE MATRIX
1770 REM   R = LOAD VECTOR
1780 REM   W = UNKNOWNS
1781 REM
1785 REM   BACK = SWITCH FOR BACKSUBSTITUTION
1790 REM
1800 REM ... Invert E(II) ...
1810 FOR K=1 TO M : PE=E(K,K)
1820 IF PE=0 THEN PRINT "PIVOT 0 AT";K : END
1830 E(K,K)=1 : FOR I=1 TO M
1840 E(I,K)=E(I,K)/PE : IF I=K THEN GOTO 1875
1850 FOR J=1 TO M : IF J=K THEN GOTO 1870
1860 E(I,J)=E(I,J)-E(I,K)*E(K,J)
1870 NEXT J
1875 NEXT I
1880 FOR J=1 TO M : E(K,J)=-E(K,J)/PE : NEXT J
1890 E(K,K)=1/PE : NEXT K
1900 REM
1910 IF BACK=1 THEN GOTO 2080
1920 REM
1930 REM ... E(EI) = E(IE) * E(II)^-1 ...
1940 FOR I=M+1 TO N : FOR J=1 TO M
1950 H=0 : FOR K=1 TO M : H=H+E(K,I)*E(K,J) : NEXT K
1960 E(I,J)=H : NEXT J : NEXT I
1970 REM
1980 REM ... E(EE) = E(EE) - E(EI)*E(IE) ...
1990 FOR I=M+1 TO N : FOR J=M+1 TO N
2000 H=0 : FOR K=1 TO M : H=H+E(I,K)*E(K,J) : NEXT K
2010 E(I,J)=E(I,J)-H : NEXT J : NEXT I
2020 REM
2030 REM ... R(E) = R(E) - E(EI)*R(I) ...
2040 FOR I=M+1 TO N
2050 H=0 : FOR K=1 TO M : H=H+E(I,K)*R(K) : NEXT K
2060 R(I)=R(I)-H : NEXT I
2070 RETURN
2080 REM ... Backsubstitution ...
2090 FOR I=1 TO M
2100 H=0 : FOR K=M+1 TO N : H=H+E(I,K)*W(K) : NEXT K
2110 R(I)=R(I)-H : NEXT I
2120 FOR I=1 TO M
2130 H=0 : FOR K=1 TO M : H=H+E(I,K)*R(K) : NEXT K
2140 W(I)=H : NEXT I : RETURN

To be continued ...
-
* Han de Bruijn; Applications&Graphics | "A little bit of Physics * No
* TUD Computing Centre; P.O. Box 354   | would be NO idleness in  * Oil
* 2600 AJ  Delft; The Netherlands.     | Mathematics" (HdB).      * for
* E-mail: Han.deBruijn@RC.TUDelft.NL --| Fax: +31 15 78 37 87 ----* Blood