2
$\begingroup$

I am trying to assemble the global stiffness matrix and global force vector from the local stiffness matrix and local force vector using a function Forcestiffness Assembly as follows: where

FA[no of elements, total structural degrees of freedom] : Global force vector.

force[element degree of freedom] : local force vector.

iel : element number.

nodes[number of nodes in an element, total no of nodes] : Nodal connectivity matrix.

ndof : degress of freedom per node.

index[ element degree of freedom] : converts the element node number to global degree of freedom.

edof : number of degrees of freedom per element.

StiffnessE[ edof, edof]: Local Stiffness matrix.

nnel: number of node in an element.

KKA[no of elements, total structural degrees of freedom, total structural degrees of freedom] : Global Stiffness Matrix.

ForceStiffnessAssembly[FA_, force_, iel_, nodes_, ndof_, index_, StiffnessE_, edof_, nnel_, KKA_] := 
  Block[{ik, start, ij, ji, ii, jj},

                ik = 0;
    Do[
            start = (nodes[[iel, ij]] - 1)*ndof;
            Do[
                ik = ik + 1;
                index[[ik]] = start + ji;
                            , {ji, 1, ndof}];
                         , {ij, 1, nnel}];
           Do[
            ii = index[[ij]];
FA[[iel, ii]] = force[[ij]];
            Do[
                jj = index[[ji]];
                KKA[[iel, ii, jj]] = StiffnessE[[ij, ji]];
                      , {ji, 1, edof}];
                          , {ij, 1, edof}];

   Return[{FA, KKA}];
   ];

Finally iam add this way to get assembled global matrices

KKF = Sum[KKA[[i, All, All]], {i, 1, nelem}];
Fv = Sum[FA[[i, All]], {i, 1, nelem}];

But I am getting error as Set::shape: Lists {FA,KKA} and ForceStiffnessAssembly[<<1>>] are not the same shape. Please help

$\endgroup$
  • $\begingroup$ An example dataset that reproduced the error would be helpful... $\endgroup$ – Henrik Schumacher Dec 9 '19 at 7:00
  • 1
    $\begingroup$ Also, in your example code, the function ForceStiffnessAssembly is not called at all. So this does not at all reproduce the error. But my guess: ForceStiffnessAssembly[<<1>>] indicates that you call ForceStiffnessAssembly with only one argument for which the expression stays unevaluated (and has length different from 2). Something like {a, b} = f[c] with f undefined. $\endgroup$ – Henrik Schumacher Dec 9 '19 at 7:27
8
$\begingroup$

Here is a simple code to do that. Set up a basic mesh:

meshCoords = N[{{0, 0}, {1, 0}, {1, 1}, {0, 1}}];
incidents = {{1, 2, 4}, {2, 3, 4}};

We assume a single equation:

(* single equation *)
dof = Length[meshCoords];
elementSize = Last[Dimensions[incidents]];

Now, we find all positions where this mesh will introduce values in the system matrix:

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

Create some symbolic FEM entries; you'd need to change this to suite you needs.

(* symbolic elements for fun*)
elementValues = Array[#, {elementSize, elementSize}] & /@ {a, b};
Dimensions[elementValues]
{2, 3, 3}

So we have two triangle elements of 3 rows and 3 columns each. Next, we assemble those. This will create a sparse array where duplicate positions given in pos will be added up.

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

Do the assembly:

stiffness = matrixAssembly[elementValues, pos, dof] 

Look at the symbolic stiffness matrix:

MatrixForm[stiffness]

enter image description here

MatrixPlot[ stiffness]

I have given a talk about this a few years back; you can still find the talk here. You might also be interested in the FEM Programming tutorial.

$\endgroup$
  • $\begingroup$ If I have a rectangular mesh with quad elements (4 noded element) 65 in number, number of nodes = 84. At each node I have 5 degree of freedom. Total system dof =420. meshcoords={84,2}, incidents={65,4}. Then can you explain in this case. $\endgroup$ – S.B.MD.Khaja Moinuddin Dec 11 '19 at 13:43
  • $\begingroup$ @S.B.MD.KhajaMoinuddin, what exactly is the problem? Start with changing the above to use quad elements. If you want some thing more specific you'd need to improve your question by adding that information - as code, not as text. $\endgroup$ – user21 Dec 11 '19 at 14:52

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