C++实现有限元计算 矩阵装配Assembly类

本系列文章致力于实现“手搓有限元，干翻Ansys的目标”，基本框架为前端显示使用QT实现交互，后端计算采用Visual Studio C++。

1、Assembly

Assembly类用来实现矩阵的装配，将单元的全局刚度矩阵、位移列矩阵、外载列矩阵装配组装为整体全局刚度矩阵、位移列矩阵和外载列矩阵，作为后续求解位移列阵的基础。

Assembly类组成架构

1.2、public function

1.2.1、构造函数与析构函数

矩阵装配过程中不需要进行初始化，所以构造函数与析构函数只做定义，不做实现（空实现）。

Assembly.h函数声明：

//***********************构造函数析构函数***********************//
/*
函数名称:		无参构造函数
*/
Assembly();/*
函数名称:		析构函数
*/
~Assembly();

Assembly.cpp函数实现：（空实现）

1.2.2、单元全局刚度矩阵组装

将单元的全局刚度矩阵装配为整体全局刚度矩阵，单元包括上一期的二维杆单元Bar2D2Node（点击此处跳转Bar2D2Node介绍）、以及将来会介绍的三维杆单元Bar3D2Node、二维截面为矩形的梁单元Rect_Beam2D2Node等等。

Assembly.h函数声明：

//***************************组装矩阵***************************//
/*
刚度矩阵的组装Bar2D2Node
bar:				二维杆单元数组
numPoint:			节点数量
numElement:			单元数量
*/
Matrix AssemblyMatK(std::vector<Bar2D2Node>bar, int numPoint, int numElement);/*
刚度矩阵的组装Bar3D2Node
bar:				三维杆单元数组
numPoint:			节点数量
numElement:			单元数量
*/
Matrix AssemblyMatK(std::vector<Bar3D2Node>bar, int numPoint, int numElement);/*
刚度矩阵的组装Rect_Beam2D2Node
beam:				二维矩形截面梁单元数组
numPoint:			节点数量
numElement:			单元数量
*/
Matrix AssemblyMatK(std::vector<Rect_Beam2D2Node>beam, int numPoint, int numElement);

Assembly.cpp函数实现：

//刚度矩阵的组装Bar2D2Node
Matrix Assembly::AssemblyMatK(std::vector<Bar2D2Node>bar, int numPoint, int numElement)
{//建立组装刚度矩阵 由于是2D问题 矩阵维度为2*numPointthis->m_SumK.CopyMat(Matrix(2 * numPoint, 2 * numPoint));for (int iEle = 0; iEle < numElement; iEle++){//获取点IDint pointID0 = bar[iEle].GetPoint0()->GetID();int pointID1 = bar[iEle].GetPoint1()->GetID();//ID List [2i, 2i+1, 2j, 2j+1]int pointIndex[4] = { 2 * pointID0, 2 * pointID0 + 1, 2* pointID1, 2* pointID1 + 1};//2D问题 逐个组装刚度矩阵for (int iRow = 0; iRow < 4; iRow++){for (int iCol = 0; iCol < 4; iCol++){double Kele = bar[iEle].GetK().GetMatrixEle(iRow, iCol) + this->m_SumK.GetMatrixEle(pointIndex[iRow], pointIndex[iCol]);this->m_SumK.SetMatrixEle(pointIndex[iRow], pointIndex[iCol], Kele);}}}return this->m_SumK;
}//刚度矩阵的组装Bar3D2Node
Matrix Assembly::AssemblyMatK(std::vector<Bar3D2Node>bar, int numPoint, int numElement)
{//建立组装刚度矩阵 由于是3D问题 矩阵维度为3*numPointthis->m_SumK.CopyMat(Matrix(3 * numPoint, 3 * numPoint));for (int iEle = 0; iEle < numElement; iEle++){//获取点IDint pointID0 = bar[iEle].GetPoint0()->GetID();int pointID1 = bar[iEle].GetPoint1()->GetID();//ID List [3i, 3i+1, 3i+2, 3j, 3j+1, 3j+2]int pointIndex[6] = { 3 * pointID0 + 0, 3 * pointID0 + 1, 3 * pointID0 + 2,3 * pointID1 + 0, 3 * pointID1 + 1, 3 * pointID1 + 2};//3D问题 逐个组装刚度矩阵for (int iRow = 0; iRow < 6; iRow++){for (int iCol = 0; iCol < 6; iCol++){//刚度矩阵累加double Kele = bar[iEle].GetK().GetMatrixEle(iRow, iCol) + this->m_SumK.GetMatrixEle(pointIndex[iRow], pointIndex[iCol]);//矩阵元素装填this->m_SumK.SetMatrixEle(pointIndex[iRow], pointIndex[iCol], Kele);}}}return this->m_SumK;
}//刚度矩阵的组装Rect_Beam2D2Node
Matrix Assembly::AssemblyMatK(std::vector<Rect_Beam2D2Node>beam, int numPoint, int numElement)
{//建立组装刚度矩阵 由于是2D问题 矩形梁每个节点自由度为3 矩阵维度为3*numPointthis->m_SumK.CopyMat(Matrix(3 * numPoint, 3 * numPoint));for (int iEle = 0; iEle < numElement; iEle++){//获取点IDint pointID0 = beam[iEle].GetPoint0()->GetID();int pointID1 = beam[iEle].GetPoint1()->GetID();//ID List [3i, 3i+1, 3i+2, 3j, 3j+1, 3j+2]int pointIndex[6] = { 3 * pointID0 + 0, 3 * pointID0 + 1, 3 * pointID0 + 2,3 * pointID1 + 0, 3 * pointID1 + 1, 3 * pointID1 + 2 };//3D问题 逐个组装刚度矩阵for (int iRow = 0; iRow < 6; iRow++){for (int iCol = 0; iCol < 6; iCol++){//刚度矩阵累加double Kele = beam[iEle].GetK().GetMatrixEle(iRow, iCol) + this->m_SumK.GetMatrixEle(pointIndex[iRow], pointIndex[iCol]);//矩阵元素装填this->m_SumK.SetMatrixEle(pointIndex[iRow], pointIndex[iCol], Kele);}}}return this->m_SumK;
}

