linbcg.cpp

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

#include <iostream>
using namespace std;
#define EPS 1.0E-14
#include "Plate.h"

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

void TPlate::linbcg(){

        // Arrays Variables Declaration;
        double *r, *rr;
        double *z, *zz;
        double *p, *pp;
        int iter;
        int itol;
        double bnrm;
        double znrm;
        int itmax = 100;
        double bknum;
        double bk;
        double bkden = 1.0;
        double akden;
        double ak;
        double err;
        double zm1nrm;
        double dxnrm;
        double xnrm;
        double tol = 1.0E-10;

        // Create 1D Arrays;
        r = Allocate_1D_Matrix(r, NEQ);
        rr = Allocate_1D_Matrix(rr, NEQ);
        z = Allocate_1D_Matrix(z, NEQ);
        zz = Allocate_1D_Matrix(zz, NEQ);
        p = Allocate_1D_Matrix(p, NEQ);
        pp = Allocate_1D_Matrix(pp, NEQ);

        // Set 1D Arrays to Zero;
        r = SetZero_1D_Matrix(r, NEQ);
        rr = SetZero_1D_Matrix(rr, NEQ);
        z = SetZero_1D_Matrix(z, NEQ);
        zz = SetZero_1D_Matrix(zz, NEQ);
        p = SetZero_1D_Matrix(p, NEQ);
        pp = SetZero_1D_Matrix(pp, NEQ);

        cout << "itol (select 1 ~ 4): ";
        cin >> itol;
        if (itol < 1 || itol > 4)
        exit(1);

        iter = 0;

        atimes (x, r, 0);

        for (int j = 0; j < NEQ; j++){
            cout << "r[" << j << "] = " << r[j] << endl;
        }

        for (int j = 0; j < NEQ; j++){
            r[j] = b[j] - r[j];
            rr[j] = r[j];
        }

        if (itol == 1){
           bnrm = snrm(b, itol);
           asolve(r, z, 0);
        }
        else if (itol == 2){
           asolve(b, z, 0);
           bnrm = snrm(z, itol);
           asolve(r, z, 0);
        }
        else if (itol == 3 || itol == 4){
           asolve(b, z, 0);
           bnrm = snrm(z, itol);
           asolve(r, z, 0);
           znrm = snrm(z, itol);
        }
        
        while (iter < itmax){
              ++iter;

              asolve(rr, zz, 1);
              bknum = 0.0;
              
              for (int j = 0; j < NEQ; j++){
                bknum += z[j] * rr[j];
              }

              if (iter == 1){
                for (int j = 0; j < NEQ; j++){
                    p[j] = z[j];
                    pp[j] = zz[j];
                }
              } else {
                bk = bknum / bkden;
                for (int j = 0; j < NEQ; j++){
                    p[j] = bk * p[j] + z[j];
                    pp[j] = bk * pp[j] + zz[j];
                }
              }

              bkden = bknum;
              atimes(p, z, 0);
              akden = 0.0;

              for (int j = 0; j < NEQ; j++){
                  akden += z[j] * pp[j];
              }

              ak = bknum / akden;                       
              atimes(pp, zz, 1);

              for (int j = 0; j < NEQ; j++){
                  x[j] += ak * p[j];   
                  r[j] -= ak * z[j];
                  rr[j] -= ak * zz[j];
              }

              asolve(r, z, 0);

              if (itol == 1){
                 err = snrm(r, itol) / bnrm;
              }
              else if (itol == 2){
                 err = snrm(z, itol) / bnrm;
              }
              else if (itol == 3 || itol == 4){

                 zm1nrm = znrm;
                 znrm = snrm(z, itol);

                 if (fabs(zm1nrm - znrm) > EPS * znrm){
                    dxnrm = fabs(ak) * snrm(p, itol);
                    err = znrm / fabs(zm1nrm - znrm) * dxnrm;
                 } else {
                    err = znrm / bnrm;
                    continue;
                 }

                 xnrm = snrm(x, itol);

                 if (err <= 0.5 * xnrm){
                    err /= xnrm;
                 } else {
                    err = znrm / bnrm;
                    continue;
                 }
              }
              if (err <= tol) break;
        }

        // Deallocate Arrays;
        De_Allocate_1D_Matrix(r);
        De_Allocate_1D_Matrix(rr);
        De_Allocate_1D_Matrix(z);
        De_Allocate_1D_Matrix(zz);
        De_Allocate_1D_Matrix(p);
        De_Allocate_1D_Matrix(pp);
}

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

void TPlate::atimes(double *x, double *r, const int itrnsp){

        if (itrnsp == 0){
           dsprsax(sa, ija, x, r);
        } else {
           dsprstx(sa, ija, x, r);
        }
}

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

void TPlate::dsprsax(double *sa, int *ija, double *x, double *b){

     for (int i = 0; i < NEQ; i++){
         b[i] = sa[i] * x[i];
         for (int k = ija[i]; k < ija[i+1]; k++){
             b[i] += sa[k] * x[ija[k]];
         }
     }
}

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

void TPlate::dsprstx(double *sa, int *ija, double *x, double *b){

        int j = 0;
        for (int i = 0; i < NEQ; i++){
            b[i] = sa[i] * x[i];
        }
        for (int i = 0; i < NEQ; i++){
            for (int k = ija[i]; k < ija[i+1]; k++){
                j = ija[k];
                b[j] += sa[k] * x[i];
            }
        }
}

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

double TPlate::snrm(double *sx, const int itol){

        int isamax;
        double ans;
        if (itol <= 3){
           ans = 0.0;
           for (int i = 0; i < NEQ; i++){
               ans += sx[i] * sx[i];
           }
           return sqrt(ans);
        } else {
           isamax = 0;
           for (int i = 0; i < NEQ; i++){
               if (fabs(sx[i]) > fabs(sx[isamax])){
                  isamax = i;
               }
           }
           return fabs(sx[isamax]);
        }
}

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

void TPlate::asolve(double *b, double *x, const int itrnsp){

        for (int i = 0; i < NEQ; i++){
            x[i] = (sa[i] != 0.0 ? b[i]/ sa[i]:b[i]);
        }
}

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