Calculating Tetrahedral Stiffness Matrix for FEM

I'm writing some FEM simulations and have the mesh sorted.

I'm now putting the stiffness matrix together and feel that this is such a generic task that surely it's already been done - but after a quick google and look through the file exchange it appears not! There are a few things that mention trusses and things but nothing that takes an obvious input and gives a straightforward stiffness matrix.

Does anyone know of any good code that takes a mesh tetrahedralisation matrix T(Ne,4) and set of vertices V(Nv,3) and produces the linear stiffness matrix K(Nv,Nv)? Obviously higher orders would also be great but linear is enough to get started.

[Jason Nicholson]

I am not aware of any automatic global stiffness matrix generator. I think you have to role your own. That is all I have. Maybe, someone who studies this will comment. However, most FEA code writer don't work in MATLAB except for high level interface. FEA is written in C and C++.

[Nathan Welch]

Thanks - yeah this seems to be the standard method. Also quite a nice intro to FEM to work out how to do this I suppose!

Answer 1

[Nathan Welch]

As there doesn't seem to be any obvious sources, I've created my own. It calculates the stiffness matrix for linear elements of a tetrahedral mesh. As with all online code - use at your own risk. If anyone can describe how to better calculate this, esp how to expand to 2nd order or (the dream) arbitrary order, I'd still be very much interested to hear it.

function K = LinearStiffnessMatrix3(nodes, elements)

%get sizes

Nn = size(nodes,1);

Ne = size(elements,1);

%create sparse matrix by storing (i,j) indices and the values

inds = 0;

Si = zeros(16*Ne,1); %ith indices of stiffness matrix

Sj = zeros(16*Ne,1); %jth indices of stiffness matrix

Sk = zeros(16*Ne,1); %value of (i,j) element of stiffness matrix

%4 indices vs 4 indices

ii = 1:4; ji = 1:4; [ii, ji] = ndgrid(ii,ji);

ii = ii(:); ji = ji(:);

% Loop through every element in the mesh - would be nice to do this without

% looping - or this maybe why it's best to use a compiled language

for i = 1:Ne

eleNodes = elements(i,1:4);

%uses simple integration of linear terms method

B = [nodes(eleNodes,:)'; [1, 1, 1, 1]];

C = inv(B);

ve = abs(det(B));

Kl2 = C(:,1:3)*C(:,1:3)'.*abs(ve)/6;

%store all 16 terms

inds = inds(end) + (1:16);

Si(inds) = eleNodes(ii);

Sj(inds) = eleNodes(ji);

Sk(inds) = Kl2(:);

end

%Build stiffness matrix - MATLABs sparse operator naturally adds together

%all terms with the same i and j - i.e. all contributions to the ith node

%from nodes j

K = sparse(Si,Sj,Sk,Nn,Nn);

end

[Biswabhanu Puhan]

Thank you Mr. Welch, did you find the formulation for 2nd order element (C3D10 in abaqus)? Kindly please let me know.