1.2.3、位移列矩阵组装

将各个单元的节点位移装配为整体节点位移。注意这里只区分为二维节点Point2D（点击此处可跳转到Point2D的介绍）和三维节点Point3D。因为对于不同单元组成的节点无非就二维节点与三维节点两种情况。

Assembly.h函数声明：

/*
线位移、角位移-矩阵组装Point2D
p:					二维节点数组
numPoint:			节点数量
ElementStyle:		单元类型枚举
*/
Matrix AssemblyMatU(std::vector<Point2D> p, int numPoint, ElementStyle eleStyle);/*
线位移、角位移-矩阵组装Point3D
p:					三维节点数组
numPoint:			节点数量
ElementStyle:		单元类型枚举
*/
Matrix AssemblyMatU(std::vector<Point3D> p, int numPoint, ElementStyle eleStyle);

Assembly.cpp函数实现：

//线位移、角位移 - 矩阵组装Point2D
Matrix Assembly::AssemblyMatU(std::vector<Point2D> p, int numPoint, ElementStyle eleStyle)
{//创建结果矩阵Matrix* resMat;switch (eleStyle){case BAR_2D2NODE://位移整体矩阵 2D 两个自由度(杆单元 存在两个自由度)resMat = new Matrix(2 * numPoint, 1);//遍历节点for (int i = 0; i < numPoint; i++){resMat->SetMatrixEle(2 * i + 0, 0, p[i].GetU());			//UresMat->SetMatrixEle(2 * i + 1, 0, p[i].GetV());			//V}break;case RECT_BEAM_2D2NODE://位移整体矩阵 2D 三个自由度(梁单元 存在弯曲变形)resMat = new Matrix(3 * numPoint, 1);//遍历节点for (int i = 0; i < numPoint; i++){resMat->SetMatrixEle(3 * i + 0, 0, p[i].GetU());			//UresMat->SetMatrixEle(3 * i + 1, 0, p[i].GetV());			//VresMat->SetMatrixEle(3 * i + 2, 0, p[i].GetTheatZ());		//TheatZ}break;default:resMat = new Matrix(3 * numPoint, 1);PRINTF_LOCAL;break;}return *resMat;
}//线位移、角位移-矩阵组装Point3D
Matrix Assembly::AssemblyMatU(std::vector<Point3D> p, int numPoint, ElementStyle eleStyle)
{//创建结果矩阵Matrix* resMat;switch (eleStyle){case BAR_3D2NODE://位移整体矩阵 3DresMat = new Matrix(3 * numPoint, 1);//遍历节点 逐个填取位移数据for (int i = 0; i < numPoint; i++){resMat->SetMatrixEle(3 * i + 0, 0, p[i].GetU());			//UresMat->SetMatrixEle(3 * i + 1, 0, p[i].GetV());			//VresMat->SetMatrixEle(3 * i + 2, 0, p[i].GetW());			//W}break;default:resMat = new Matrix(3 * numPoint, 1);PRINTF_LOCAL;break;}return *resMat;
}

