等几何分析编程记录 --- 未完待续

---

前言

记录代码

eigen库中常用代码

矩阵置0 求逆

MatrixXd B(3, 18);
	B.setZero(); //矩阵置0  
B.inverse()  //矩阵求逆

一、调用gismo细化后生成新的ele单元

ele单元：每个单元的全局控制点编号
gismo细化：节点插入（具体查看gismo库，也可以查看gismo中NURBS的相关函数的使用）

1.1 理论分析

以两个单元的双二次张量基型B样条平面为例：

1.2 用一个单元的双二次曲面验证一下

：

两次加密后

1.3 程序和结果

程序：

#include<iostream>
#include <gismo.h>
#include <Eigen/Dense>
using namespace Eigen;
using namespace gismo;
using namespace std;
#include <vector>

//将节点向量中的重节点值去掉（只保存一个）
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;
}

int main()
{
   
    // Nurbs module   网址： https://gismo.github.io/group__Nurbs.html
    // 该模块包含b样条和NURBS基函数
    // G+Smo在处理张量积型B样条方面非常强大:
    // 例1：
    // 建立一个二元二次张量积基，并计算了几个点上的基函数  首先，我们需要节点向量:
    gsKnotVector<> kv0(0, 1, 1, 3);    //节点的左端点  右端点  内部节点个数  左/右端点的节点重数   
    gsKnotVector<> kv1(0, 1, 0, 3);    //start,end,interior knots, start/end multiplicites of knots

    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;

    // 生成张量基型Nurbs曲面
    gsTensorNurbs<2, double> nurbs(kv1, kv0, Node);       // nurbs(V,U,Node)  先V方向上节点向量  后U方向节点向量
    // 调用细化函数
    nurbs.uniformRefine();
    
    //  细化后调用节点向量
    vector<double> U = nurbs.knots(1);       // 细化后U方向节点向量 
    vector<double> V = nurbs.knots(0);       // 细化后V方向节点向量

    gsInfo << nurbs.knots(1);   // 细化后V方向节点向量
    cout << endl;
    gsInfo << nurbs.knots(0);   // 细化后U方向节点向量
    cout << "---------分割线----------" << endl;
    //  细化后调用节点坐标
    Node = nurbs.coefs();            //  细化后调用节点坐标矩阵
    cout << "控制点坐标为: " << endl;
    gsInfo << nurbs.coefs() << endl;

    // 将节点向量中的重节点值去掉（只保存一个）
    vector<double> val_U = val_vector(U);
    vector<double> val_V = val_vector(V);

    // 基函数个数、控制点个数
    int p = 2;                 //  基函数次数
    int num_BasisFun_U = U.size() - p - 1;      // U方向方向基函数个数 = U方向控制点个数
    int num_of_point_U = U.size() - p - 1;      // U方向控制点个数
    int num_BasisFun_V = V.size() - p - 1;
    int num_of_point_V = V.size() - p - 1;

    // 单元个数
    int num_of_ele_U = val_U.size() - 1;             // U方向区间个数 = U方向单元个数
    int num_of_ele_V = val_V.size() - 1;
    int num_of_ele = num_of_ele_U * num_of_ele_V;    // 细化后的单元总数

    // 生成每个单元对应的全局编号
    // 通过使用一个矩阵，矩阵的对应位置的值 = 全局编号的值，更好的处理划分单元的全局编号（好像有点说不清楚）
    MatrixXi A(num_of_point_V, num_of_point_U);
    int sum = -1;
    for (int i = 0; i < num_of_point_U; i++)
    {
   
        for (int j = 0; j < num_of_point_V; j++)
        {
   
            sum += 1;
            A(j, i) = sum;
        }
    }
    // 单元信息
    MatrixXi ele(num_of_ele, 9);
    ele.setZero();
    for (int k = 0; k < num_of_ele; k++)
    {
   
        int i = k % num_of_ele_V;
        int j = k / num_of_ele_V;
        ele.row(k) << A(i, j), A(i + 1, j), A(i + 2, j),
            A(i, j + 1), A(i + 1, j + 1), A(i + 2, j + 1),
            A(i, j + 2), A(i + 1, j + 2), A(i + 2, j + 2);
    }

    cout << "---------分割线----------" << endl;
    gsInfo << "A = "<< endl << A << endl;

    gsInfo << "ele = " << endl << ele << endl;
    
   

    cout << "---------分割线----------" << endl;
    cout << num_of_point_U << endl;


}

结果：

二、将细化后的单元信息和节点信息代入IGA程序中遇到的问题

时间：2023.2.14，改bug一天，差点心态爆炸。情人节快乐 直接代入发现出差，找了以圈下来发现是特别shabi的错误，难受

2.1 程序bug

2.1.1 数据类型一定要注意

报错的页面

2.1.2 初始输入数据要和程序匹配

IGA_Plate程序处理的单元信息是从1开始的；而我在第一部分理论分析是从0开始的，不匹配

理论分析的ele矩阵

代入IGA_Plate程序中的ele单元应为：

报错的页面

2.1.3 有点傻的错误—总刚的行数和列数 num_GK 给错了

二个单元的ele.rows() = 2，所以之前没报错

应改为：

报错的页面

2.1.4 组总刚出错 — 向量下标超出范围

向量下标超出范围

之前case by case 求积分上下限有错
之前版本

改正后应为

废话不多说，直接上程序吧

2.2 主程序文件

找bug过程把好多IGA里面的函数又写在主程序文件，一行行的找错误（说多了都是泪）：

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