5
$\begingroup$

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?

$\endgroup$
  • 2
    $\begingroup$ 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. $\endgroup$ – C. E. Aug 19 '13 at 8:52
  • $\begingroup$ 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. $\endgroup$ – Nasser Aug 19 '13 at 10:00
  • $\begingroup$ 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. $\endgroup$ – rm -rf Dec 31 '13 at 1:18
6
$\begingroup$

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.

$\endgroup$
  • 2
    $\begingroup$ Alternatively, one can use Internal`WithLocalSettings[]: matrixAssembly[values_, pos_, dim_] := With[{spopt = SystemOptions["SparseArrayOptions"]}, Internal`WithLocalSettings[SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> 1}], SparseArray[pos -> Flatten[values], dim], SetSystemOptions[spopt]]] $\endgroup$ – J. M. will be back soon Jul 5 '16 at 12:48
6
+50
$\begingroup$

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

enter image description here enter image description here enter image description here

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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – C. E. Aug 19 '13 at 8:51
  • $\begingroup$ @Anon ups, you are right. But I think it is a mistake in question :p $\endgroup$ – Kuba Aug 19 '13 at 8:51
  • $\begingroup$ @Cogicero sorry, haven't got pinged earlier. What about x[[2,3]]? $\endgroup$ – Kuba Mar 7 '16 at 20:25
  • $\begingroup$ @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]]. $\endgroup$ – Kuba Mar 8 '16 at 7:18
  • $\begingroup$ @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. $\endgroup$ – Cogicero Mar 8 '16 at 12:01
3
$\begingroup$

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.

$\endgroup$
2
$\begingroup$

(* 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} ];

Mathematica graphics

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

Mathematica graphics

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Kuba Aug 19 '13 at 10:02
  • $\begingroup$ @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. $\endgroup$ – Nasser Aug 19 '13 at 10:04
  • $\begingroup$ I think you are right, I think so. But we can't be sure :) $\endgroup$ – Kuba Aug 19 '13 at 10:08
  • 1
    $\begingroup$ @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) $\endgroup$ – Nasser Aug 19 '13 at 10:10
  • $\begingroup$ Good point. I could have checked this. $\endgroup$ – Kuba Aug 19 '13 at 10:13
0
$\begingroup$

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}}

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.