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);
}
}
Comments
Please sign in to leave a comment.