using MathNet.Numerics.LinearAlgebra;

using System;

​

public class SAP2000PlateElement

{

    // NodeCoordinates: assumed as a 4x2 array [node, (x,y)]

    public double[,] NodeCoordinates { get; }

    public double Thickness { get; }

    public double YoungsModulus { get; }

    public double PoissonsRatio { get; }

​

    // A small drilling stiffness penalty (this value may be tuned)

    public double DrillingStiffnessFactor { get; } = 1e-4;

​

    public SAP2000PlateElement(double[,] nodeCoordinates, double thickness, double youngsModulus, double poissonsRatio)

    {

        NodeCoordinates = nodeCoordinates;

        Thickness = thickness;

        YoungsModulus = youngsModulus;

        PoissonsRatio = poissonsRatio;

    }

​

    /// <summary>

    /// Computes the combined element stiffness matrix that includes membrane, bending, shear,

    /// and a small drilling stiffness for stabilization.

    /// </summary>

    public Matrix<double> ComputeStiffnessMatrix()

    {

        // Six DOF per node (u, v, w, θx, θy, θz)

        var dofPerNode = 6;

        var totalDof = 4 * dofPerNode;

        var ke = Matrix<double>.Build.Dense(totalDof, totalDof);

​

        // Gauss integration (2x2)

        double[][] gaussPoints =

        [

            [-1.0 / Math.Sqrt(3), -1.0 / Math.Sqrt(3)],

            [1.0 / Math.Sqrt(3), -1.0 / Math.Sqrt(3)],

            [1.0 / Math.Sqrt(3),  1.0 / Math.Sqrt(3)],

            [-1.0 / Math.Sqrt(3),  1.0 / Math.Sqrt(3)]

        ];

        var gaussWeight = 1.0;

​

        foreach (var point in gaussPoints)

        {

            var xi = point[0];

            var eta = point[1];

​

            // --- Compute shape function derivatives and shape functions ---

            var dN_nat = ComputeShapeFunctionDerivatives(xi, eta);  // 2x4 matrix (d/dxi, d/deta)

            var N = ComputeShapeFunctions(xi, eta);             // length 4

​

            // Jacobian: transformation from (ξ,η) to (x,y)

            var J = ComputeJacobian(dN_nat);

            var detJ = J.Determinant();

            var invJ = J.Inverse();

​

            // Global derivatives (dN/dx, dN/dy)

            var dN_global = invJ * dN_nat;  // 2x4

​

            // Differential area in global coordinates.

            var dA = detJ * gaussWeight;

​

            // --- Assemble individual stiffness contributions ---

            // 1. Membrane contribution (in-plane: u, v → DOFs 0 and 1)

            var B_mem = ComputeMembraneBMatrix(dN_global);

            var D_mem = ComputeMembraneMaterialMatrix();

            ke += AssembleFullMatrix(B_mem.Transpose() * D_mem * B_mem * dA, 0);

​

            // 2. Bending contribution (rotations: θx, θy → DOFs 3 and 4)

            var B_bend = ComputeBendingBMatrix(dN_global);

            var D_bend = ComputeBendingMaterialMatrix();

            ke += AssembleFullMatrix(B_bend.Transpose() * D_bend * B_bend * dA, 3);

​

            // 3. Shear contribution (transverse displacement w (DOF 2) and rotations θx, θy (DOF 3 & 4))

            var B_shear = ComputeShearBMatrix(dN_global, N);

            var D_shear = ComputeShearMaterialMatrix();

            ke += AssembleFullMatrix(B_shear.Transpose() * D_shear * B_shear * dA, 2);

​

            // Note: The function AssembleFullMatrix embeds the element contribution

            // into the global 4-node, 6-dof-per-node stiffness. For each contribution, only

            // the relevant DOFs in each node are affected.

        }

​

        // 4. Drilling stiffness contribution is added as a penalty term on DOF 5 for each node.

        for (var i = 0; i < 4; i++)

        {

            var idx = i * dofPerNode + 5;

            ke[idx, idx] += DrillingStiffnessFactor * YoungsModulus * Thickness;

        }

​

        return ke;

    }

​

    #region Shape Functions and Jacobian

    /// <summary>

    /// Bilinear shape functions for a quadrilateral element.

    /// </summary>

    private double[] ComputeShapeFunctions(double xi, double eta)

    {

        var N = new double[4];

        N[0] = 0.25 * (1 - xi) * (1 - eta);

        N[1] = 0.25 * (1 + xi) * (1 - eta);

        N[2] = 0.25 * (1 + xi) * (1 + eta);

        N[3] = 0.25 * (1 - xi) * (1 + eta);

        return N;

    }

​

    /// <summary>

    /// Derivatives of the bilinear shape functions with respect to natural coordinates (ξ,η).

