前言

只是为方便学习，不做其他用途。
时间：2023.02.20
细化后的结果和仿真结果相比误差更大了
将老师给的B样条基函数换成了gismo库中的基函数（本来想的是可能老师的基函数程序出问题了，结果是我想多了），两个单元的等几何结果一致
本文用于存储：用gismo库中的基函数处理等几何（没有细化）

一、需要添加gismo库程序才能运行

1.1 主程序文件

#include<iostream>
#include <gismo.h>
#include <Eigen/Dense>
using namespace Eigen;
using namespace gismo;
using namespace std;
#include <vector>
#include "nurbs_uniformRefine.h"
#include "IGA_Plate.h"//初始物理量
int IGA_Plate::p = 2;
double IGA_Plate::E = 10;            // 弹性模量
double IGA_Plate::mu = 0.3;          // 泊松比
double IGA_Plate::t = 0.01;          // 厚度// 细化前的初始物理量
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 k_Refine::pre_Node = Node_Matrix();           //节点信息矩阵
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::Node = Node_Matrix();
MatrixXd IGA_Plate::ele = ele_Matrix();       //注意：单元编号是整型数据  数据类型要统一
// 数据类型最后我都统一改成了double型了  本来 ele 想要用 int型， 发现得改好多之前写的IGA里面的东西，就用int型了gsKnotVector<> IGA_Plate::U(0, 1, 1, 3);           //节点向量
gsKnotVector<> IGA_Plate::V(0, 1, 0, 3);// 输出向量的各个值
void Print_Vector(vector<double> U)
{cout << "  向量值 " << " = " << endl;for (int i = 0; i < U.size(); i++){cout << "  " << U[i] << "  " ;}cout << endl;
}
// 输出矩阵的各个值
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()
{IGA_Plate iga_2d;// 组装刚度矩阵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();// 载荷向量 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(3) = 0; disp(5) = 0; disp(4) = 0;int len_disp = 2 * iga_2d.Node.rows();        //节点位移矩阵长度disp = iga_2d.Solve_1_Model(Gk, disp, F);        //置“1”法求解方程组cout << "\n 节点位移 disp = \n" << disp << endl;}

1.2 头文件 “IGA_Plate.h”

#pragma once
#include<iostream>
#include <gismo.h>
#include <Eigen/Dense>
using namespace Eigen;
using namespace std;
using namespace gismo;
#include <vector>class  IGA_Plate
{public://初始物理量static double E, mu, t;          //弹性模量  泊松比  厚度static int p;                    //基函数次数static MatrixXd Node;            //节点信息矩阵static MatrixXd ele;             //单元信息static gsKnotVector<> U, V;      //节点向量public:int size_GK = 2 * Node.rows();                  //总刚的大小  2*控制点个数（Node行数）int num_Ke = ele.rows();                       //单刚的个数  int num_ele_node = ele.cols();                 //单刚中节点的个数  9int size_Ke = 2 * (ele.cols());                //单刚矩阵的大小 18  不变int num_BasisFun_V = V.size() - p - 1;         //V 方向基函数个数  节点个数int num_BasisFun_U = U.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(gsKnotVector<> kv, int i, double u){// i:表示生成的第几个基函数gsBSplineBasis<double>base(kv);     // 生成单参数B样条基函数gsMatrix< double > val_u(1, 1);val_u << u;gsMatrix<double> result = base.evalSingle(i, val_u);    // Nip在u点的值double Nip = result(0, 0);return Nip;}//求单个基函数的导数double DersOneBasisFun_double_First_Ders(gsKnotVector<> kv, int i, double u){// i:表示生成的第几个基函数gsBSplineBasis<double>base(kv);     // 生成单参数B样条基函数gsMatrix< double > val_u(1, 1);val_u << u;gsMatrix<double> result = base.derivSingle(i, val_u);    // Nip在u点的导数值double der_Nip = result(0, 0);return der_Nip;}//将第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);                          // 用于存储相应基函数的导数der_eta_xi.setZero();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)中的iint val1 = (num - 1) / num_BasisFun_V;int val2 = (num - 1) % num_BasisFun_V;        // val2: 表示V方向的第几个基函数double U_xi = OneBasisFun_withreturn(U, val1, u);     //U方向基函数double V_eta = OneBasisFun_withreturn(V, val2, v);double der_U_xi = DersOneBasisFun_double_First_Ders(U, val1, u);    //U方向基函数导数double der_V_eta = DersOneBasisFun_double_First_Ders(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(int k, double u, double v){// Nk：表示的第k+1个单元// 计算单元应变矩阵BMatrixXd 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)中的iint val2 = (num - 1) % num_BasisFun_V;        // val2: 表示V方向的第几个基函数double  U_xi = OneBasisFun_withreturn(U, val1, u);    //U方向基函数double V_eta = OneBasisFun_withreturn(V, val2, v);double  der_U_xi = DersOneBasisFun_double_First_Ders(U, val1, u);   //U方向基函数导数double der_V_eta = DersOneBasisFun_double_First_Ders(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;// num_BasisFun_V - 2:  V方向单元个数int num_U = k / (num_BasisFun_V - 2);         // U方向的第几个区间int num_V = k % (num_BasisFun_V - 2);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(size_GK, size_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个单元节点编号用向量存储起来VectorXd 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++){int a = Dof(n1);int b = Dof(n2); Gk(a, b) = Gk(a, b) + Ke(n1, 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){// 置“1”法求解方程组// K：总刚   U：初始化的节点位移向量  P:节点力int len_U = U.rows();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;}
};

1.3 结果

原先程序计算结果和调用gismo中基函数结果一致（施加载荷都如下图所示）：

师弟用其他软件给我仿真的边界是前三个点位移不动，所以位移边界一开始我给错了，应该为下图：

用gismo中B样条基函数替代自己写的基函数相关推荐

在3dmax软件中添加样条的方法和详细步骤
在3dmax软件中添加样条的方法和详细步骤! 在3dmax软件中添加样条的方法和详细步骤!三dsMax是一款三建模.动画和渲染软件.借助3dsMax,可以创造一个宏伟的游戏世界,布彩的场景,实现设计可 ...
python字符串替换源码_Python实现字符串中某个字母的替代功能
Python实现字符串中某个字母的替代功能今晚想实现这样一个功能:将输入字符串中的字母 "i" 变成字母 "p".当时想的很简单,直接用for循环遍历,然后替 ...
【Kotlin】Kotlin 中使用 Lambda 表达式替代对象表达式原理分析 ( 尾随 Lambda - Trailing Lambda 语法 | 接口对象表达式 = 接口#函数类型对象 )
文章目录一.尾随 Lambda - Trailing Lambda 语法二.Kotlin 中使用 Lambda 表达式替代对象表达式原理 1.Lambda 替换对象表达式 2.原理分析 3.示例分 ...
替代微软txt文本编辑器_如何在Microsoft Excel中向对象添加替代文本
替代微软txt文本编辑器 Alternative text (alt text) allows screen readers to capture the description of an obje ...
如何在Microsoft Word中向对象添加替代文本
Alternative text (alt text) allows screen readers to capture the description of an object and read i ...
谈谈我对前端组件化中“组件”的理解，顺带写个Vue与React的demo
谈谈我对前端组件化中"组件"的理解,顺带写个Vue与React的demo 前言前端已经过了单兵作战的时代了,现在一个稍微复杂一点的项目都需要几个人协同开发,一个战略级别的APP的 ...
python循环语句-python中的for循环语句怎么写
python中的for循环语句怎么写? Python for 循环语句 Python for循环可以遍历任何序列的项目,如一个列表或者一个字符串. for循环的语法格式如下: 1 2 for iter ...
idea中构造器和toString方法覆写的快捷键
idea中构造器和toString方法覆写的快捷键写上一篇博文已经是3个月前的事情了,这中间因为各种杂事耽搁,好久没有学习我的Java了.要利用这个暑假好好的学习一波了. 废话不多说,进入主题. 老 ...
科技论文中的分析与综合-如何写好科技论文之我见(七）
科技论文中的分析与综合-----如何写好科技论文之我见(七) 闵应骅分析与综合这两术语大家经常用.但是,真要说它们的定义,那可是哲学范围里的事.形式逻辑里面就有分析与综合.我在初中教几何的时候,就常 ...

用gismo中B样条基函数替代自己写的基函数

文章目录

前言

一、需要添加gismo库程序才能运行

1.1 主程序文件

1.2 头文件 “IGA_Plate.h”

1.3 结果

用gismo中B样条基函数替代自己写的基函数相关推荐

最新文章

热门文章