MatLib 1.0 alpha

#include "../../MatLib/MathLib.h" #include <iostream> using namespace std; using namespace MathLib; using namespace MathLib::Matrix; int solve() { //和求解问题相关的常数 const Int32 nPoint = 12; const Int32 nDOF = 2; const Int32 nTotalDOF = nPoint * nDOF; const Int32 nFixedPoint = 2; const Int32 nElement = 5; const Int32 nNode = 4; const Int32 nEDOF = nNode * nDOF; // const Int32 nDime = 2; //材料常数 const Float young = 1500; const Float poiss = 0.01; //高斯积分常数 const Int32 nGauss = 2; const FullMatrix GaussWeight(2, 1, 1); FullMatrix GaussPoint(2, 1); GaussPoint.assigner() = -0.577350269189626, 0.577350269189626; //局部坐标到全局坐标的映射 FullMatrix LocalNodes(nElement, nNode, 0); { for (int part = 1; part <= nElement; ++part) { Int32 last = (part - 1) * 2; LocalNodes(part, 1) = last + 1; LocalNodes(part, 2) = last + 3; LocalNodes(part, 3) = last + 4; LocalNodes(part, 4) = last + 2; } } //全局坐标 FullMatrix Coord(12, 2, 0); { const Int32 Height = 1; Coord.assigner() = 0, 0, 0, Height, 1, 0, 2, Height, 2, 0, 4, Height, 4, 0, 5, Height, 7, 0, 6, Height, 10, 0, 10, Height; } //本质边界条件 const FullMatrix FixedValue(nTotalDOF, 1, 0); BasicMatrix<Int32> FixedPoint(1, nFixedPoint, 0); BasicMatrix<Int32> BoundDOF(nFixedPoint, nDOF, 0); BasicMatrix<Int32> IsFixed(1, nTotalDOF); { FixedPoint(1, 1) = 1; BoundDOF(1, 1) = 0; BoundDOF(1, 2) = 0; IsFixed(1, 1) = 1; IsFixed(1, 2) = 1; FixedPoint(1, 2) = 2; BoundDOF(2, 1) = 0; BoundDOF(2, 2) = 1; IsFixed(1, 3) = 1; IsFixed(1, 4) = 0; } //D矩阵 const Float E0 = young / (1 - poiss * poiss); const Float mu0 = poiss / (1 - poiss); FullMatrix DMatrix(3, 3, 0); { DMatrix.assigner() = 1, mu0, 0, mu0, 1, 0, 0, 0, 0.5 * (1 - mu0); DMatrix *= E0 / (1 - mu0 * mu0); } //总刚矩阵和总体载荷 FullMatrix GStiff(nTotalDOF, nTotalDOF, 0); FullMatrix GLoad(nTotalDOF, 1, 0); for (Int32 CurrElem = 1; CurrElem <= nElement; ++CurrElem) { //单元形状参数 FullMatrix XNode(nNode); FullMatrix YNode(nNode); for (Int32 i = 1; i <= nNode; ++i) { Int32 idx = (Int32)LocalNodes(CurrElem, i); XNode(i) = Coord(idx, 1); YNode(i) = Coord(idx, 2); } Float a1 = (-XNode(1) + XNode(2) + XNode(3) - XNode(4)) / 4; Float a2 = ( XNode(1) - XNode(2) + XNode(3) - XNode(4)) / 4; Float a3 = (-XNode(1) - XNode(2) + XNode(3) + XNode(4)) / 4; Float b1 = (-YNode(1) + YNode(2) + YNode(3) - YNode(4)) / 4; Float b2 = ( YNode(1) - YNode(2) + YNode(3) - YNode(4)) / 4; Float b3 = (-YNode(1) - YNode(2) + YNode(3) + YNode(4)) / 4; Float J0 = a1 * b3 - a3 * b1; Float J1 = a1 * b2 - a2 * b1; Float J2 = a2 * b3 - a3 * b2; //单刚矩阵 FullMatrix EStiff(nEDOF, nEDOF, 0); for (Int32 iGauss = 1; iGauss <= nGauss; ++iGauss) { Float x = GaussPoint(iGauss); for (Int32 jGauss = 1; jGauss <= nGauss; ++jGauss) { Float y = GaussPoint(jGauss); FullMatrix Deriv4(2, 4); { FullMatrix Right(4, 4); FullMatrix Left(1, 4); Left.assigner() = 0, 1, 0, y; Right.assigner() = 1, 1, 1, 1, -1, 1, 1,-1, -1,-1, 1, 1, 1,-1, 1,-1; Deriv4(1, 1, 1, 4) = 0.25 * Left * Right; Left.assigner() = 0, 0, 1, x; Deriv4(2, 2, 1, 4) = 0.25 * Left * Right; } Float djacb = J0 + J1 * x + J2 * y; Float dvolu = djacb * GaussWeight(iGauss) * GaussWeight(jGauss); FullMatrix xjaci(2, 2); { xjaci.assigner() = b3 + b2 * x, -(b1 + b2 * y), -(a3 + a2 * x), a1 + a2 * y; xjaci /= djacb; } FullMatrix bmatrix(3, 8); { FullMatrix cartd(2, 4); cartd = xjaci * Deriv4; bmatrix.assigner() = cartd(1, 1), 0, cartd(1, 2), 0, cartd(1, 3), 0, cartd(1, 4), 0, 0, cartd(2, 1), 0, cartd(2, 2), 0, cartd(2, 3), 0, cartd(2, 4), cartd(2, 1), cartd(1, 1), cartd(2, 2), cartd(1, 2), cartd(2, 3), cartd(1, 3), cartd(2, 4), cartd(1, 4); } EStiff += dvolu * Tran(bmatrix) * DMatrix * bmatrix; } } //单元载荷 FullMatrix EForce(nEDOF, 1, 0); { if(CurrElem == nElement) { EForce.assigner() = 0, 0, 1000, 0, -1000, 0, 0, 0; } } //单刚矩阵到总刚矩阵的组装 BasicMatrix<Int32> index(8); for (Int32 inode = 1; inode <= nNode; ++inode) { Int32 idx = (Int32)LocalNodes(CurrElem, inode); for (Int32 idof = 1; idof <= nDOF; ++idof) { index(nDOF * (inode - 1) + idof) = nDOF * (idx - 1) + idof; } } for (Int32 i = 1; i <= nEDOF; ++i) { Int32 idxi = index(i); GLoad(idxi) = GLoad(idxi) + EForce(i); for (Int32 j = 1; j <= nEDOF; ++j) { Int32 idxj = index(j); GStiff(idxi, idxj) += EStiff(i, j); } } } //使用本质边界条件 for (Int32 i = 1; i <= nFixedPoint; ++i) for (Int32 j = 1; j <= nDOF; ++j) { Int32 idx = (FixedPoint(1, i) - 1) * nDOF + j; if (IsFixed(1, idx)) { Float value = FixedValue(idx); for (Int32 itotv = 1; itotv <= nTotalDOF; ++itotv) { GLoad(itotv) -= value * GStiff(itotv, idx); GStiff(itotv, idx) = 0; GStiff(idx, itotv) = 0; } GStiff(idx, idx) = 1; GLoad(idx, 1) = value; } } //求解 FullMatrix res(GStiff.row(), 1); GaussSolve(res, GStiff, GLoad); cout.precision(16); cout<< res; cout << res(21) << '/t' << res(22) << endl; return 0; }

 

 给我发邮件可以获得本库的代码

mail:bailiang@vip.qq.com