1

Edit:

Okay, my question was viewed quite a few times but not answered so I assume I phrased/presented in a confusing way. I'll try to rephrase my problem to make it more accessible because luckily I figured out what caused my function to be slow as well as a rather simple fix and this might be useful for anyone that runs into a similar problem in the future.

Here's the slow code:

(*//////*)
(*Assembly of the global [K]-Matrix*)
(*//////*)

(*"empty" [K]-Matrix*)
kGlobal=SparseArray[{},{nodes,nodes}];  

(*Function that adds one element "e" to the global[K]-Matrix*)
  kAssembly[e_]:=(
   n1=conMatrix[[e,1]];
   n2=conMatrix[[e,2]];
   n3=conMatrix[[e,3]];

   kGlobal[[n1,n1]]=kGlobal[[n1,n1]]+kElement[[e,1,1]];
   kGlobal[[n1,n2]]=kGlobal[[n1,n2]]+kElement[[e,2,1]];
   kGlobal[[n1,n3]]=kGlobal[[n1,n3]]+kElement[[e,3,1]];

   kGlobal[[n2,n1]]=kGlobal[[n2,n1]]+kElement[[e,2,1]];
   kGlobal[[n2,n2]]=kGlobal[[n2,n2]]+kElement[[e,2,2]];
   kGlobal[[n2,n3]]=kGlobal[[n2,n3]]+kElement[[e,3,2]]; 

   kGlobal[[n3,n1]]=kGlobal[[n3,n1]]+kElement[[e,3,1]];
   kGlobal[[n3,n2]]=kGlobal[[n3,n2]]+kElement[[e,3,2]];
   kGlobal[[n3,n3]]=kGlobal[[n3,n3]]+kElement[[e,3,3]];
                  );

(*Executing the function for every element to assemble the entire global [K]-matrix*)   
Do[kAssembly[e],{e,1,elements}];

So what am I doing here:

I have a created a mesh beforehand and this mesh is built by nodes that are connected into Triangluar elements. Those connections are stored in the connectivity matrix conMatrix. Each row of the connectivity matrix has 3 values - the global numbers of its nodes, the first global node represents the elements first node the 2nd the 2nd node and so on. it looks something like this

$$conMatrix=\begin{bmatrix}1&3&2\\3&4&2\\3&5&4\\&...&\end{bmatrix}$$

Also beforehand I calculated the Stiffness - so called [K]-Matrix for each element, that's kElement. It represents relations between the nodes (for example k12 is the relation between node 1&2 in that element). kElementconsist out of the symmetric [K]-matrices of each element. Because it's symmetric, you can save 3 calculations per element and I didn't calculate it entirely:

$$kElement=\begin{bmatrix} \begin{bmatrix}k11&&\\k21&k22&\\k31&k32&k33\end{bmatrix}_{Element1},&\begin{bmatrix}k11&&\\k21&k22&\\k31&k32&k33\end{bmatrix}_{Element2},&...\end{bmatrix}$$

During the so called Assembly, the relations between each node on a global scale get described. In the conMatrix Example I gave, it is apparent, that global node 3 is part of Elements 1,2,3. That means that global $K_{11} =K_{22,Elem.1}+K_{11,Elem.2}+K_{11,Elem.3} $. Same goes for every other connection - for example elements can also share lines.

This assembly process is what kAssembly is doing. The function first extracts the elements nodes out of conMatrix, then it adds the values at the right position in kGlobal which is a SparseArray. The function delivers the correct result and works reasonably fast in small meshes. However the SparseArray hits dimensions > 20,000 x 20,000 it starts to become incredibly slow. 25,000 x 25,000 (with 176,000 non-zero values) took 275 seconds (on a really fast computer).

After some tinkering with SparseArray I realized that if you gradually add a lot of nonzero values to a SparseArray one at a time every following addition to the SparseArray becomes slower and slower. I did a little test to emphasize that by continously redcording the time after another 500 values were added to a SpraseArray - you can see the time increases slowly but in a exponential fashion and then there's some sort of bend (not sure what's going on):

enter image description here

p4di
  • 155
  • 8
  • 2
    Hard to say anything with incomplete code. We need all the definitions of the functions and variables involved in order to run your code. – MarcoB Jul 03 '18 at 20:00
  • there are pastebin links to the complete code (step 1: mesh generation, step 2: rest) in the original post. I'll try to edit the rest of my code in though now. – p4di Jul 04 '18 at 08:16

1 Answers1

0

To put it short:

it was the way I assembled my SparseArray. If you gradually add a lot of nonzero values to a SparseArray every further addition to the SparseArray becomes slower and slower. This is also described here.

Longer answer:

So it would be much better to Assemble the SparseArrayin one single step in a {node,node} -> value fashion. my problem though was that I would have had the same {node,node} expressions with different values (that need to be summed up) associated to them. SparseArray however ignores duplicate mentions and simply takes the first. So in a case {3,3}->1,{3,3}->3 it would add 1 at position {3,3} instead of the desired 4.

In Mr.Wizards answer to another problem, I've found a solution to this problem though.

The option SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> 1}]; does automatically sum up duplicate entries.

With that in mind, I built 2 functions to create the global node "coordinates"/indexes as well as the associated values for every element:

index[e_]:=(
            {n1,n2,n3}=conMatrix[[e,1;;3]];
            {{n1,n1},{n1,n2},{n1,n3},{n2,n1},{n2,n2},{n2,n3},{n3,n1},{n3,n2},{n3,n3}})    

values[e_]:=(
              k11=kElement[[e,1,1]];
             {k21,k22}=kElement[[e,2,1;;2]];
             {k31,k32,k33}=kElement[[e,3,1;;3]];
              k12=k21;k13=k31;k23=k32;
             {k11,k12,k13, k21,k22,k23, k31,k32,k33})

This allows me to apply those functions to every Element with Tables and then flatten them into a list:

ind=Flatten[Table[index[e],{e,1,elements}],1];]
val=Flatten[Table[values[e],{e,1,elements}],1];]

with all values and indexes known, all that's left is building the SparseArray:

SetSystemOptions["SparseArrayOptions"->{"TreatRepeatedEntries"->1}];
kGlobal=SparseArray[ind->val];]

The lists ind and valtook 0.2 & 0.3 seconds the SparseArray took 0.08 second. Significant speed-up compared to the 275 seconds before.

p4di
  • 155
  • 8