Creating overlapped block matrices

Question

I'm doing an FEM assignment using Mathematica.

(EK1 = {{a11, a12}, {a21, a22}}) // MatrixForm
(EK2 = {{b22, b23}, {b32, b33}}) // MatrixForm

I don't know how the best way to create new matrix like this:

(K = {{a11, a12, 0}, {a21, a22 + b22, b23}, {0, b32, b33}}) // MatrixForm

MATLAB does it quite easily but how about Mathematica?

Answer 1

I gave a presentation on one (efficient) way to do it at the 2009 tech conference.

The essence is this function:

matrixAssembly[ values_, pos_, dim_] := Block[{matrix, p},
  System`SetSystemOptions[ 
   "SparseArrayOptions" -> {"TreatRepeatedEntries" -> 1}];
  matrix = SparseArray[ pos -> Flatten[ values], dim];
  System`SetSystemOptions[ 
   "SparseArrayOptions" -> {"TreatRepeatedEntries" -> 0}];
  Return[ matrix]]

The option "SparseArrayOptions" -> {"TreatRepeatedEntries" -> 1 will add up elements at the same position.

pos can be computed like this:

pos = Flatten[ Map[ Outer[ List, #, #] &, incidents], 2];

Where, as an example, incidents of element 1 and 3 look like this:

incidents[[{1, 3}]]
{{23, 82, 64, 83}, {83, 82, 67, 66}}

In the above mentioned presentation you'll also find a method how to apply Dirichlet conditions, time integrating the system and perform some model order reduction.

Also note that since version 10.0 Mathematica has a Finite Element Programming tutorial that shows how to use the new in V10 FEM building blocks on a low level.

Answer 2

You can use ArrayPad:

ArrayPad[EK2, {1, 0}] + ArrayPad[EK1, {0, 1}]

and general approach (Edit: overlap specification added):

n = {3, 3};
arrays = Array[#, n] & /@ {a, b, c, d, e}; (*arrays to work with*)

app[a1_, a2_, overlap_: 1] := With[{dim = Dimensions@a1},
                                   ArrayPad[a1, {0, n[[1]] - overlap}] + 
                                   ArrayPad[a2, Transpose@{dim - overlap, {0, 0}}]];

In V10+ once can use Fold[app, arrays] instead of Fold[app, First@arrays,Rest@arrays]

Fold[app, arrays] // MatrixForm
Fold[app[##, 2] &, arrays] // MatrixForm
Fold[app[##, 3] &, arrays] // MatrixForm

It is not so much general, only for constant n$x$n set but I don't know what should be an expected result for set of arrays of different dimensions.

Answer 3

I wrote a small function once to build block-wise diagonal matrices. The zero matrix is just 0 and identity is 1. You could use it in this case like this:

diagonalize[list_] := ArrayFlatten@(DiagonalMatrix[Array[x, Length@list]] /. 
    Table[x[i] -> list[[i]], {i, Length@list}]);
diagonalize[{EK1, 0}] + diagonalize[{0, EK2}] // MatrixForm

Oh, this is assuming that the second zero in your K is actually a21. Otherwise I don't see a pattern here.

Answer 4

(* I am assuming this is meant to be a global stiffness matrix construction question and OP made type for element (2,1) *)

There are many ways to build global stiffness matrix. For your case, here is a quick hack. But if you simply google the topic of building global stiffness matrix from elements stiffness matrices, you'll find many many methods on it in finite elements books.

Notice, I gave element names inside the local stiffness as b11,b12, etc.... different than what you had to be consistent with the first matrix.

EK1 = {{a11, a12}, {a21, a22}}; (*local stiffness matrix*)
EK2 = {{b11, b12}, {b21, b22}}; (*local stiffness matrix*)
n = 2; (*size of local stiffness matrix *)
locals = {EK1, EK2}; (*local stiffness matrices*)
m = Length[locals]; (*number of elements*)
r = (n m) - (m - 1); (*size of global stiffness matrix*)
K1 = Table[0, {r}, {r}];
z = 1;

Do[ K1[[z ;; z + n - 1, z ;; z + n - 1]] = locals[[i]];
    If[i > 1, K1[[z, z]] += locals[[i - 1]][[-1, -1]]];
    z = z + n - 1,
  {i, 1, m} ];

For example, for 5 elements, you'll get this:

Answer 5

Using SparseArrays

EK1 = {{a11, a12}, {a21, a22}};
EK2 = {{b22, b23}, {b32, b33}};

{sEK1, sEK2} = 
SparseArray[#, {3, 3}, 
   0] &@(Thread[#[[2 ;;]] -> #[[1]][Sequence @@ #[[2 ;;]]], 
       2] &@(ToExpression@
       Characters@ToString[#]) & /@ (Flatten@#)) & /@ {EK1, EK2};

(sEK1 + sEK2) // Normal

{{a[1, 1], a[1, 2], 0}, {a[2, 1], a[2, 2] + b[2, 2], b[2, 3]}, {0, b[3, 2], b[3, 3]}}

Or, more simply:

EK1alt = SparseArray[{{1, 1} -> a11, {1, 2} -> a12, {2, 1} -> a21, {2, 2} -> a22}, {3, 3}, 0];
EK2alt = SparseArray[{{2, 2} -> b22, {2, 3} -> b23, {3, 2} -> b32, {3, 3} -> b33}, {3, 3}, 0];

(EK1alt + EK2alt) // Normal

{{a11, a12, 0}, {a21, a22 + b22, b23}, {0, b32, b33}}