readme

需要添加eigen库，详细见eigen库的使用
两个程序：一个是求解IGA的类（.h文件）IGA_Plate.h ；另一个是主程序
将eigen库添加后，直接运行主程序可出结果

初稿没有达到通用性要求
需要继续修改

问题描述

编程遇到的问题

结果

程序

IGA_Plate.h 文件

#pragma once#include<iostream>
#include <Eigen/Dense>
using namespace Eigen;
using namespace std;
#include <vector>class  IGA_Plate
{public://初始物理量static double E;          //弹性模量static double mu;         //泊松比static double t;          //厚度static MatrixXd Node;     //节点信息矩阵，第一列为节点编号，2-4列分别为x，y（，z）坐标static MatrixXd ele;      //单元信息，第一列为单元编号，后面各列为单元上的节点编号static vector<double> U;   //节点向量static vector<double> V;static int p;              //基函数次数public://int num_GK = ele.rows() * Node.rows();      //总刚的大小  2*12//int num_Ke = ele.rows();                       //单刚的个数  2//int num_ele_node = ele.cols() - 1;             //单刚中节点的个数  9//int size_Ke = ele.rows() * (ele.cols() - 1);   //单刚矩阵的大小MatrixXd D_Matrix(){MatrixXd D(3, 3);D << 1, mu, 0,mu, 1, 0,0, 0, (1 - mu) / 2;D *= (E / (1 - mu * mu));return D;}MatrixXd D = D_Matrix();//求基函数的值 （函数求值）double OneBasisFun_withreturn(int p, int m, vector<double>& U, int i, double u){// p: p次B样条    m：节点向量的个数-1   // U: 节点向量[u0,u1,...,um]   i：区间段[ui ui+1,...,ui+p+1]    // u：求p次B样条基函数 Ni,p 在u点的函数值double Nip;vector<double> N(p + 1);if ((i == 0 && u == U[0]) || (i == m - p - 1 && u == U[m])){Nip = 1.0;return Nip;}if (u < U[i] || u >= U[i + p + 1]){Nip = 0.0;return Nip;}int jj, j;for (jj = 0; jj <= p; jj++) {if (u >= U[i + jj] && u < U[i + jj + 1])N[jj] = 1.0;elseN[jj] = 0.0;}double saved;int k;for (k = 1; k <= p; k++){if (N[0] == 0.0) saved = 0.0;else saved = ((u - U[i]) * N[0]) / (U[i + k] - U[i]);for (j = 0; j < p - k + 1; j++){double Uleft = U[i + j + 1];double Uright = U[i + j + k + 1];if (N[j + 1] == 0.0){N[j] = saved; saved = 0.0;}else{double temp = N[j + 1] / (Uright - Uleft);N[j] = saved + (Uright - u) * temp;saved = (u - Uleft) * temp;}}}Nip = N[0];return Nip;}//求单个基函数的导数double DersOneBasisFun_double_First_Ders(int p, vector<double>& U, int i, double u){// 原先是DersOneBasisFun_double_First_Ders(int p,int m,vector<double> &U,int i,double u,int n)// 但是我看参数 n 没有用到，就将 int n   删除了// m：节点向量的个数-1  没有用到，就将 int m   删除了// p: p次B样条     // U: 节点向量[u0,u1,...,um]   i：区间段[ui ui+1,...,ui+p+1]    // u：求p次B样条基函数 Ni,p 在u点的导数值double DersOnebasisfun_value;if (p == 3) {double Onebasisfun_i_pminus1 = OneBasisFun_withreturn(p - 1, U.size() - 1, U, i, u);double Onebasisfun_iplus1_pminus1 = OneBasisFun_withreturn(p - 1, U.size() - 1, U, i + 1, u);double Uipi = 1.0 / (U[i + p] - U[i]);double Uip1i1 = 1.0 / (U[i + p + 1] - U[i + 1]);if ((U[i + p] - U[i]) == 0)Uipi = 0.0;if ((U[i + p + 1] - U[i + 1]) == 0)Uip1i1 = 0.0;DersOnebasisfun_value = p * (Onebasisfun_i_pminus1 * Uipi - Onebasisfun_iplus1_pminus1 * Uip1i1);if (u == 1.0){u = 0.99999;double Onebasisfun_i_pminus1 = OneBasisFun_withreturn(p - 1, U.size() - 1, U, i, u);double Onebasisfun_iplus1_pminus1 = OneBasisFun_withreturn(p - 1, U.size() - 1, U, i + 1, u);double Uipi = 1.0 / (U[i + p] - U[i]);double Uip1i1 = 1.0 / (U[i + p + 1] - U[i + 1]);if ((U[i + p] - U[i]) == 0)Uipi = 0.0;if ((U[i + p + 1] - U[i + 1]) == 0)Uip1i1 = 0.0;double DersOnebasisfun_value_zero = p * (Onebasisfun_i_pminus1 * Uipi - Onebasisfun_iplus1_pminus1 * Uip1i1);DersOnebasisfun_value = -DersOnebasisfun_value_zero;//cout<< "Hehehehhe" <<  " " << i << " " << Onebasisfun_i_pminus1  << " " << Uipi  << " " <<  Onebasisfun_iplus1_pminus1  << " " << Uip1i1 << " " << DersOnebasisfun_value_zero << "\n";}}if (p == 2) {double Onebasisfun_i_pminus1 = OneBasisFun_withreturn(p - 1, U.size() - 1, U, i, u);double Onebasisfun_iplus1_pminus1 = OneBasisFun_withreturn(p - 1, U.size() - 1, U, i + 1, u);double Uipi = 1.0 / (U[i + p] - U[i]);double Uip1i1 = 1.0 / (U[i + p + 1] - U[i + 1]);if ((U[i + p] - U[i]) == 0)Uipi = 0.0;if ((U[i + p + 1] - U[i + 1]) == 0)Uip1i1 = 0.0;DersOnebasisfun_value = p * (Onebasisfun_i_pminus1 * Uipi - Onebasisfun_iplus1_pminus1 * Uip1i1);if (u == 1.0){u = 0.99999;double Onebasisfun_i_pminus1 = OneBasisFun_withreturn(p - 1, U.size() - 1, U, i, u);double Onebasisfun_iplus1_pminus1 = OneBasisFun_withreturn(p - 1, U.size() - 1, U, i + 1, u);double Uipi = 1.0 / (U[i + p] - U[i]);double Uip1i1 = 1.0 / (U[i + p + 1] - U[i + 1]);if ((U[i + p] - U[i]) == 0)Uipi = 0.0;if ((U[i + p + 1] - U[i + 1]) == 0)Uip1i1 = 0.0;double DersOnebasisfun_value_zero = p * (Onebasisfun_i_pminus1 * Uipi - Onebasisfun_iplus1_pminus1 * Uip1i1);DersOnebasisfun_value = -DersOnebasisfun_value_zero;//cout<< "Hehehehhe" <<  " " << i << " " << Onebasisfun_i_pminus1  << " " << Uipi  << " " <<  Onebasisfun_iplus1_pminus1  << " " << Uip1i1 << " " << DersOnebasisfun_value_zero << "\n";}}if (p == 1){if (i == 0)DersOnebasisfun_value = -1;else {DersOnebasisfun_value = 1;}}return DersOnebasisfun_value;}//将第k个单元的节点编号  对应的坐标用9*2的向量存储Eigen::MatrixXd coor(Eigen::MatrixXd Node, Eigen::MatrixXd ele, int k){// k从0开始  0对应第一个单元  1对应第二个单元...MatrixXd coor(9, 2);                     //节点坐标MatrixXd _ele = ele.block(k, 1, 1, 9);   // 将第i个单元节点编号用向量存储起来for (int i = 0; i < 9; i++){int j = _ele(i);             // j: 代表第几个节点编号coor(i, 0) = Node(j - 1, 1);coor(i, 1) = Node(j - 1, 2);}return coor;}//计算单元的雅克比矩阵Eigen::MatrixXd Jacobi(Eigen::MatrixXd Node, Eigen::MatrixXd ele, int k, vector<double> U, vector<double> V, int p, double u, double v){// Node: 节点坐标信息    ele：节点信息   k：表示的第k+1个单元// U：U方向的节点向量    p：B样条基函数的次数   // 计算雅克比矩阵 （每个单元有且只有一个雅克比矩阵）MatrixXd J(2, 2);       //雅克比矩阵J.setZero();MatrixXd _ele = ele.block(k, 1, 1, 9);   // 将第k个单元节点编号用向量存储起来MatrixXd coor1(9, 2);coor1 = coor(Node, ele, k);              //将第k个单元的节点编号  对应的坐标用9*2的向量存储MatrixXd der_eta_xi(2, 9);               // 用于存储相应基函数的导数int num;int Length_U = U.size();int Length_V = V.size();for (int i = 0; i < 9; i++){num = _ele(i);             // num : 代表第几个节点编号// 将节点编号和基函数对应起来      // 1 ---  N0,2 * M0,2   2 ---- N0,2 * M1,2   3 ---- N0,2 * M2,2// 4 ---  N1,2 * M0,2   5 ---- N1,2 * M1,2   6 ---- N1,2 * M2,2   7 ---  N2,2 * M0,2...... int val1 = (num - 1) / 3;        // val1: 表示U方向的第几个基函数   对应求导函数DersOneBasisFun_double_First_Ders(p, U, i, u)中的iint val2 = (num - 1) % 3;        // val2: 表示V方向的第几个基函数double U_xi = OneBasisFun_withreturn(p, Length_U - 1, U, val1, u);    //U方向基函数double V_eta = OneBasisFun_withreturn(p, Length_V - 1, V, val2, v);double der_U_xi = DersOneBasisFun_double_First_Ders(p, U, val1, u);   //U方向基函数导数double der_V_eta = DersOneBasisFun_double_First_Ders(p, V, val2, v);der_eta_xi(0, i) = der_U_xi * V_eta;der_eta_xi(1, i) = U_xi * der_V_eta;}J = der_eta_xi * coor1;return J;}// 计算单元应变矩阵BEigen::MatrixXd B_matrix(Eigen::MatrixXd Node, Eigen::MatrixXd ele, int k, vector<double> U, vector<double> V, int p, double u, double v){// Node: 节点坐标信息    ele：节点信息   k：表示的第k+1个单元// U：U方向的节点向量    p：B样条基函数的次数   // 计算单元应变矩阵BMatrixXd B(3, 18);B.setZero();MatrixXd J = Jacobi(Node, ele, k, U, V, p, u, v);int num;int Length_U = U.size();int Length_V = V.size();MatrixXd der_Ni_xy(2, 1);MatrixXd der_Ni_xi_eta(2, 1);for (int i = 0; i < 9; i++){MatrixXd _ele = ele.block(k, 1, 1, 9);   // 将第k个单元节点编号用向量存储起来num = _ele(i);             // num : 代表第几个节点编号// 1 ---  N0,2 * M0,2   2 ---- N0,2 * M1,2   3 ---- N0,2 * M2,2// 4 ---  N1,2 * M0,2   5 ---- N1,2 * M1,2   6 ---- N1,2 * M2,2   7 ---  N2,2 * M0,2...... int val1 = (num - 1) / 3;        // val1: 表示U方向的第几个基函数   对应求导函数DersOneBasisFun_double_First_Ders(p, U, i, u)中的iint val2 = (num - 1) % 3;        // val2: 表示V方向的第几个基函数double  U_xi = OneBasisFun_withreturn(p, Length_U - 1, U, val1, u);   //U方向基函数double V_eta = OneBasisFun_withreturn(p, Length_V - 1, V, val2, v);double  der_U_xi = DersOneBasisFun_double_First_Ders(p, U, val1, u);  //U方向基函数导数double der_V_eta = DersOneBasisFun_double_First_Ders(p, V, val2, v);double der_Ni_xi = der_U_xi * V_eta;double der_Ni_eta = U_xi * der_V_eta;                //Ni对η(eta)的导数der_Ni_xi_eta(0, 0) = der_Ni_xi;der_Ni_xi_eta(1, 0) = der_Ni_eta;der_Ni_xy = J.inverse() * der_Ni_xi_eta;            //求出Ni在物理域上的导数B(0, 2 * i) = der_Ni_xy(0, 0);B(1, 2 * i + 1) = der_Ni_xy(1, 0);B(2, 2 * i) = der_Ni_xy(1, 0);B(2, 2 * i + 1) = der_Ni_xy(0, 0);}return B;}//将节点向量中的重节点值去掉（只保存一个）vector<double> val_vector(vector<double> vec){// 因为要用到高斯积分，第一步要将每个单元的积分上下限 和节点向量的值 对应起来vector<double> val_vec;                          //新定义一个向量，用于存储节点数值val_vec.push_back(vec[0]);for (int i = 0; i < vec.size() - 1; i++){if (vec[i] != vec[i + 1]){val_vec.push_back(vec[i + 1]);          //用于存储节点向量U中  不同节点的值（重节点的值只保存一个）}}return val_vec;}// 计算单元刚度矩阵 Eigen::MatrixXd IGA_Stiffness(Eigen::MatrixXd Node, Eigen::MatrixXd ele, int k, vector<double> U, vector<double> V, int p, MatrixXd D, double t){// Node: 节点坐标信息    ele：节点信息   k：表示的第k+1个单元// U：U方向的节点向量    p：B样条基函数的次数   t:厚度// 计算单元刚度矩阵// 因为要用到高斯积分，第一步要将每个单元的积分上下限 和节点向量的值 对应起来MatrixXd Ke(18, 18);Ke.setZero();MatrixXd B(3, 18);B.setZero();MatrixXd J(2, 2);J.setZero();// 将节点向量中的重节点值去掉（重节点信息只保存一个）vector<double> val_U;val_U = val_vector(U);                // 调用函数 val_vector(vector<double> vec)   将节点向量中的重节点值去掉vector<double> val_V;val_V = val_vector(V);// 求U和V方向的积分上下限double U_a, U_b, V_a, V_b;double u, v;int num_U = k % val_U.size();         // U方向的第几个区间int num_V = k / val_U.size();U_a = val_U[num_U];                   // U方向上的积分下限（参数域区间左端点）U_b = val_U[num_U + 1];               // U方向上的积分上限（参数域区间右端点）V_a = val_V[num_V];V_b = val_V[num_V + 1];// Gauss积分   U和V方向都使用两点Gauss公式// 高斯点double Gauss[4][2];Gauss[0][0] = -0.57735;Gauss[0][1] = -0.57735;Gauss[1][0] = 0.577350;Gauss[1][1] = -0.57735;Gauss[2][0] = 0.57735;Gauss[2][1] = 0.57735;Gauss[3][0] = -0.57735;Gauss[3][1] = 0.57735;for (int j = 0; j < 4; j++){u = (U_b - U_a) * Gauss[j][0] / 2 + (U_b + U_a) / 2;v = (V_b - V_a) * Gauss[j][1] / 2 + (V_b + V_a) / 2;B = B_matrix(Node, ele, k, U, V, p, u, v);J = Jacobi(Node, ele, k, U, V, p, u, v);Ke += t * (U_b - U_a) / 2 * (V_b - V_a) / 2 * B.transpose() * D * B * J.determinant();}return Ke;}// 组装总体刚度矩阵 Eigen::MatrixXd IGA_Assembly(Eigen::MatrixXd Node, Eigen::MatrixXd ele, vector<double> U, vector<double> V, int p, MatrixXd D, double t){// Node: 节点坐标信息    ele：节点信息   k：表示的第k+1个单元// U：U方向的节点向量    p：B样条基函数的次数   t:厚度// 组装刚度矩阵int num_GK = ele.rows() * Node.rows();      //总刚的大小  2*12MatrixXd Gk(num_GK, num_GK);Gk.setZero();int num_Ke = ele.rows();                       //单刚的个数  2int num_ele_node = ele.cols() - 1;             //单刚中节点的个数  9int size_Ke = ele.rows() * (ele.cols() - 1);   //单刚矩阵的大小MatrixXd Ke(size_Ke, size_Ke);                           //单元刚度矩阵Ke.setZero();for (int k = 0; k < num_Ke; k++){Ke = IGA_Stiffness(Node, ele, k, U, V, p, D, t);     // 第k个单元刚度矩阵MatrixXd _ele = ele.block(k, 1, 1, 9);               // 将第k个单元节点编号用向量存储起来VectorXi Dof(num_ele_node * 2);                      // 定义一个  单刚中节点的个数*2  的向量(9*2)  用于组总刚for (int r = 0; r < num_ele_node; r++){Dof(2 * r) = _ele(r) * 2 - 2;Dof(2 * r + 1) = _ele(r) * 2 - 1;}for (int n1 = 0; n1 < num_ele_node * 2; n1++){for (int n2 = 0; n2 < num_ele_node * 2; n2++){Gk(Dof(n1), Dof(n2)) = Gk(Dof(n1), Dof(n2)) + Ke(n1, n2);}}}return Gk;}VectorXd Solve_1_Model(MatrixXd K, VectorXd U, VectorXd P, int len_U){//置“1”法求解方程组//K：总刚   U：初始化的节点位移向量  P:节点力for (int j = 0; j < len_U; j++){if (U(j) == 0){K.block(j, 0, 1, K.cols()).setZero();    //将总刚矩阵K的行 置“0”  K.cols()：K矩阵的列长K.block(0, j, K.rows(), 1).setZero();    //将总刚矩阵K的列 置“0”  K.rows()：K矩阵的行长K(j, j) = 1;P(j) = 0;          //节点力置“0”}}// 节点位移U = K.lu().solve(P);         //LU分解求解线性方程组K * U = Preturn U;}};

主程序文件

#include<iostream>
#include <Eigen/Dense>
using namespace Eigen;
using namespace std;
#include <vector>
#include "IGA_Plate.h"//初始物理量
double IGA_Plate::E = 10;            // 弹性模量
double IGA_Plate::mu = 0.3;          // 泊松比
double IGA_Plate::t = 0.01;          // 厚度//节点信息，第一列为节点编号，2-4列分别为x，y，z坐标
MatrixXd Node_Matrix()
{MatrixXd Node(12, 3);Node << 1, 0, 0,2, 0, 1,3, 0, 2,4, 0.5, 0,5, 0.5, 1,6, 0.5, 2,7, 1.5, 0,8, 1.5, 1,9, 1.5, 2,10, 2, 0,11, 2, 1,12, 2, 2;return Node;
}
MatrixXd IGA_Plate::Node = Node_Matrix();MatrixXd ele_Matrix()
{MatrixXd ele(2, 10);ele << 1, 1, 2, 3, 4, 5, 6, 7, 8, 9,2, 4, 5, 6, 7, 8, 9, 10, 11, 12;return ele;
}
MatrixXd IGA_Plate::ele = ele_Matrix();
// 节点向量
vector<double> IGA_Plate::U = { 0,0,0,0.5,1,1,1 };
vector<double> IGA_Plate::V = { 0,0,0,1,1,1 };int IGA_Plate::p = 2;// 输出矩阵的各个值
void Print(MatrixXd K)
{for (int j = 0; j < K.rows(); j++){for (int i = 0; i < K.cols(); i++){cout << K(j, i) << "  ";}cout << endl;}
}int main(void)
{//初始物理量IGA_Plate iga_2d;iga_2d.D_Matrix();// 组装刚度矩阵int num_GK = iga_2d.ele.rows() * iga_2d.Node.rows();      //总刚的大小  2*12MatrixXd Gk(num_GK, num_GK);Gk.setZero();Gk = iga_2d.IGA_Assembly(iga_2d.Node, iga_2d.ele, iga_2d.U, iga_2d.V, iga_2d.p, iga_2d.D, iga_2d.t);// 载荷向量 VectorXd F(num_GK);F.setZero();F(18) = 1.0 / 3;F(20) = 1.0 / 3;F(22) = 1.0 / 3;VectorXd disp(num_GK);         // 位移 displacement disp.setOnes();disp(0) = 0; disp(1) = 0; disp(2) = 0; disp(4) = 0;int len_disp = 2 * iga_2d.Node.rows();       //节点位移矩阵长度disp = iga_2d.Solve_1_Model(Gk, disp, F, len_disp);cout << "-----------分割线-------------" << endl;cout << "\n 节点位移 disp = \n" << disp << endl;
}

第二版（小修改增加通用性）

IGA_Plate.h 文件

#pragma once
#include<iostream>
#include <Eigen/Dense>
using namespace Eigen;
using namespace std;
#include <vector>class  IGA_Plate
{public://初始物理量static double E;          //弹性模量static double mu;         //泊松比static double t;          //厚度static int p;              //基函数次数static MatrixXd Node;     //节点信息矩阵，第一列为节点编号，2-4列分别为x，y（，z）坐标static MatrixXd ele;      //单元信息，第一列为单元编号，后面各列为单元上的节点编号static vector<double> U;   //节点向量static vector<double> V;public:int num_GK = ele.rows() * Node.rows();         //总刚的大小  2*12int num_Ke = ele.rows();                       //单刚的个数  2int num_ele_node = ele.cols() - 1;             //单刚中节点的个数  9int size_Ke = ele.rows() * (ele.cols() - 1);   //单刚矩阵的大小 18int num_BasisFun_V = V.size() - p - 1;         //V 方向基函数个数public:MatrixXd D_Matrix(){MatrixXd D(3, 3);D << 1, mu, 0,mu, 1, 0,0, 0, (1 - mu) / 2;D *= (E / (1 - mu * mu));return D;}MatrixXd D = D_Matrix();//求基函数的值 （函数求值）double OneBasisFun_withreturn(int p, int m, vector<double>& U, int i, double u){// p: p次B样条    m：节点向量的个数-1   // U: 节点向量[u0,u1,...,um]   i：区间段[ui ui+1,...,ui+p+1]    // u：求p次B样条基函数 Ni,p 在u点的函数值double Nip;vector<double> N(p + 1);if ((i == 0 && u == U[0]) || (i == m - p - 1 && u == U[m])){Nip = 1.0;return Nip;}if (u < U[i] || u >= U[i + p + 1]){Nip = 0.0;return Nip;}int jj, j;for (jj = 0; jj <= p; jj++) {if (u >= U[i + jj] && u < U[i + jj + 1])N[jj] = 1.0;elseN[jj] = 0.0;}double saved;int k;for (k = 1; k <= p; k++){if (N[0] == 0.0) saved = 0.0;else saved = ((u - U[i]) * N[0]) / (U[i + k] - U[i]);for (j = 0; j < p - k + 1; j++){double Uleft = U[i + j + 1];double Uright = U[i + j + k + 1];if (N[j + 1] == 0.0){N[j] = saved; saved = 0.0;}else{double temp = N[j + 1] / (Uright - Uleft);N[j] = saved + (Uright - u) * temp;saved = (u - Uleft) * temp;}}}Nip = N[0];return Nip;}//求单个基函数的导数double DersOneBasisFun_double_First_Ders(int p, vector<double>& U, int i, double u){// 原先是DersOneBasisFun_double_First_Ders(int p,int m,vector<double> &U,int i,double u,int n)// 但是我看参数 n 没有用到，就将 int n   删除了// m：节点向量的个数-1  没有用到，就将 int m   删除了// p: p次B样条     // U: 节点向量[u0,u1,...,um]   i：区间段[ui ui+1,...,ui+p+1]    // u：求p次B样条基函数 Ni,p 在u点的导数值double DersOnebasisfun_value;if (p == 3) {double Onebasisfun_i_pminus1 = OneBasisFun_withreturn(p - 1, U.size() - 1, U, i, u);double Onebasisfun_iplus1_pminus1 = OneBasisFun_withreturn(p - 1, U.size() - 1, U, i + 1, u);double Uipi = 1.0 / (U[i + p] - U[i]);double Uip1i1 = 1.0 / (U[i + p + 1] - U[i + 1]);if ((U[i + p] - U[i]) == 0)Uipi = 0.0;if ((U[i + p + 1] - U[i + 1]) == 0)Uip1i1 = 0.0;DersOnebasisfun_value = p * (Onebasisfun_i_pminus1 * Uipi - Onebasisfun_iplus1_pminus1 * Uip1i1);if (u == 1.0){u = 0.99999;double Onebasisfun_i_pminus1 = OneBasisFun_withreturn(p - 1, U.size() - 1, U, i, u);double Onebasisfun_iplus1_pminus1 = OneBasisFun_withreturn(p - 1, U.size() - 1, U, i + 1, u);double Uipi = 1.0 / (U[i + p] - U[i]);double Uip1i1 = 1.0 / (U[i + p + 1] - U[i + 1]);if ((U[i + p] - U[i]) == 0)Uipi = 0.0;if ((U[i + p + 1] - U[i + 1]) == 0)Uip1i1 = 0.0;double DersOnebasisfun_value_zero = p * (Onebasisfun_i_pminus1 * Uipi - Onebasisfun_iplus1_pminus1 * Uip1i1);DersOnebasisfun_value = -DersOnebasisfun_value_zero;//cout<< "Hehehehhe" <<  " " << i << " " << Onebasisfun_i_pminus1  << " " << Uipi  << " " <<  Onebasisfun_iplus1_pminus1  << " " << Uip1i1 << " " << DersOnebasisfun_value_zero << "\n";}}if (p == 2) {double Onebasisfun_i_pminus1 = OneBasisFun_withreturn(p - 1, U.size() - 1, U, i, u);double Onebasisfun_iplus1_pminus1 = OneBasisFun_withreturn(p - 1, U.size() - 1, U, i + 1, u);double Uipi = 1.0 / (U[i + p] - U[i]);double Uip1i1 = 1.0 / (U[i + p + 1] - U[i + 1]);if ((U[i + p] - U[i]) == 0)Uipi = 0.0;if ((U[i + p + 1] - U[i + 1]) == 0)Uip1i1 = 0.0;DersOnebasisfun_value = p * (Onebasisfun_i_pminus1 * Uipi - Onebasisfun_iplus1_pminus1 * Uip1i1);if (u == 1.0){u = 0.99999;double Onebasisfun_i_pminus1 = OneBasisFun_withreturn(p - 1, U.size() - 1, U, i, u);double Onebasisfun_iplus1_pminus1 = OneBasisFun_withreturn(p - 1, U.size() - 1, U, i + 1, u);double Uipi = 1.0 / (U[i + p] - U[i]);double Uip1i1 = 1.0 / (U[i + p + 1] - U[i + 1]);if ((U[i + p] - U[i]) == 0)Uipi = 0.0;if ((U[i + p + 1] - U[i + 1]) == 0)Uip1i1 = 0.0;double DersOnebasisfun_value_zero = p * (Onebasisfun_i_pminus1 * Uipi - Onebasisfun_iplus1_pminus1 * Uip1i1);DersOnebasisfun_value = -DersOnebasisfun_value_zero;//cout<< "Hehehehhe" <<  " " << i << " " << Onebasisfun_i_pminus1  << " " << Uipi  << " " <<  Onebasisfun_iplus1_pminus1  << " " << Uip1i1 << " " << DersOnebasisfun_value_zero << "\n";}}if (p == 1){if (i == 0)DersOnebasisfun_value = -1;else {DersOnebasisfun_value = 1;}}return DersOnebasisfun_value;}//将第k个单元的节点编号  对应的坐标用9*2的向量存储Eigen::MatrixXd coor(Eigen::MatrixXd Node, Eigen::MatrixXd ele, int k){// k从0开始  0对应第一个单元  1对应第二个单元...MatrixXd coor(num_ele_node, Node.cols()-1);               // 每个单元的节点坐标存储为一个矩阵   节点坐标二维 Node.cols()-1 = 2MatrixXd _ele = ele.block(k, 1, 1, num_ele_node);         // 将第k个单元节点编号用向量存储起来for (int i = 0; i < num_ele_node; i++){int j = _ele(i);             // j: 代表第几个节点编号coor(i, 0) = Node(j - 1, 1);coor(i, 1) = Node(j - 1, 2);}return coor;}//将节点向量中的重节点值去掉（只保存一个）vector<double> val_vector(vector<double> vec){// 因为要用到高斯积分，第一步要将每个单元的积分上下限 和节点向量的值 对应起来vector<double> val_vec;                          //新定义一个向量，用于存储节点数值val_vec.push_back(vec[0]);for (int i = 0; i < vec.size() - 1; i++){if (vec[i] != vec[i + 1]){val_vec.push_back(vec[i + 1]);          //用于存储节点向量U中  不同节点的值（重节点的值只保存一个）}}return val_vec;}//计算单元的雅克比矩阵Eigen::MatrixXd Jacobi(Eigen::MatrixXd Node, Eigen::MatrixXd ele, int k, vector<double> U, vector<double> V, int p, double u, double v){// Node: 节点坐标信息    ele：节点信息   k：表示的第k+1个单元// U：U方向的节点向量    p：B样条基函数的次数   // 计算雅克比矩阵 （每个单元有且只有一个雅克比矩阵）MatrixXd J(2, 2);       //雅克比矩阵J.setZero();MatrixXd _ele = ele.block(k, 1, 1, num_ele_node);   // 将第k个单元节点编号用向量存储起来MatrixXd coor1(num_ele_node, 2);coor1 = coor(Node, ele, k);                         //将第k个单元的节点编号  对应的坐标用9*2的向量存储MatrixXd der_eta_xi(2, 9);                          // 用于存储相应基函数的导数int num;int Length_U = U.size();int Length_V = V.size();for (int i = 0; i < num_ele_node; i++){num = _ele(i);             // num : 代表k单元中第i个节点编号// 将节点编号和基函数对应起来      // 1 ---  N0,2 * M0,2   2 ---- N0,2 * M1,2   3 ---- N0,2 * M2,2// 4 ---  N1,2 * M0,2   5 ---- N1,2 * M1,2   6 ---- N1,2 * M2,2   7 ---  N2,2 * M0,2...... // val1: 表示U方向的第几个基函数   对应求导函数DersOneBasisFun_double_First_Ders(p, U, i, u)中的i// val_U.size()int val1 = (num - 1) / num_BasisFun_V;int val2 = (num - 1) % num_BasisFun_V;        // val2: 表示V方向的第几个基函数double U_xi = OneBasisFun_withreturn(p, Length_U - 1, U, val1, u);    //U方向基函数double V_eta = OneBasisFun_withreturn(p, Length_V - 1, V, val2, v);double der_U_xi = DersOneBasisFun_double_First_Ders(p, U, val1, u);   //U方向基函数导数double der_V_eta = DersOneBasisFun_double_First_Ders(p, V, val2, v);der_eta_xi(0, i) = der_U_xi * V_eta;der_eta_xi(1, i) = U_xi * der_V_eta;}J = der_eta_xi * coor1;return J;}// 计算单元应变矩阵BEigen::MatrixXd B_matrix(Eigen::MatrixXd Node, Eigen::MatrixXd ele, int k, vector<double> U, vector<double> V, int p, double u, double v){// Node: 节点坐标信息    ele：节点信息   k：表示的第k+1个单元// U：U方向的节点向量    p：B样条基函数的次数   // 计算单元应变矩阵BMatrixXd B(3, 2 * num_ele_node);B.setZero();MatrixXd J = Jacobi(Node, ele, k, U, V, p, u, v);int num;int Length_U = U.size();int Length_V = V.size();MatrixXd der_Ni_xy(2, 1);MatrixXd der_Ni_xi_eta(2, 1);for (int i = 0; i < num_ele_node; i++){MatrixXd _ele = ele.block(k, 1, 1, num_ele_node);   // 将第k个单元节点编号用向量存储起来num = _ele(i);             // num : 代表第几个节点编号// 1 ---  N0,2 * M0,2   2 ---- N0,2 * M1,2   3 ---- N0,2 * M2,2// 4 ---  N1,2 * M0,2   5 ---- N1,2 * M1,2   6 ---- N1,2 * M2,2   7 ---  N2,2 * M0,2...... int val1 = (num - 1) / num_BasisFun_V;        // val1: 表示U方向的第几个基函数   对应求导函数DersOneBasisFun_double_First_Ders(p, U, i, u)中的iint val2 = (num - 1) % num_BasisFun_V;        // val2: 表示V方向的第几个基函数double  U_xi = OneBasisFun_withreturn(p, Length_U - 1, U, val1, u);   //U方向基函数double V_eta = OneBasisFun_withreturn(p, Length_V - 1, V, val2, v);double  der_U_xi = DersOneBasisFun_double_First_Ders(p, U, val1, u);  //U方向基函数导数double der_V_eta = DersOneBasisFun_double_First_Ders(p, V, val2, v);double der_Ni_xi = der_U_xi * V_eta;double der_Ni_eta = U_xi * der_V_eta;                //Ni对η(eta)的导数der_Ni_xi_eta(0, 0) = der_Ni_xi;der_Ni_xi_eta(1, 0) = der_Ni_eta;der_Ni_xy = J.inverse() * der_Ni_xi_eta;            //求出Ni在物理域上的导数B(0, 2 * i) = der_Ni_xy(0, 0);B(1, 2 * i + 1) = der_Ni_xy(1, 0);B(2, 2 * i) = der_Ni_xy(1, 0);B(2, 2 * i + 1) = der_Ni_xy(0, 0);}return B;}// 计算单元刚度矩阵 Eigen::MatrixXd IGA_Stiffness(Eigen::MatrixXd Node, Eigen::MatrixXd ele, int k, vector<double> U, vector<double> V, int p, MatrixXd D, double t){// Node: 节点坐标信息    ele：节点信息   k：表示的第k+1个单元// U：U方向的节点向量    p：B样条基函数的次数   t:厚度// 计算单元刚度矩阵// 因为要用到高斯积分，第一步要将每个单元的积分上下限 和节点向量的值 对应起来MatrixXd Ke(size_Ke, size_Ke);Ke.setZero();MatrixXd B(3, 2 * num_ele_node);B.setZero();MatrixXd J(2, 2);J.setZero();// 将节点向量中的重节点值去掉（重节点信息只保存一个）vector<double> val_U;val_U = val_vector(U);                // 调用函数 val_vector(vector<double> vec)   将节点向量中的重节点值去掉vector<double> val_V;val_V = val_vector(V);// 求U和V方向的积分上下限double U_a, U_b, V_a, V_b;double u, v;int num_U = k % val_U.size();         // U方向的第几个区间int num_V = k / val_U.size();U_a = val_U[num_U];                   // U方向上的积分下限（参数域区间左端点）U_b = val_U[num_U + 1];               // U方向上的积分上限（参数域区间右端点）V_a = val_V[num_V];V_b = val_V[num_V + 1];// Gauss积分   U和V方向都使用两点Gauss公式// 高斯点double Gauss[4][2];Gauss[0][0] = -0.57735;Gauss[0][1] = -0.57735;Gauss[1][0] = 0.577350;Gauss[1][1] = -0.57735;Gauss[2][0] = 0.57735;Gauss[2][1] = 0.57735;Gauss[3][0] = -0.57735;Gauss[3][1] = 0.57735;for (int j = 0; j < 4; j++){u = (U_b - U_a) * Gauss[j][0] / 2 + (U_b + U_a) / 2;v = (V_b - V_a) * Gauss[j][1] / 2 + (V_b + V_a) / 2;B = B_matrix(Node, ele, k, U, V, p, u, v);J = Jacobi(Node, ele, k, U, V, p, u, v);Ke += t * (U_b - U_a) / 2 * (V_b - V_a) / 2 * B.transpose() * D * B * J.determinant();}return Ke;}// 组装总体刚度矩阵 Eigen::MatrixXd IGA_Assembly(){// Node: 节点坐标信息    ele：节点信息   k：表示的第k+1个单元// U：U方向的节点向量    p：B样条基函数的次数   t:厚度// 组装刚度矩阵MatrixXd Gk(num_GK, num_GK);Gk.setZero();MatrixXd Ke(size_Ke, size_Ke);                           //单元刚度矩阵Ke.setZero();for (int k = 0; k < num_Ke; k++){Ke = IGA_Stiffness(Node, ele, k, U, V, p, D, t);     // 第k个单元刚度矩阵MatrixXd _ele = ele.block(k, 1, 1, 9);               // 将第k个单元节点编号用向量存储起来VectorXi Dof(num_ele_node * 2);                      // 定义一个  单刚中节点的个数*2  的向量(9*2)  用于组总刚for (int r = 0; r < num_ele_node; r++){Dof(2 * r) = _ele(r) * 2 - 2;Dof(2 * r + 1) = _ele(r) * 2 - 1;}for (int n1 = 0; n1 < num_ele_node * 2; n1++){for (int n2 = 0; n2 < num_ele_node * 2; n2++){Gk(Dof(n1), Dof(n2)) = Gk(Dof(n1), Dof(n2)) + Ke(n1, n2);}}}return Gk;}VectorXd Solve_1_Model(MatrixXd K, VectorXd U, VectorXd P, int len_U){//置“1”法求解方程组//K：总刚   U：初始化的节点位移向量  P:节点力for (int j = 0; j < len_U; j++){if (U(j) == 0){K.block(j, 0, 1, K.cols()).setZero();    //将总刚矩阵K的行 置“0”  K.cols()：K矩阵的列长K.block(0, j, K.rows(), 1).setZero();    //将总刚矩阵K的列 置“0”  K.rows()：K矩阵的行长K(j, j) = 1;P(j) = 0;          //节点力置“0”}}// 节点位移U = K.lu().solve(P);         //LU分解求解线性方程组K * U = Preturn U;}};

主程序文件

#include<iostream>
#include <Eigen/Dense>
using namespace Eigen;
using namespace std;
#include <vector>
#include "IGA_Plate.h"//初始物理量
double IGA_Plate::E = 10;            // 弹性模量
double IGA_Plate::mu = 0.3;          // 泊松比
double IGA_Plate::t = 0.01;          // 厚度//节点信息，第一列为节点编号，2-4列分别为x，y，z坐标
MatrixXd Node_Matrix()
{MatrixXd Node(12, 3);Node << 1, 0, 0,2, 0, 1,3, 0, 2,4, 0.5, 0,5, 0.5, 1,6, 0.5, 2,7, 1.5, 0,8, 1.5, 1,9, 1.5, 2,10, 2, 0,11, 2, 1,12, 2, 2;return Node;
}
MatrixXd IGA_Plate::Node = Node_Matrix();MatrixXd ele_Matrix()
{MatrixXd ele(2, 10);ele << 1, 1, 2, 3, 4, 5, 6, 7, 8, 9,2, 4, 5, 6, 7, 8, 9, 10, 11, 12;return ele;
}
MatrixXd IGA_Plate::ele = ele_Matrix();// 节点向量
vector<double> IGA_Plate::U = { 0,0,0,0.5,1,1,1 };
vector<double> IGA_Plate::V = { 0,0,0,1,1,1 };int IGA_Plate::p = 2;// 输出矩阵的各个值
void Print(MatrixXd K)
{for (int j = 0; j < K.rows(); j++){for (int i = 0; i < K.cols(); i++){cout << K(j, i) << "  ";}cout << endl;}
}int main(void)
{//初始物理量IGA_Plate iga_2d;iga_2d.D_Matrix();// 组装刚度矩阵int num_GK = iga_2d.ele.rows() * iga_2d.Node.rows();      //总刚的大小  2*12MatrixXd Gk(num_GK, num_GK);Gk.setZero();Gk = iga_2d.IGA_Assembly();cout << "-----------分割线4-------------" << endl;// 载荷向量 VectorXd F(num_GK);F.setZero();F(18) = 1.0 / 3;F(20) = 1.0 / 3;F(22) = 1.0 / 3;VectorXd disp(num_GK);         // 位移 displacement disp.setOnes();disp(0) = 0; disp(1) = 0; disp(2) = 0; disp(4) = 0;int len_disp = 2 * iga_2d.Node.rows();       //节点位移矩阵长度disp = iga_2d.Solve_1_Model(Gk, disp, F, len_disp);cout << "-----------分割线3-------------" << endl;cout << "\n 节点位移 disp = \n" << disp << endl;cout << "-----------分割线6-------------" << endl;cout << iga_2d.val_vector(iga_2d.U).size() << endl;cout << iga_2d.U.size() << endl;
}

平面问题IGA程序初稿相关推荐

空间刚架matlab_基本平面刚架MATLAB程序
% 平面刚架 MATLAB 程序 % 2003.9.16 2007.2.28 2008.4.1 2009.10 2011.10 2013.9 2014.09 %******************** ...
空间刚架matlab_2016基本平面刚架各种荷载MATLAB程序
% 平面刚架 MATLAB 程序 % 2003.9.16 2007.2.28 2008.4.1 2009.10 2011.10 2013.9 2014.09 2016.03 %************ ...
为特使建立控制平面的指南-部署权衡
部署控制平面组件构建并设计了控制平面后,您将需要确切确定如何部署其组件. 在这里,您可以选择将控制平面与数据平面共置一处以集中控制平面. 这里还有一个中间立场:部署与控制平面位于同一位置的某些组件, ...
matlab三角形单元,平面三角形单元常应变单元matlab程序的编制.doc
平面三角形单元常应变单元matlab程序的编制.doc 1三角形常应变单元程序的编制与使用有限元法是求解微分方程边值问题的一种通用数值方法,该方法是一种基于变分法(或变分里兹法)而发展起来的求解微分方 ...
5800计算器公路三维全能程序
5800计算器公路三维全能程序.(丢掉图纸轻松测量)说明清晰! 各位,我的9860程序记忆编写完毕,跟这个程序功能差不多,但是更好用,加入了隧道计算功能,能计算超欠挖,渐变隧道,不限制圆心数目,程序 ...
精简5800三维程序
精简5800三维程序正反算选择程序:ZS-FS Deg: "1LZ＝>XY,2XY＝>LZ"?U: //正反算选择,正算选1, ...
精简5800三维程序]
精简5800三维程序正反算选择程序:ZS-FS Deg: "1LZ＝>XY,2XY＝>LZ"?U: //正反算选择,正算选1, ...
平板 matlab,MATLAB实现平板平面度数据处理
计算机光盘软件与应用 2011年第16期工程技术 Computer CD Software and Applications MATLAB实现平板平而度数据处理李明贵 (贵州省机械电子产品质量监督 ...
网内计算：可编程数据平面和技术特定应用综述
网内计算:可编程数据平面和技术特定应用综述摘要--与云计算相比,边缘计算提供了更靠近终端设备的处理,降低了用户体验的延迟.最新的In-Network Computing范例采用可编程网络元素在数据达 ...

平面问题IGA程序初稿

readme

问题描述

编程遇到的问题

结果

程序

IGA_Plate.h 文件

主程序文件

第二版（小修改增加通用性）

IGA_Plate.h 文件

主程序文件

平面问题IGA程序初稿相关推荐

最新文章

热门文章