平面问题IGA程序初稿有详细对问题的描述
程序
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, mu, t; //弹性模量 泊松比 厚度
static int p; //基函数次数
static MatrixXd Node; //节点信息矩阵，第一列为节点编号，2-4列分别为x，y（，z）坐标
static MatrixXd ele; //单元信息，第一列为单元编号，后面各列为单元上的节点编号
static vector<double> U, V; //节点向量
public:
int num_GK = 2 * Node.rows(); //总刚的大小 2*12
int num_Ke = ele.rows(); //单刚的个数 2
int num_ele_node = ele.cols(); //单刚中节点的个数 9
int size_Ke = 2 * (ele.cols()); //单刚矩阵的大小 18
int 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;
else
N[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(int k)
{
// k从0开始 0对应第一个单元 1对应第二个单元...
MatrixXd coor(num_ele_node, Node.cols()); // 每个单元的节点坐标存储为一个矩阵 节点坐标二维 Node.cols()
MatrixXd _ele = ele.block(k, 0, 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, 0);
coor(i, 1) = Node(j - 1, 1);
}
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(int k, double u, double v)
{
// k：表示的第k+1个单元
// 计算雅克比矩阵 （每个单元有且只有一个雅克比矩阵）
MatrixXd J(2, 2); //雅克比矩阵
J.setZero();
MatrixXd _ele = ele.block(k, 0, 1, num_ele_node); // 将第k个单元节点编号用向量存储起来
MatrixXd coor1(num_ele_node, 2);
coor1 = coor(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
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;
}
// 计算单元应变矩阵B
Eigen::MatrixXd B_matrix(int k, double u, double v)
{
// Nk：表示的第k+1个单元
// 计算单元应变矩阵B
MatrixXd B(3, 2 * num_ele_node);
B.setZero();
MatrixXd J = Jacobi(k, 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, 0, 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)中的i
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);
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(int k)
{
// t:厚度 k：表示的第k+1个单元
// 计算单元刚度矩阵
// 因为要用到高斯积分，第一步要将每个单元的积分上下限 和节点向量的值 对应起来
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(k, u, v);
J = Jacobi(k, 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(k); // 第k个单元刚度矩阵
MatrixXd _ele = ele.block(k, 0, 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;
}
//置“1”法求解方程组
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 = P
return 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, 2);
Node << 0, 0,
0, 1,
0, 2,
0.5, 0,
0.5, 1,
0.5, 2,
1.5, 0,
1.5, 1,
1.5, 2,
2, 0,
2, 1,
2, 2;
return Node;
}
MatrixXd IGA_Plate::Node = Node_Matrix();
// 输出矩阵的各个值
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;
}
}
MatrixXd ele_Matrix()
{
MatrixXd ele(2, 9);
ele << 1, 2, 3, 4, 5, 6, 7, 8, 9,
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;
int main(void)
{
//初始物理量
IGA_Plate iga_2d;
Print(iga_2d.coor(0));
// 组装刚度矩阵
int num_GK = iga_2d.ele.rows() * iga_2d.Node.rows(); //总刚的大小 2*12
MatrixXd Gk(num_GK, num_GK);
Gk.setZero();
Gk = iga_2d.IGA_Assembly();
// 载荷向量
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); //置“1”法求解方程组
cout << "\n 节点位移 disp = \n" << disp << endl;
}
7336
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