1.2.4、外载列矩阵组装

将各个单元的节点外载装配为整体节点外载。

Assembly.h函数声明：

/*
力、力矩-矩阵组装Point2D
p:					二维节点数组
numPoint:			节点数量
ElementStyle:		单元类型枚举
*/
Matrix AssemblyMatF(std::vector<Point2D> p, int numPoint, ElementStyle eleStyle);/*
力、力矩-矩阵组装Point3D
p:					三维节点数组
numPoint:			节点数量
ElementStyle:		单元类型枚举
*/
Matrix AssemblyMatF(std::vector<Point3D> p, int numPoint, ElementStyle eleStyle);

Assembly.cpp函数实现：

//力、力矩-矩阵组装Point2D
Matrix Assembly::AssemblyMatF(std::vector<Point2D> p, int numPoint, ElementStyle eleStyle)
{Matrix* resMat;switch (eleStyle){case BAR_2D2NODE://结果矩阵 2D 两个自由度resMat = new Matrix(2 * numPoint, 1);//遍历节点for (int i = 0; i < numPoint; i++){resMat->SetMatrixEle(2 * i + 0, 0, p[i].GetFx());			//UresMat->SetMatrixEle(2 * i + 1, 0, p[i].GetFy());			//V}break;case RECT_BEAM_2D2NODE://结果矩阵 2D 三个自由度resMat = new Matrix(3 * numPoint, 1);//遍历节点for (int i = 0; i < numPoint; i++){resMat->SetMatrixEle(3 * i + 0, 0, p[i].GetFx());			//UresMat->SetMatrixEle(3 * i + 1, 0, p[i].GetFy());			//VresMat->SetMatrixEle(3 * i + 2, 0, p[i].GetMz());			//Mz}break;default:resMat = new Matrix(3 * numPoint, 1);PRINTF_LOCAL;break;}return *resMat;
}//力、力矩-矩阵组装Point3D
Matrix Assembly::AssemblyMatF(std::vector<Point3D> p, int numPoint, ElementStyle eleStyle)
{//创建结果矩阵Matrix* resMat;switch (eleStyle){case BAR_3D2NODE://结果矩阵 3DresMat = new Matrix(3 * numPoint, 1);//遍历节点for (int i = 0; i < numPoint; i++){resMat->SetMatrixEle(3 * i + 0, 0, p[i].GetFx());			//x方向外载resMat->SetMatrixEle(3 * i + 1, 0, p[i].GetFy());			//y方向外载resMat->SetMatrixEle(3 * i + 2, 0, p[i].GetFz());			//z方向外载}break;default:resMat = new Matrix(3 * numPoint, 1);break;}return *resMat;
}

1.3、private variable

私有变量，只能在类内进行调用。

private:Matrix m_SumK;			//组装后的刚度矩阵

1.4、全部源码

Assembly.h函数声明：

#ifndef _ASSEMBLY_H_
#define _ASSEMBLY_H_#include "Bar2D2Node.h"
#include "Bar3D2Node.h"
#include "Beam2D2Node.h"
#include "Rectangle_Beam2D2Node.h"
#include "GlobalVars.h"
#include <vector>class Assembly
{
public://***********************构造函数析构函数***********************///*函数名称:		无参构造函数*/Assembly();/*函数名称:		析构函数*/~Assembly();//***************************组装矩阵***************************///*刚度矩阵的组装Bar2D2Nodebar:				二维杆单元数组numPoint:			节点数量numElement:			单元数量*/Matrix AssemblyMatK(std::vector<Bar2D2Node>bar, int numPoint, int numElement);/*刚度矩阵的组装Bar3D2Nodebar:				三维杆单元数组numPoint:			节点数量numElement:			单元数量*/Matrix AssemblyMatK(std::vector<Bar3D2Node>bar, int numPoint, int numElement);/*刚度矩阵的组装Rect_Beam2D2Nodebeam:				二维矩形截面梁单元数组numPoint:			节点数量numElement:			单元数量*/Matrix AssemblyMatK(std::vector<Rect_Beam2D2Node>beam, int numPoint, int numElement);/*线位移、角位移-矩阵组装Point2Dp:					二维节点数组numPoint:			节点数量ElementStyle:		单元类型枚举*/Matrix AssemblyMatU(std::vector<Point2D> p, int numPoint, ElementStyle eleStyle);/*线位移、角位移-矩阵组装Point3Dp:					三维节点数组numPoint:			节点数量ElementStyle:		单元类型枚举*/Matrix AssemblyMatU(std::vector<Point3D> p, int numPoint, ElementStyle eleStyle);/*力、力矩-矩阵组装Point2Dp:					二维节点数组numPoint:			节点数量ElementStyle:		单元类型枚举*/Matrix AssemblyMatF(std::vector<Point2D> p, int numPoint, ElementStyle eleStyle);/*力、力矩-矩阵组装Point3Dp:					三维节点数组numPoint:			节点数量ElementStyle:		单元类型枚举*/Matrix AssemblyMatF(std::vector<Point3D> p, int numPoint, ElementStyle eleStyle);private:Matrix m_SumK;			//组装后的刚度矩阵
};#endif

Assembly.cpp函数实现：

#include "Assembly.h"Assembly::Assembly()
{
}Assembly::~Assembly()
{
}//刚度矩阵的组装Bar2D2Node
Matrix Assembly::AssemblyMatK(std::vector<Bar2D2Node>bar, int numPoint, int numElement)
{//建立组装刚度矩阵 由于是2D问题 矩阵维度为2*numPointthis->m_SumK.CopyMat(Matrix(2 * numPoint, 2 * numPoint));for (int iEle = 0; iEle < numElement; iEle++){//获取点IDint pointID0 = bar[iEle].GetPoint0()->GetID();int pointID1 = bar[iEle].GetPoint1()->GetID();//ID List [2i, 2i+1, 2j, 2j+1]int pointIndex[4] = { 2 * pointID0, 2 * pointID0 + 1, 2* pointID1, 2* pointID1 + 1};//2D问题 逐个组装刚度矩阵for (int iRow = 0; iRow < 4; iRow++){for (int iCol = 0; iCol < 4; iCol++){double Kele = bar[iEle].GetK().GetMatrixEle(iRow, iCol) + this->m_SumK.GetMatrixEle(pointIndex[iRow], pointIndex[iCol]);this->m_SumK.SetMatrixEle(pointIndex[iRow], pointIndex[iCol], Kele);}}}return this->m_SumK;
}//刚度矩阵的组装Bar3D2Node
Matrix Assembly::AssemblyMatK(std::vector<Bar3D2Node>bar, int numPoint, int numElement)
{//建立组装刚度矩阵 由于是3D问题 矩阵维度为3*numPointthis->m_SumK.CopyMat(Matrix(3 * numPoint, 3 * numPoint));for (int iEle = 0; iEle < numElement; iEle++){//获取点IDint pointID0 = bar[iEle].GetPoint0()->GetID();int pointID1 = bar[iEle].GetPoint1()->GetID();//ID List [3i, 3i+1, 3i+2, 3j, 3j+1, 3j+2]int pointIndex[6] = { 3 * pointID0 + 0, 3 * pointID0 + 1, 3 * pointID0 + 2,3 * pointID1 + 0, 3 * pointID1 + 1, 3 * pointID1 + 2};//3D问题 逐个组装刚度矩阵for (int iRow = 0; iRow < 6; iRow++){for (int iCol = 0; iCol < 6; iCol++){//刚度矩阵累加double Kele = bar[iEle].GetK().GetMatrixEle(iRow, iCol) + this->m_SumK.GetMatrixEle(pointIndex[iRow], pointIndex[iCol]);//矩阵元素装填this->m_SumK.SetMatrixEle(pointIndex[iRow], pointIndex[iCol], Kele);}}}return this->m_SumK;
}//刚度矩阵的组装Rect_Beam2D2Node
Matrix Assembly::AssemblyMatK(std::vector<Rect_Beam2D2Node>beam, int numPoint, int numElement)
{//建立组装刚度矩阵 由于是2D问题 矩形梁每个节点自由度为3 矩阵维度为3*numPointthis->m_SumK.CopyMat(Matrix(3 * numPoint, 3 * numPoint));for (int iEle = 0; iEle < numElement; iEle++){//获取点IDint pointID0 = beam[iEle].GetPoint0()->GetID();int pointID1 = beam[iEle].GetPoint1()->GetID();//ID List [3i, 3i+1, 3i+2, 3j, 3j+1, 3j+2]int pointIndex[6] = { 3 * pointID0 + 0, 3 * pointID0 + 1, 3 * pointID0 + 2,3 * pointID1 + 0, 3 * pointID1 + 1, 3 * pointID1 + 2 };//3D问题 逐个组装刚度矩阵for (int iRow = 0; iRow < 6; iRow++){for (int iCol = 0; iCol < 6; iCol++){//刚度矩阵累加double Kele = beam[iEle].GetK().GetMatrixEle(iRow, iCol) + this->m_SumK.GetMatrixEle(pointIndex[iRow], pointIndex[iCol]);//矩阵元素装填this->m_SumK.SetMatrixEle(pointIndex[iRow], pointIndex[iCol], Kele);}}}return this->m_SumK;
}//线位移、角位移 - 矩阵组装Point2D
Matrix Assembly::AssemblyMatU(std::vector<Point2D> p, int numPoint, ElementStyle eleStyle)
{//创建结果矩阵Matrix* resMat;switch (eleStyle){case BAR_2D2NODE://位移整体矩阵 2D 两个自由度(杆单元 存在两个自由度)resMat = new Matrix(2 * numPoint, 1);//遍历节点for (int i = 0; i < numPoint; i++){resMat->SetMatrixEle(2 * i + 0, 0, p[i].GetU());			//UresMat->SetMatrixEle(2 * i + 1, 0, p[i].GetV());			//V}break;case RECT_BEAM_2D2NODE://位移整体矩阵 2D 三个自由度(梁单元 存在弯曲变形)resMat = new Matrix(3 * numPoint, 1);//遍历节点for (int i = 0; i < numPoint; i++){resMat->SetMatrixEle(3 * i + 0, 0, p[i].GetU());			//UresMat->SetMatrixEle(3 * i + 1, 0, p[i].GetV());			//VresMat->SetMatrixEle(3 * i + 2, 0, p[i].GetTheatZ());		//TheatZ}break;default:resMat = new Matrix(3 * numPoint, 1);PRINTF_LOCAL;break;}return *resMat;
}//线位移、角位移-矩阵组装Point3D
Matrix Assembly::AssemblyMatU(std::vector<Point3D> p, int numPoint, ElementStyle eleStyle)
{//创建结果矩阵Matrix* resMat;switch (eleStyle){case BAR_3D2NODE://位移整体矩阵 3DresMat = new Matrix(3 * numPoint, 1);//遍历节点 逐个填取位移数据for (int i = 0; i < numPoint; i++){resMat->SetMatrixEle(3 * i + 0, 0, p[i].GetU());			//UresMat->SetMatrixEle(3 * i + 1, 0, p[i].GetV());			//VresMat->SetMatrixEle(3 * i + 2, 0, p[i].GetW());			//W}break;default:resMat = new Matrix(3 * numPoint, 1);PRINTF_LOCAL;break;}return *resMat;
}//力、力矩-矩阵组装Point2D
Matrix Assembly::AssemblyMatF(std::vector<Point2D> p, int numPoint, ElementStyle eleStyle)
{Matrix* resMat;switch (eleStyle){case BAR_2D2NODE://结果矩阵 2D 两个自由度resMat = new Matrix(2 * numPoint, 1);//遍历节点for (int i = 0; i < numPoint; i++){resMat->SetMatrixEle(2 * i + 0, 0, p[i].GetFx());			//UresMat->SetMatrixEle(2 * i + 1, 0, p[i].GetFy());			//V}break;case RECT_BEAM_2D2NODE://结果矩阵 2D 三个自由度resMat = new Matrix(3 * numPoint, 1);//遍历节点for (int i = 0; i < numPoint; i++){resMat->SetMatrixEle(3 * i + 0, 0, p[i].GetFx());			//UresMat->SetMatrixEle(3 * i + 1, 0, p[i].GetFy());			//VresMat->SetMatrixEle(3 * i + 2, 0, p[i].GetMz());			//Mz}break;default:resMat = new Matrix(3 * numPoint, 1);PRINTF_LOCAL;break;}return *resMat;
}//力、力矩-矩阵组装Point3D
Matrix Assembly::AssemblyMatF(std::vector<Point3D> p, int numPoint, ElementStyle eleStyle)
{//创建结果矩阵Matrix* resMat;switch (eleStyle){case BAR_3D2NODE://结果矩阵 3DresMat = new Matrix(3 * numPoint, 1);//遍历节点for (int i = 0; i < numPoint; i++){resMat->SetMatrixEle(3 * i + 0, 0, p[i].GetFx());			//x方向外载resMat->SetMatrixEle(3 * i + 1, 0, p[i].GetFy());			//y方向外载resMat->SetMatrixEle(3 * i + 2, 0, p[i].GetFz());			//z方向外载}break;default:resMat = new Matrix(3 * numPoint, 1);break;}return *resMat;
}