考虑重力荷载的二维/三维有限元分析C++编程实现

//===========================================================================//   Program No.11 General two-(plane strain) or three-dimensional an alysis//             of elastic solids(optional gravity loading).//===========================================================================#include "geom.h"#include "main.h"#include <fstream>#include <iostream>#include <string>#include <utility>#include <vector>#include <iomanip>using namespace std;
int main() {  int fixed_freedoms;//固定位移数  int loaded_nodes;//受载节点数  int ndim;//维数  int ndof;//单元自由度  int nels;//单元数  int neq;//网格自由度  int nip;//单元积分节点数  int nn;                            //节点数  int nod;//单元节点数  int nodof;//节点自由度  int nprops = 3;//材料特性数  int np_types;//材料类型数  int nr;                            //约束节点数  int nst;//应力（应变）项数目  double det;//雅可比矩阵行列式  double penalty = 1.0e20;           //罚因子  string element;//单元类型  string fin_name, fout_name;//读入文件名，写出文件名  pair<string, string> argv;//文件名返回值
//========================= input and initialisation =========================  argv = getname();  fin_name = argv.first;  ifstream fin(fin_name, ios::in);  fout_name = argv.second;  ofstream fout(fout_name, ios::out);
  fin >> element >> nod;  fin >> nels >> nn >> nip >> nodof >> nst >> ndim >> np_types;  ndof = nod * nodof;               //单元自由度  vector<double> eld(ndof);//单元节点位移（单元自由度）  vector<double> eps(nst);//应变项（应力/应变项数目）  vector<int> etype(nels, 1);//单元特性类型向量（单元数）  vector<double> fun(nod);//形函数（单元节点数）  vector<int> g(ndof);//单元定位向量（单元自由度）  vector<double> gc(ndim);//积分点坐标（维数）  vector<int> num(nod);//单元节点编号向量（单元节点数）  vector<double> sigma(nst);//应力项（应力/应变项数）  vector<double> weights(nip);//数值积分加权系数（单元积分节点数）  vector<vector<double> > bee(nst, vector<double>(ndof));//应变-位移矩阵（应力/应变项数, 单元自由度）  vector<vector<double> > coord(nod, vector<double>(ndim));//单元节点坐标（单元节点数, 维数）  vector<vector<double> > dee(nst, vector<double>(nst));//应力-应变矩阵（应力/应变项数, 应力/应变项数）  vector<vector<double> > der(nod, vector<double>(ndim));//局部坐标下形函数的偏导数（单元节点数, 维数）  vector<vector<double> > deriv(nod, vector<double>(ndim));//整体坐标下形函数的偏导数（单元节点数, 维数）  vector<vector<double> > g_coord(nn, vector<double>(ndim));//所有单元的节点坐标（节点数, 维数）  vector<vector<int> > g_g(nels, vector<int>(ndof));//总单元定位矩阵（单元数，单元自由度）  vector<vector<int> > g_num(nels, vector<int>(nod));//总单元节点编号矩阵（单元数, 单元节点数）  vector<vector<double> > jac(ndim, vector<double>(ndim));//雅可比矩阵（维数, 维数）  vector<vector<double> > km(ndof, vector<double>(ndof));//单元刚度矩阵（单元自由度, 单元自由度）  vector<vector<int> > nf(nn, vector<int>(nodof, 1));//节点自由度矩阵（节点数，节点自由度）  vector<vector<double> > points(nip, vector<double>(ndim));//单元积分点坐标（单元的积分点数, 维数）  vector<vector<double> > prop(np_types, vector<double>(nprops));//单元特征矩阵（材料类型数, 材料特性数）
  for (int i = 0; i < np_types; i++) {//从文件中读入:单元特征矩阵    for (int j = 0; j < nprops; j++)      fin >> prop[i][j];  }  if (np_types > 1) {//从文件中读入:单元特性类型向量    for (int i = 0; i < nels; i++)      fin >> etype[i];  }  for (int i = 0; i < nn; i++) {    for (int j = 0; j < ndim; j++) {      fin >> g_coord[i][j];       //从文件中读入:各节点的坐标    }  }  for (int i = 0; i < nels; i++) {    for (int j = 0; j < nod; j++) {      fin >> g_num[i][j];//从文件中读入:各单元的节点编号    }  }  fin >> nr;//从文件中读入:约束节点数  vector<int> k1(nr);//约束节点编号向量（约束节点数）  if (nr > 0) {//从文件中读入:约束节点编号及自由度    for (int i = 0; i < nr; i++) {      fin >> k1[i];      for (int j = 0; j < nodof; j++) {        int m = k1[i] - 1;        fin >> nf[m][j];      }    }  }  auto result0 = formnf(nf);  neq = result0.first;//返回：网格自由度数目  nf = result0.second;//生成：节点自由度矩阵  vector<int> kdiag(neq);//对角项局部向量（网格自由度）  vector<double> loads(neq);//总体荷载/位移向量（网格自由度）  vector<double> gravlo(neq);//总体重力荷载向量（网格自由度）
//=============== loop the elements to find global array sizes ================  for (int i = 0; i < nels; i++) {    num = g_num[i];//第i个单元的节点编号向量    g = num_to_g(num, nf);//生成：单元定位向量    kdiag = fkdiag(g, kdiag);//生成：对角项局部向量    g_g[i] = g;//第i个单元的定位向量  }  for (int i = 1; i < neq; i++)//生成：对角项局部向量    kdiag[i] = kdiag[i] + kdiag[i - 1];  fout << "There are " << neq << " equations and the skyline storage is "       << kdiag[neq - 1] << endl;  cout << "There are " << neq << " equations and the skyline storage is "       << kdiag[neq - 1] << endl;
//================ element stiffness integration and assembly ================  vector<double> kv(kdiag[neq - 1]);//初始化：整体刚度向量  auto result1 = sample(element, nip, ndim);  points = result1.first;//生成：单元积分点坐标  weights = result1.second;//生成：加权系数
  for (int i = 0; i < nels; i++) {    int m = etype[i] - 1;    dee = deemat(prop[m][0], prop[m][1], nst);       //生成：应力-应变矩阵[D]    num = g_num[i];//第i个单元的节点编号向量    g = g_g[i];//第i个单元的定位向量    for (int j = 0; j < nod; j++) {      int k = num[j] - 1;      coord[j] = g_coord[k];//第i个单元的节点坐标矩阵    }
    fill(eld.begin(), eld.end(), 0.0);    fill(km.begin(), km.end(), vector<double>(ndof, 0.0));    vector<vector<double> > btdb(ndof, vector<double>(ndof, 0.0));
    for (int ii = 0; ii < nip; ii++) {      fun = shape_fun(ii, nod, points);//生成：形函数      der = shape_der(ii, nod, points);//生成：局部坐标下形函数的偏导数      jac = matmul(transpose(der), coord);       //生成：雅可比矩阵      det = determinant(jac);//生成：雅可比矩阵行列式      jac = invert(jac);//生成：雅可比矩阵的逆      deriv = matmul(der, jac);//生成：整体坐标下形函数的偏导数      bee = beemat(nst, ndof, deriv);//生成：应变-位移矩阵[B]      btdb = matmul(matmul(transpose(bee), dee), bee);      for (int jj = 0; jj < ndof; jj++) {        for (int kk = 0; kk < ndof; kk++)          btdb[jj][kk] *= det * weights[ii];      }      for (int jj = 0; jj < ndof; jj++) {        for (int kk = 0; kk < ndof; kk++)          km[jj][kk] += btdb[jj][kk];      }//生成：单元刚度矩阵      for (int jj = nodof - 1; jj < ndof; jj += nodof)         eld[jj] += fun[jj / nodof] * det * weights[ii];    }    kv = fsparv(km, g, kdiag, kv);//生成：整体刚度矩阵    for (int k = nodof - 1; k < ndof; k += nodof) {      if (g[k]) {        int n = g[k] - 1;        gravlo[n] -= eld[k] * prop[m][2];//生成：总体网格的重力荷载向量      }    }  }  fin >> loaded_nodes;//从文件读入：受载节点数  vector<int> k2(loaded_nodes);//初始化：受载节点编号向量  if (loaded_nodes) {//从文件读入：受载节点编号，相应荷载向量    for (int i = 0; i < loaded_nodes; i++) {      fin >> k2[i];      int m1 = k2[i] - 1;      for (int j = 0; j < nodof; j++) {        if (nf[m1][j]) {          int m2 = nf[m1][j] - 1;          fin >> loads[m2];        }        else          fin.ignore(5);      }    }  }  for (int k = 0; k < neq; k++)    loads[k] += gravlo[k];//考虑重力荷载  if (ndim == 3)//生成：三维模型的前处理文件Ensight Gold格式    mesh_ensi(fin_name, g_coord, g_num, element, etype, nf, loads, 1, 1, 0.0, true);  else    mesh(g_coord, g_num, fin_name);//生成：二维模型的初始网格图(*_mesh.ps)fin >> fixed_freedoms;//从文件读入：固定位移的数量  if (fixed_freedoms) {    vector<int> node(fixed_freedoms);//固定节点向量    vector<int> no(fixed_freedoms);//固定自由度编号向量    vector<int> sense(fixed_freedoms);//固定自由度的检测向量    vector<double> value(fixed_freedoms);//给定位移向量    for (int i = 0; i < fixed_freedoms; i++) {      fin >> node[i] >> sense[i] >> value[i];      int j = node[i] - 1;      int k = sense[i] - 1;      no[i] = nf[j][k];      int m = no[i] - 1;      int n = kdiag[m] - 1;      kv[n] += penalty;      loads[m] = kv[n] * value[i];    }  }
//============================ equation solution =============================  kv = sparin(kv, kdiag);//整体刚度矩阵kv的Cholesky三角分解形式  loads = spabac(kv, loads, kdiag);//求解有限元方程，得到总体 位移向量  if (ndim == 3) {    fout << endl << "  Node    x-disp        y-disp        z-disp" << endl;    cout << endl << "  Node    x-disp        y-disp        z-disp" << endl;  }  else {    fout << endl << "  Node    x-disp        y-disp" << endl;    cout << endl << "  Node    x-disp        y-disp" << endl;  }  for (int i = 0; i < nn; i++) {    fout << setw(4) << i + 1;    cout << setw(4) << i + 1;    for (int j = 0; j < nodof; j++) {      if (nf[i][j]) {        int k = nf[i][j] - 1;        fout << setw(14) << scientific << setprecision(4) << loads[k];        cout << setw(14) << scientific << setprecision(4) << loads[k];      }      else {        fout << "    0.0000e+00";        cout << "    0.0000e+00";      }    }    fout << endl;    cout << endl;  }
//================ recover stresses at element Gauss - points ================//nip = 1;//points.resize(nip, vector<double>(ndim));//weights.resize(nip);//auto result2 = sample(element, nip, ndim);//points = result2.first;//生成：单元积分点坐标//weights = result2.second;//生成：加权系数  fout << endl << " The integration point ( nip = " << nip << " ) stresses are:" << endl;  cout << endl << " The integration point ( nip = " << nip << " ) stresses are:" << endl;  if (ndim == 3) {    fout << " Element  x-coord       sig_x         tau_xy" << endl;    cout << " Element  x-coord       sig_x         tau_xy" << endl;    fout << "          y-coord       sig_y         tau_yz" << endl;    cout << "          y-coord       sig_y         tau_yz" << endl;    fout << "          z-coord       sig_z         tau_zx" << endl;    cout << "          z-coord       sig_z         tau_zx" << endl;  }  else {    fout << " Element  x-coord       y-coord        sig_x         sig_y        tau_xy" << endl;    cout << " Element  x-coord       y-coord        sig_x         sig_y        tau_xy" << endl;  }
  for (int i = 0; i < nels; i++) {    if (np_types > 1) {      int m = etype[i] - 1;      dee = deemat(prop[m][0], prop[m][1], nst);    }//生成：应力-应变矩阵[D]    else      dee = deemat(prop[0][0], prop[0][1], nst);    num = g_num[i];//生成：第i个单元的节点编号向量    for (int j = 0; j < nod; j++) {      int k = num[j] - 1;      coord[j] = g_coord[k];//生成：第i个单元的节点坐标矩阵    }    g = g_g[i];//生成：第i个单元的定位向量    for (int k = 0; k < ndof; k++) {      if (g[k]) {        int m = g[k] - 1;        eld[k] = loads[m];//生成：第i个单元的位移向量      }      else        eld[k] = 0.0;    }    for (int ii = 0; ii < nip; ii++) {      fout << setw(4) << i + 1;      cout << setw(4) << i + 1;      fun = shape_fun(ii, nod, points);//生成：形函数      der = shape_der(ii, nod, points);//生成：局部坐标下形函数的偏导数      gc = matmul(fun, coord);//生成：积分点坐标      jac = matmul(transpose(der), coord);       //生成：雅可比矩阵      jac = invert(jac);//生成：雅可比矩阵的逆      deriv = matmul(der, jac);//生成：整体坐标下形函数的偏导数      bee = beemat(nst, ndof, deriv);//生成：应变-位移矩阵[B]      eps = matmul(bee, eld);//生成：单元应变向量      sigma = matmul(dee, eps);//生成：单元应力向量      if (ndim == 3) {        fout << setw(14) << scientific << setprecision(4) << gc[0]             << setw(14) << scientific << setprecision(4) << sigma[0]             << setw(14) << scientific << setprecision(4) << sigma[3] << endl;        cout << setw(14) << scientific << setprecision(4) << gc[0]             << setw(14) << scientific << setprecision(4) << sigma[0]             << setw(14) << scientific << setprecision(4) << sigma[3] << endl;        fout << "    " << setw(14) << scientific << setprecision(4) << gc[1]             << setw(14) << scientific << setprecision(4) << sigma[1]             << setw(14) << scientific << setprecision(4) << sigma[4] << endl;        cout << "    " << setw(14) << scientific << setprecision(4) << gc[1]             << setw(14) << scientific << setprecision(4) << sigma[1]             << setw(14) << scientific << setprecision(4) << sigma[4] << endl;        fout << "    " << setw(14) << scientific << setprecision(4) << gc[2]             << setw(14) << scientific << setprecision(4) << sigma[2]             << setw(14) << scientific << setprecision(4) << sigma[5] << endl;        cout << "    " << setw(14) << scientific << setprecision(4) << gc[2]             << setw(14) << scientific << setprecision(4) << sigma[2]             << setw(14) << scientific << setprecision(4) << sigma[5] << endl;      }      else {        for (int n = 0; n < ndim; n++) {          fout << setw(14) << scientific << setprecision(4) << gc[n];          cout << setw(14) << scientific << setprecision(4) << gc[n];        }        for (int n = 0; n < nst; n++) {          fout << setw(14) << scientific << setprecision(4) << sigma[n];          cout << setw(14) << scientific << setprecision(4) << sigma[n];        }        fout << endl;        cout << endl;      }    }  }
  if (ndim == 3)    dis msh_ensi(fin_name, 1, nf, loads);//生成：三维模型的后处理文件Ensight Gold格式  else {    dis msh(loads, nf, 0.05, g_coord, g_num, fin_name);//生成：二维模型的变形网格图(*_dis.ps)    vecmsh(loads, nf, 0.05, 0.1, g_coord, g_num, fin_name);//生成：二维模型的位移矢量图(*_vec.ps)  }  fin.close();  fout.close();}