等几何程序保存

---

3.14下午三点

1 mian文件

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

// -----------------------------------------------------------------------------------------
// 基本初始物理量
double IGA_Plate::E = 1e6;            // 弹性模量
double IGA_Plate::mu = 0.3;           // 泊松比
double IGA_Plate::t = 0.01;           // 厚度
// -----------------------------------------------------------------------------------------

// ------------------------------------------------------------------------------------------------------------------
//细化前初始物理量
//节点坐标
MatrixXd Node_Matrix()
{
   
    //节点编号一行行排列，从左到右排列
    MatrixXd Node(12, 2);
    Node << 0, 0,
        0.5, 0,
        1.5, 0,
        2, 0,
        0, 1,
        0.5, 1,
        1.5, 1,
        2, 1,
        0, 2,
        0.5, 2,
        1.5, 2,
        2, 2;
    return Node;
}
//单元信息
MatrixXd ele_Matrix()
{
   
    //单元内包含的节点的全局编号，一行行排列
    MatrixXd ele(2, 9);
    ele << 1, 2, 3, 5, 6, 7, 9, 10, 11,
        2, 3, 4, 6, 7, 8, 10, 11, 12;
    return ele;
}
MatrixXd k_Refine::pre_Node = Node_Matrix();           //节点信息矩阵
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(MatrixXd K)
{
   
    for (int j = 0; j < K.rows(); j++)
    {
   
        for (int i = 0; i < K.cols(); i++)
        {
   
            cout << K(j, i) << "  ";
        }
        cout << endl;
    }
}
// 输出向量的各个值
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;   //  定义一个类对象
MatrixXd IGA_Plate::Node = span_2_3.Node;//节点信息矩阵
MatrixXd IGA_Plate::ele = span_2_3.ele;//单元信息       
vector<double> IGA_Plate::U = span_2_3.U;//节点向量U
vector<double> IGA_Plate::V = span_2_3.V;//节点向量V
MatrixXd IGA_Plate::BC = span_2_3.return_BC();//位移约束数据
// 细化次数
int k_Refine::k = 0;
// ------------------------------------------------------------------------------------------------------------------



int main()
{
   
    IGA_Plate iga_2d;
    Gauss gauss; //定义高斯类
    k_Refine k_refine;
    MatrixXd Node = iga_2d.Node;

    MatrixXd Node1 = k_refine.Node;

    cout << "\n 节点位移不变量 disp = \n" << Node.row(0) << endl;
    cout << "---------分割线1----------" << endl;
    cout << "\n Node = \n" << Node.row(0) << endl;
    Print(Node);

    vector<double> U1 = k_refine.U;
    cout << "\n U1 = " << endl;
    Print_Vector(U1);

    vector<double> U2 = k_refine.V;
    cout << "\n U2 = " << endl;
    Print_Vector(U2);

    MatrixXd ele1 = k_refine.ele;
    cout << "\n ele1 = \n" << ele1 << endl;

    MatrixXd BC(8, 2);
    BC.setZero();
    BC = k_refine.return_BC();
    cout << "\n BC = \n" << BC << endl;

    VectorXd disp1 = iga_2d.return_disp();   // 位移 displacement 
    cout << "\n disp = " << endl;
    Print(disp1);


    gsInfo << "here1" << endl;;
    cout << "---------分割线1----------" << endl;
    

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

}

2 Gauss.h文件

#pragma once
#include<iostream>
#include <Eigen/Dense>
using namespace Eigen;
using namespace std;
#include <vector>

class  Gauss
{
   
public:
	int n;        // 高斯点数
	MatrixXd gp;    // 标准高斯区间的高斯点
	MatrixXd wgt;   // 标准高斯区间的高斯权重



public:
	// 计算高斯映射的雅可比矩阵
	double jacobian_gauss_mapping(vector<double> elDoma)
	{
   
		// 计算高斯映射的雅可比矩阵
		// 输入
		//     elDoma : 传入的单元参数域     
		//         1维:   elDoma = [u1, u2] 
		//         2维:   elDoma = [u1, u2, v1, v2] 
		//         3维:   elDoma = [u1, u2, v1, v2, w1, w2] 
		// 输出: 
		//   jacobian --- 雅可比值
		// 
		// ----------------------------------------------------------------
		double jacobian;
		if (elDoma.size() == 2)   // 1维: 
		{
   
			jacobian = (elDoma[0, 1] - elDoma[0, 0]) / 2;
		}

		if (elDoma.size() == 4)   // 2维: 
		{
   
			double Jxi = (elDoma[0, 1] - elDoma[0, 0]) / 2;
			double Jeta = (elDoma[0, 3] - elDoma[0, 2]) / 2;

			jacobian = Jxi * Jeta;
		}

		if (elDoma.size() == 6)   // 3维: 
		{
   
			double Jxi = (elDoma[0, 1] - elDoma[0, 0]) / 2;
			double Jeta = (elDoma[0, 3] - elDoma[0, 2]) / 2;
			double Jzeta = (elDoma[0, 5] - elDoma[0, 4]) / 2;

			jacobian = Jxi * Jeta * Jzeta;
		}
		return jacobian;
	}


	//得到一维和二维高斯积分点及其权值   三维待写
	//先写两点高斯公式   必须先调用这个函数，才能得到类中 gp 和 wgt 的值
	void gauss_quadrature(vector<double> Gauss_num)
	{
   
		// ----------------------------------------------------------------
		//[gp, wgt] = gauss_quadrature(gp_x)
		// 输入
		//     elDoma : 传入的单元参数域   
		//         gp_x : x方向高斯点数     gp_y : y方向高斯点数...
		//         1维:   Gauss_num = {gp_x}     
		//         例子:  如调用两点高斯公式 vector<double> Gauss_num = { 2 } 
		//                如调用三点高斯公式 vector<double> Gauss_num = { 3 }
		//         2维:   Gauss_num = {gp_x,gp_y} 
		//         例子:  如调用U和V方向都用两点高斯公式 vector<double> Gauss_num = { 2 ,2 } 
		//         例子:  如调用两点高斯公式 vector<double> Gauss_num = { 3,3 } 
		//         3维:   Gauss_num = {gp_x,gp_y,gp_z} 
		// 输出: 
		//   gp   ------- 标准高斯区间的高斯点
		//   wgt  ------- 标准高斯区间的高斯权值
		//   需要调用这个类函数后  才能调用类中的这两个值
		// ----------------------------------------------------------------

		if (Gauss_num.size() == 1)