■Assembly作業■
三角形要素による要素　通称"剛性"マトリクス[k]の計算方法が分かったところで、このマトリクス[k]をグローバル"剛性"マトリクス[K]に足し込む(assembly)方法を勉強しておきましょう。
まず、要素毎に[k]を計算します。これは、3×3のマトリクスで、ここではk_KLで表すことにします。この[k]を節点番号に注意を払いながらグローバル"剛性"マトリクス[K]に足し込めばよいことになります。Assemblyの実演については、 計算例をご覧下さい。ここでは、連立方程式の解法で紹介した対称Multi-Diagonal Matrixを[K]として使うことにします。
次に要素のK節点に対応したグローバル節点番号をiで表します。そして、要素のL節点に対応したグローバル節点番号からiを差し引きそして１を加算した値をjと置きます。ここで、注意したいことは、[K]の１列が[k]の対角要素に対応します。すると、k_KLはK_ijに足し込めば良いことになります。ただし、j>0という条件が有ります。
Fortranで書くと以下のようになります。ここにNEは要素数、IELは要素番号、NDは要素の節点番号、SKは[k]、RKは[K]のことです。下のプログラムでは、要素毎剛性マトリクスSK(I,J)を計算し、その後SK(I,J)をグローバル剛性マトリクスRK(I,J)に足し込んでいます。

C------ INTEGRATION OF STIFFNESS TERM
      DO IEL = 1 ,NE
      B(1,1) = YCOORD(NODEX(IEL,2)) - YCOORD(NODEX(IEL,3))
      B(1,2) = YCOORD(NODEX(IEL,3)) - YCOORD(NODEX(IEL,1))
      B(1,3) = YCOORD(NODEX(IEL,1)) - YCOORD(NODEX(IEL,2))
      B(2,1) = XCOORD(NODEX(IEL,3)) - XCOORD(NODEX(IEL,2))
      B(2,2) = XCOORD(NODEX(IEL,1)) - XCOORD(NODEX(IEL,3))
      B(2,3) = XCOORD(NODEX(IEL,2)) - XCOORD(NODEX(IEL,1))
      AREA = ( B(2,3)*B(1,2) - B(1,3)*B(2,2) ) / 2.D0
C-------B(1,J)=DN/DX----B(2,J)=DN/DY
      DO I = 1 , ND
      B(1,I) = B(1,I) / (2.D0*AREA)
      B(2,I) = B(2,I) / (2.D0*AREA)
      END DO
      DO I = 1 , ND
      DO J = 1 , ND
      SK(I,J) = ( (B(1,I)*EXX+B(2,I)*EXY) * B(1,J)
     *        +   (B(1,I)*EXY+B(2,I)*EYY) * B(2,J) )*AREA
      END DO
      END DO
C----ASSEMBLY OF SK(I,J) INTO RH(I,J)
      DO K = 1 , ND
      I = NODEX(IEL,K)
      DO L = 1 , ND
      J = NODEX(IEL,L) - I + 1
      IF ( J .GE. 1 ) RK (I,J) = RK(I,J) + SK(K,L)
      END DO
      END DO
      END DO