Creating overlapped block matrices

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?

• Is there a pattern to this? With the b23 where it is and not a21 on the left side of that row I can't think of one. – C. E. Aug 19 '13 at 8:52
• Anon, there are stiffness matrices. There is a coupling between elements and where they join, or overlap, the stiffness is added. But this all depends on the elements and the structure and can get very complicated. – Nasser Aug 19 '13 at 10:00
• Viet Tran, you haven't responded to Anon's comment or any of the answerers. It seems like the 0 in position {2,1} was a typo (as everyone seems to agree) and I've edited the question to change it to a21, so that the answers make sense. If you disagree, please leave a comment, else I'll assume that you are in agreement with this change. – rm -rf Dec 31 '13 at 1:18

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},
SystemSetSystemOptions[
"SparseArrayOptions" -> {"TreatRepeatedEntries" -> 1}];
matrix = SparseArray[ pos -> Flatten[ values], dim];
SystemSetSystemOptions[
"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.

• Alternatively, one can use InternalWithLocalSettings[]: matrixAssembly[values_, pos_, dim_] := With[{spopt = SystemOptions["SparseArrayOptions"]}, InternalWithLocalSettings[SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> 1}], SparseArray[pos -> Flatten[values], dim], SetSystemOptions[spopt]]] – J. M. will be back soon Jul 5 '16 at 12:48

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[] - 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.

• This has a21 where there is a zero in the OP's example. I think that might be a typo on the OP's part though. – C. E. Aug 19 '13 at 8:51
• @Anon ups, you are right. But I think it is a mistake in question :p – Kuba Aug 19 '13 at 8:51
• @Cogicero sorry, haven't got pinged earlier. What about x[[2,3]]? – Kuba Mar 7 '16 at 20:25
• @Cogicero Unless I'm mistaken you are using wrong syntax, you can't take element from A with A[2,3] but you have to use e.g. Part: A[[2,3]]. – Kuba Mar 8 '16 at 7:18
• @Kuba Each matrix has a value before you run the above code, but after running your code I still get the picture above with elements like a[1,1], etc. I can't use the global Matrix that I end up with, because instead of a[1,1] I should have "1" since that's row 1, col1 of a 3*3 unit matrix. Instead of a[1,2] appearing in my final answer it should be 0 since that's row 1 col 2 of a 3*3 unit matrix. Mathematica still treats the matrices as if they are arbitrary and unknown. I really hope you understand my problem now. It's a bit frustrating, and I am not sure how to GET the actual Parts to show. – Cogicero Mar 8 '16 at 12:01

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.

(* 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: • I posted similar answer but @Anon has pointed out that it does not fit the question. Take a look at position {2, 1} of OP's result. – Kuba Aug 19 '13 at 10:02
• @Kuba. I see. Oh well. I just guess the OP is either wrong (made typo), or they should then say what does this pattern represents. I assumed it is a global stiffness matrix. If not, then I will also delete this answer. – Nasser Aug 19 '13 at 10:04
• I think you are right, I think so. But we can't be sure :) – Kuba Aug 19 '13 at 10:08
• @Kuba I also read at the top OP saying FEM assignment using Mathematica. and that is why I thought right away it is stiffness matrix. (FEM== finite elements) – Nasser Aug 19 '13 at 10:10
• Good point. I could have checked this. – Kuba Aug 19 '13 at 10:13

Using SparseArrays

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

{sEK1, sEK2} =
SparseArray[#, {3, 3},
0] &@(Thread[#[[2 ;;]] -> #[][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}}