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); 
        } 
    }