# Custom element Stiffness Matrix

To build a custom stiffness matrix for a particular FEM element, just derive your own element class from the relevant standard one and override the CalcStiffnessMatrix() method. You need only to fill the upper triangle of the K matrix, the StiffnessComputation() call will fill the lower one.

Don't forget to provide also the stressMatrix, needed for stress computation.

``````public class MyHexa8 : Hexa8
{

public MyHexa8(int nodeIndex1, int nodeIndex2, int nodeIndex3, int nodeIndex4, int nodeIndex5, int nodeIndex6, int nodeIndex7, int nodeIndex8, Material mat)
: base(nodeIndex1, nodeIndex2, nodeIndex3, nodeIndex4, nodeIndex5, nodeIndex6, nodeIndex7, nodeIndex8, mat)
{

}

public override void CalcStiffness(Point3D[] nodes, int elemIndex)
{
const int size = 24;

K = new double[size, size];
stressMatrix = new double[8, size, 6];

for (int i = 0; i < size; i++)
{
for (int j = 0; j < size; j++)
{
if (j <= i)
{
K[i, j] = 1;
}
}
}

for (int i0 = 0; i0 < stressMatrix.GetLength(0); i0++)
for (int i1 = 0; i1 < stressMatrix.GetLength(1); i1++)
for (int i2 = 0; i2 < stressMatrix.GetLength(2); i2++)
{
stressMatrix[i0, i1, i2] = 1;
}

StiffnessComputation(nodes);
}
}
``````

The following code shows how to reproduce Eyeshot's element Tria3 defining a custom element MyTria3.

``````public class MyTria3 : Tria3
{
//new public MyMaterial Material;
public MyTria3(int nodeIndex1, int nodeIndex2, int nodeIndex3, Material mat, Double thickness, elementType elType)
: base(nodeIndex1, nodeIndex2, nodeIndex3, mat, thickness, elType)
{

}

public override void CalcStiffness(Point3D[] nodes, int elemIndex)
{
const int size = 6;

K = new double[size, size];
stressMatrix = new double[1, size, 3];

// declare the matrices, using mathnet package for matrix operations
MathNet.Numerics.LinearAlgebra.Matrix Bm;
MathNet.Numerics.LinearAlgebra.Matrix BmT;
MathNet.Numerics.LinearAlgebra.Matrix Kem = MathNet.Numerics.LinearAlgebra.Matrix.Build.Dense(6, 6, 1);
MathNet.Numerics.LinearAlgebra.Matrix Em = MathNet.Numerics.LinearAlgebra.Matrix.Build.Dense(3, 3, 0);

//compute area of the triangle
double[] AreaT = new double[3];

AreaT[0] = nodes[Connection[1]].X * nodes[Connection[2]].Y - nodes[Connection[2]].X * nodes[Connection[1]].Y;
AreaT[1] = nodes[Connection[2]].X * nodes[Connection[0]].Y - nodes[Connection[0]].X * nodes[Connection[2]].Y;
AreaT[2] = nodes[Connection[0]].X * nodes[Connection[1]].Y - nodes[Connection[1]].X * nodes[Connection[0]].Y;

double Area = 0.5 * (AreaT[0] + AreaT[1] + AreaT[2]);

// building custom B-like matrix
Bm = MathNet.Numerics.LinearAlgebra.Matrix.Build.Dense(3, 6, 0);

double doubleArea = 2 * Area;

Bm[0, 0] = (nodes[Connection[1]].Y - nodes[Connection[2]].Y) / doubleArea;
Bm[0, 2] = (nodes[Connection[2]].Y - nodes[Connection[0]].Y) / doubleArea;
Bm[0, 4] = (nodes[Connection[0]].Y - nodes[Connection[1]].Y) / doubleArea;

Bm[1, 1] = (nodes[Connection[2]].X - nodes[Connection[1]].X) / doubleArea;
Bm[1, 3] = (nodes[Connection[0]].X - nodes[Connection[2]].X) / doubleArea;
Bm[1, 5] = (nodes[Connection[1]].X - nodes[Connection[0]].X) / doubleArea;

Bm[2, 0] = (nodes[Connection[2]].X - nodes[Connection[1]].X) / doubleArea;
Bm[2, 1] = (nodes[Connection[1]].Y - nodes[Connection[2]].Y) / doubleArea;
Bm[2, 2] = (nodes[Connection[0]].X - nodes[Connection[2]].X) / doubleArea;
Bm[2, 3] = (nodes[Connection[2]].Y - nodes[Connection[0]].Y) / doubleArea;
Bm[2, 4] = (nodes[Connection[1]].X - nodes[Connection[0]].X) / doubleArea;
Bm[2, 5] = (nodes[Connection[0]].Y - nodes[Connection[1]].Y) / doubleArea;

// get D matrix from the element's material
Em = MathNet.Numerics.LinearAlgebra.Matrix.Build.DenseOfArray(this.Material.Matrix);

// create transposed B-like matrix
BmT = Bm.Transpose();

// compute the product Em*Bm, that will be used also for the stresses
MathNet.Numerics.LinearAlgebra.Matrix EmBm = (Em * Bm);

// calculate stiffness matrix
Kem = Thickness * Area * (BmT *  EmBm);

//export from mathnet matrix to K
K = Kem.ToArray();

for (int i0 = 0; i0 < stressMatrix.GetLength(0); i0++)
for (int i1 = 0; i1 < stressMatrix.GetLength(1); i1++)
for (int i2 = 0; i2 < stressMatrix.GetLength(2); i2++)
{
stressMatrix[i0, i1, i2] += EmBm[i2, i1];
}
StiffnessComputation(nodes);
}
}
``````
0 out of 0 found this helpful
Have more questions? Submit a request