    /// Returns a 2x4 matrix, where first row is dN/dxi and second row is dN/deta.

    /// </summary>

    private Matrix<double> ComputeShapeFunctionDerivatives(double xi, double eta)

    {

        var dN = Matrix<double>.Build.Dense(2, 4);

        // dN/dxi

        dN[0, 0] = -0.25 * (1 - eta);

        dN[0, 1] = 0.25 * (1 - eta);

        dN[0, 2] = 0.25 * (1 + eta);

        dN[0, 3] = -0.25 * (1 + eta);

        // dN/deta

        dN[1, 0] = -0.25 * (1 - xi);

        dN[1, 1] = -0.25 * (1 + xi);

        dN[1, 2] = 0.25 * (1 + xi);

        dN[1, 3] = 0.25 * (1 - xi);

        return dN;

    }

​

    /// <summary>

    /// Computes the 2x2 Jacobian matrix using the natural derivatives and nodal coordinates.

    /// </summary>

    private Matrix<double> ComputeJacobian(Matrix<double> dN_nat)

    {

        var J = Matrix<double>.Build.Dense(2, 2);

        for (var i = 0; i < 4; i++)

        {

            var x = NodeCoordinates[i, 0];

            var y = NodeCoordinates[i, 1];

            J[0, 0] += dN_nat[0, i] * x;

            J[0, 1] += dN_nat[0, i] * y;

            J[1, 0] += dN_nat[1, i] * x;

            J[1, 1] += dN_nat[1, i] * y;

        }

        return J;

    }

    #endregion

​

    #region B-Matrices and Material Matrices

    /// <summary>

    /// Membrane B-matrix. This formulation corresponds to a standard plane stress element.

    /// The returned matrix has 3 rows (ε_x, ε_y, γ_xy) and 8 columns (two DOFs per node).

    /// This B-matrix is embedded into the full 4-node (24 DOF) structure in AssembleFullMatrix.

    /// </summary>

    private Matrix<double> ComputeMembraneBMatrix(Matrix<double> dN_global)

    {

        // There are 4 nodes, each contributes u and v.

        var Bm = Matrix<double>.Build.Dense(3, 4 * 2);

        for (var i = 0; i < 4; i++)

        {

            var dN_dx = dN_global[0, i];

            var dN_dy = dN_global[1, i];

            var col = i * 2;

            Bm[0, col] = dN_dx;      // ε_x = du/dx

            Bm[1, col + 1] = dN_dy;      // ε_y = dv/dy

            Bm[2, col] = dN_dy;      // γ_xy = du/dy

            Bm[2, col + 1] = dN_dx;      // γ_xy = dv/dx

        }

        return Bm;

    }

​

    /// <summary>

    /// Returns the plane stress material matrix for the membrane portion.

    /// </summary>

    private Matrix<double> ComputeMembraneMaterialMatrix()

    {

        var coeff = YoungsModulus / (1 - Math.Pow(PoissonsRatio, 2));

        var Dm = Matrix<double>.Build.Dense(3, 3);

        Dm[0, 0] = coeff; Dm[0, 1] = coeff * PoissonsRatio;

        Dm[1, 0] = coeff * PoissonsRatio; Dm[1, 1] = coeff;

        Dm[2, 2] = coeff * (1 - PoissonsRatio) / 2;

        // Note: This matrix is to be multiplied by Thickness.

        return Dm;

    }

​

    /// <summary>

    /// Bending B-matrix. It relates nodal rotations to plate curvatures:

    ///   κₓ = dθₓ/dx,  κ_y = dθ_y/dy,  κ_xy = dθₓ/dy + dθ_y/dx.

    /// This matrix has 3 rows and 8 columns (2 rotational DOFs per node) and is later embedded.

    /// </summary>

    private Matrix<double> ComputeBendingBMatrix(Matrix<double> dN_global)

    {

        var Bb = Matrix<double>.Build.Dense(3, 4 * 2);

        for (var i = 0; i < 4; i++)

        {

            var dN_dx = dN_global[0, i];

            var dN_dy = dN_global[1, i];

            var col = i * 2;

            // κₓ = dθₓ/dx

            Bb[0, col] = dN_dx;

            // κ_y = dθ_y/dy

            Bb[1, col + 1] = dN_dy;

            // κ_xy = dθₓ/dy + dθ_y/dx

            Bb[2, col] = dN_dy;

            Bb[2, col + 1] = dN_dx;

        }

        return Bb;

    }

​

    /// <summary>

    /// Returns the bending material matrix.

    /// </summary>

    private Matrix<double> ComputeBendingMaterialMatrix()

