gesolver.cpp

来自「有限元计算程序用于计算平面板结构单元的应力和位移。」· C++ 代码 · 共 77 行

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

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

double* TPlate::Cal_GESolver(){

        
        int NEQP = 0, NEQM = 0;
        int LIM = 0;
        int I, J, K;
        double DUM;

        if (MBAND > 1) goto esc5;
           for (int N = 1; N <= NEQ; N++){
               U[N] = U[N] / S[N][1];
           }
           return U;

esc5:
        NEQP = NEQ + 1;
        NEQM = NEQ - 1;

        // Forward Reduction of the Coefficient Matrix [S];
        for (int N = 1; N <= NEQM; N++){
            LIM = MIN(MBAND, (NEQP-N));
            for (int L = 2; L <= LIM; L++){
                DUM = S[N][L] / S[N][1];
                I = N + L - 1;
                J = 0;
                for (int K = L; K <= LIM; K++){
                    J = J + 1;
                    S[I][J] = S[I][J] - (DUM * S[N][K]);
                }
                S[N][L] = DUM;
            }
        }

        // Forward Reduction of the Constant Vector {U};
        for (int N = 1; N <= NEQM; N++){
            LIM = MIN(MBAND, (NEQP-N));
            for (int L = 2; L <= LIM; L++){
                I = N + L - 1;
                U[I] = U[I] - (S[N][L] * U[N]);
            }
            U[N] = U[N] / S[N][1];
        }

        U[NEQ] = U[NEQ] / S[NEQ][1];

        // Backward Substitution. Formwe Unknowns {X} Overwrite {C};
        for (int N = NEQM; N >= 1; N--){
            LIM = MIN(MBAND, (NEQP-N));
            for (int L = 2; L <= LIM; L++){
                K = N + L - 1;
                U[N] = U[N] - (S[N][L] * U[K]);
            }
        }

        return U;
}

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

int TPlate::MIN(int x, int y){

        if (x < y){
           return x;
        } else {
           return y;
        }
}

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