IGA程序副本-2023.02.20

---

前言

只是为保存，怕之后出错，找不到原程序

一、主程序

#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"

// 细化前的初始物理量
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 k_Refine::pre_ele = ele_Matrix();            //单元信息
gsKnotVector<> k_Refine::pre_U(0, 1, 1, 3);           //节点向量
gsKnotVector<> k_Refine::pre_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;
}

k_Refine span_2_3;   //  定义一个类对象

//初始物理量
int IGA_Plate::p = 2;
double IGA_Plate::E = 10;            // 弹性模量
double IGA_Plate::mu = 0.3;          // 泊松比
double IGA_Plate::t = 0.01;          // 厚度

MatrixXd IGA_Plate::Node = span_2_3.Node;
MatrixXd IGA_Plate::ele = span_2_3.ele;       //注意：单元编号是整型数据  数据类型要统一
// 数据类型最后我都统一改成了double型了  本来 ele 想要用 int型， 发现得改好多之前写的IGA里面的东西，就用int型了

// 节点向量
vector<double> IGA_Plate::U = span_2_3.U;
vector<double> IGA_Plate::V = span_2_3.V;

// 等效节点载荷
VectorXd  return_F(IGA_Plate iga_2d)
{
   
    // 等效节点载荷
    int num_Point = iga_2d.num_BasisFun_V * iga_2d.num_BasisFun_U;  // 节点总个数
    int Len_F = 2 * num_Point;     // x 和 y 两个方向的力
    VectorXd F(Len_F);
    F.setZero();
    for (int j = 0; j < iga_2d.num_BasisFun_V; j++)
    {
   
        F(Len_F - 2 - 2 * j) = 1.0 / iga_2d.num_BasisFun_V;
    }
    //cout << " num_Point = " << num_Point << endl;
    //gsInfo << "\n 等效节点力 F = \n" << F << endl;
    return F;
}
// 对边界的处理  位移不变节点的处理   用向量存储，不变的点存储0，位移改变的节点存储1
VectorXd  return_disp(MatrixXd BC, IGA_Plate iga_2d)
{
   
    int size_GK = iga_2d.size_GK;
    VectorXd disp(size_GK);         // 位移 displacement 
    disp.setOnes();
    for (int i = 0; i < BC.rows(); i++)
    {
   
        int num_node = BC(i, 0);      //  第几个节点
        int type_node = BC(i, 1);     //  该节点的 x方向还是y方向位移不变  1：表示x方向   2：表示y方向
        if (type_node == 1)
        {
   
            disp(2 * num_node - 2) = 0;
        }
        if (type_node == 2)
        {
   
            disp(2 * num_node - 1) = 0;
        }
        else
        {
   
            //cout << " 输入数据有误！！！" << endl;
        }
    }
    return disp;
}

int k_Refine::k = 0;    // 细化次数
int main()
{
   
    IGA_Plate iga_2d;

    // 组装刚度矩阵
    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 = return_F(iga_2d);
    gsInfo << "\n 等效节点力 F = \n" << F << endl;

    // 位移边界条件  位移不变的点的处理
    MatrixXd BC = span_2_3.BC;
    VectorXd disp = return_disp(BC, iga_2d);
    
    cout << "\n 节点位移不变量 disp = \n" << BC << endl;
    cout << "\n 节点位移不变量 disp = \n" << disp << endl;

    disp = iga_2d.Solve_1_Model(Gk, disp, F);        //置“1”法求解方程组
    

    cout << "\n 节点位移 disp = \n" << disp << endl;
    double max = disp.maxCoeff();
    cout << "\n 最大节点位移 disp = \n" << max << 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, mu, t;          //弹性模量  泊松比  厚度
	static int p;                    //基函数次数
	static MatrixXd Node;            //节点信息矩阵
	static MatrixXd ele;             //单元信息
	static vector<double> U, V;      //节点向量

public:
	int size_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 <