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 socalled "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 superelements (: bottomup approach)
or substructures (: topdown 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 followup 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 matrixblocks, 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
* Email: Han.deBruijn@RC.TUDelft.NL  Fax: +31 15 78 37 87 * Blood