bandedstiffnessmatrix.cpp

有限元计算程序用于计算平面板结构单元的应力和位移。

#include <iostream>
using namespace std;
#include "Plate.h"

//---------------------------------------------------------------------------

double** TPlate::Banded_Stiffness_Matrix(){

        S = Allocate_2D_Matrix(S, (NEQ + 1), (MBAND + 1));
        S = SetZero_2D_Matrix(S, (NEQ + 1), (MBAND + 1));

        double DiffX, DiffY;
        double a, b;
        double Stiff = 0;
        double k11, k12, k13, k14, k15, k17, k22, k24, k28;
        int K, L;
        int *KK;
        double **SE;

        KK = Allocate_1D_Matrix(KK, (Quad * DOF));
        SE = Allocate_2D_Matrix(SE, (Quad * DOF), (Quad * DOF));

        KK = SetZero_1D_Matrix(KK, (Quad * DOF));
        SE = SetZero_2D_Matrix(SE, (Quad * DOF), (Quad * DOF));

        Stiff = E * t; 

        for (int i = 0; i < Elements; i++){

            for (int j = 0; j < Quad; j++){
                d[j] = Connect[i][j] - 1;
            }

            DiffX = XYZ[d[1]][0] - XYZ[d[0]][0];
            DiffY = XYZ[d[2]][1] - XYZ[d[1]][1];

            a = 0.5 * DiffX;    aa = a;
            b = 0.5 * DiffY;    bb = b;

            k11 = (((a*a)*(1-v)) + (2*(b*b))) / (6*a*b*(1-(v*v)));
            k12 = 1 / (8*(1-v));
            k13 = (((a*a)*(1-v)) - (4*(b*b))) / (12*a*b*(1-(v*v)));
            k14 = (3*v - 1) / (8*(1-(v*v)));
            k15 = -1 * (((a*a)*(1-v)) + (2*(b*b))) / (12*a*b*(1-(v*v)));
            k17 = ((b*b) - (a*a)*(1-v)) / (6*a*b*(1-(v*v)));
            k22 = ((2*(a*a)) + (b*b)*(1-v)) / (6*a*b*(1-(v*v)));
            k24 = ((a*a) - (b*b)*(1-v)) / (6*a*b*(1-(v*v)));
            k28 = (((b*b)*(1-v)) - (4*(a*a))) / (12*a*b*(1-(v*v)));

            SE[0][0] = Stiff * k11;
            SE[0][1] = Stiff * k12;
            SE[0][2] = Stiff * k13;
            SE[0][3] = Stiff * k14;
            SE[0][4] = Stiff * k15;
            SE[0][5] = -1 * Stiff * k12;
            SE[0][6] = Stiff * k17;
            SE[0][7] = -1 * Stiff * k14;

            SE[1][1] = Stiff * k22;
            SE[1][2] = -1 * Stiff * k14;
            SE[1][3] = Stiff * k24;
            SE[1][4] = -1 * Stiff * k12;
            SE[1][5] = (-0.5) * Stiff * k22;
            SE[1][6] = Stiff * k14;
            SE[1][7] = Stiff * k28;

            SE[2][2] = Stiff * k11;
            SE[2][3] = -1 * Stiff * k12;
            SE[2][4] = Stiff * k17;
            SE[2][5] = Stiff * k14;
            SE[2][6] = Stiff * k15;
            SE[2][7] = Stiff * k12;

            SE[3][3] = Stiff * k22;
            SE[3][4] = -1 * Stiff * k14;
            SE[3][5] = Stiff * k28;
            SE[3][6] = Stiff * k12;
            SE[3][7] = (-0.5) * Stiff * k22;

            SE[4][4] = Stiff * k11;
            SE[4][5] = Stiff * k12;
            SE[4][6] = Stiff * k13;
            SE[4][7] = Stiff * k14;

            SE[5][5] = Stiff * k22;
            SE[5][6] = -1 * Stiff * k14;
            SE[5][7] = Stiff * k24;

            SE[6][6] = Stiff * k11;
            SE[6][7] = -1 * Stiff * k12;

            SE[7][7] = Stiff * k22;

            for (int i = 7; i >= 1; i--){
                for (int j = 6; j >= 0; j--){
                    SE[i][j] = SE[j][i];
                }
            }

            // Assembly active stiffness equation in banded format
            // (pp. 51  Figure 2.10-5)

            KK[0] = ID[0][d[0]];
            KK[1] = ID[1][d[0]];
            KK[2] = ID[0][d[1]];
            KK[3] = ID[1][d[1]];
            KK[4] = ID[0][d[2]];
            KK[5] = ID[1][d[2]];
            KK[6] = ID[0][d[3]];
            KK[7] = ID[1][d[3]];

            for (int p = 0; p < (Quad * DOF); p++){
                if (KK[p] <= 0) goto esc3;
                K = KK[p];
                for (int q = 0; q < (Quad * DOF); q++){
                    if (KK[q] < K) goto esc4;
                    L = KK[q] - K + 1;
                    S[K][L] = S[K][L] + SE[p][q];
esc4:           }
esc3:       }
        }

        De_Allocate_1D_Matrix(KK);
        De_Allocate_2D_Matrix(SE, (Quad * DOF));

        return S;
}

//---------------------------------------------------------------------------