    {

        var coeff = (YoungsModulus * Math.Pow(Thickness, 3)) / (12 * (1 - Math.Pow(PoissonsRatio, 2)));

        var Db = Matrix<double>.Build.Dense(3, 3);

        Db[0, 0] = coeff;

        Db[0, 1] = PoissonsRatio * coeff;

        Db[1, 0] = PoissonsRatio * coeff;

        Db[1, 1] = coeff;

        Db[2, 2] = (1 - PoissonsRatio) * coeff / 2;

        return Db;

    }

​

    /// <summary>

    /// Shear B-matrix for Mindlin-Reissner theory.

    /// Here, the shear strains γ_xz and γ_yz are given by:

    ///   γ_xz = dw/dx - θₓ,   γ_yz = dw/dy - θ_y.

    /// For each node, the contributions come from the derivative of w (DOF 2) and the shape function

    /// (to account for the rotational DOFs at indices 3 and 4).

    /// This matrix has 2 rows and 3 columns per node (w, θₓ, θ_y).

    /// </summary>

    private Matrix<double> ComputeShearBMatrix(Matrix<double> dN_global, double[] N)

    {

        var Bs = Matrix<double>.Build.Dense(2, 4 * 3);

        for (var i = 0; i < 4; i++)

        {

            var dN_dx = dN_global[0, i];

            var dN_dy = dN_global[1, i];

            var col = i * 3;

            // γ_xz = dw/dx - θₓ

            Bs[0, col] = dN_dx;    // from deflection w (to be associated with DOF index 2)

            Bs[0, col + 1] = -N[i];    // contribution for rotation θₓ (to be associated with DOF index 3)

            // γ_yz = dw/dy - θ_y

            Bs[1, col] = dN_dy;    // from deflection w (DOF index 2)

            Bs[1, col + 2] = -N[i];    // contribution for rotation θ_y (DOF index 4)

        }

        return Bs;

    }

​

    /// <summary>

    /// Returns the shear material matrix.

    /// For Mindlin-Reissner, the transverse shear modulus is given by G = E/(2(1+ν)).

    /// A shear correction factor (typically 5/6) is applied.

    /// </summary>

    private Matrix<double> ComputeShearMaterialMatrix()

    {

        var k = 5.0 / 6.0;

        var G = YoungsModulus / (2 * (1 + PoissonsRatio));

        var coeff = k * (YoungsModulus * Thickness) / (2 * (1 + PoissonsRatio));

        // 2x2 identity matrix scaled by coeff

        var Ds = Matrix<double>.Build.DenseIdentity(2) * coeff;

        return Ds;

    }

    #endregion

​

    #region Assembly of Contributions into the Full DOF Ordering

    /// <summary>

    /// Assembles a contribution matrix (which is built for a subset of DOFs) into the full element stiffness matrix.

    /// The submatrix is assumed to be built for a single field per node (e.g. 2 DOFs for membrane, bending, or 3 DOFs for shear),

    /// and this method embeds it into the full ordering for the 4-node element (6 DOF per node).

    /// 

    /// The parameter offset specifies which DOF(s) the submatrix corresponds to:

    /// - For membrane: offset = 0 corresponds to DOFs 0-1 per node.

    /// - For bending: offset = 3 corresponds to DOFs 3-4 per node.

    /// - For shear: offset = 2 corresponds to DOF 2 plus the corresponding rotation DOFs; note that for shear, a mapping is defined below.

    /// </summary>

    private Matrix<double> AssembleFullMatrix(Matrix<double> keSub, int offset)

    {

        // Create a full contribution matrix with dimensions (24 x 24) for 4 nodes x 6 DOF.

        var dofPerNode = 6;

        var totalDof = 4 * dofPerNode;

        var Kfull = Matrix<double>.Build.Dense(totalDof, totalDof);

​

        // Depending on the offset, we know how many DOFs per node are being affected:

        // For membrane and bending, the submatrix has 2 DOFs per node.

        // For shear, the submatrix has 3 DOFs per node.

        int subDofPerNode;

        if (offset == 0 || offset == 3) subDofPerNode = 2;

        else if (offset == 2) subDofPerNode = 3;

        else throw new Exception("Invalid offset in assembly.");

​

        for (var node = 0; node < 4; node++)

        {

            for (var i = 0; i < subDofPerNode; i++)

            {

                for (var j = 0; j < subDofPerNode; j++)

                {

                    var row = node * dofPerNode + (offset + i);

                    var col = node * dofPerNode + (offset + j);

                    // Add the diagonal block from the submatrix.

                    // Here we assume that keSub is organized with each node’s sub-DOFs in sequence.

                    var subIndex = node * subDofPerNode + i;

                    var subCol = node * subDofPerNode + j;

                    Kfull[row, col] += keSub[subIndex, subCol];

                }

            }

        }

        return Kfull;

    }

    #endregion

}