MatLib 1.0 alpha

#include"../../MatLib/MathLib.h"#include<iostream>using namespace std;using namespace MathLib;using namespace MathLib::Matrix;intsolve(){//和求解问题相关的常数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 FullMatrixGaussWeight(2, 1, 1);FullMatrixGaussPoint(2, 1);GaussPoint.assigner()=-0.577350269189626, 0.577350269189626;//局部坐标到全局坐标的映射FullMatrixLocalNodes(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;}}//全局坐标FullMatrixCoord(12, 2, 0);{const Int32Height=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 FullMatrixFixedValue(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 FloatE0=young / (1 - poiss * poiss);const Floatmu0=poiss / (1 - poiss);FullMatrixDMatrix(3, 3, 0);{DMatrix.assigner()=1, mu0, 0, mu0, 1, 0, 0, 0, 0.5 * (1 - mu0);DMatrix*=E0 / (1 - mu0 * mu0);}//总刚矩阵和总体载荷FullMatrixGStiff(nTotalDOF, nTotalDOF, 0);FullMatrixGLoad(nTotalDOF, 1, 0);for (Int32 CurrElem = 1; CurrElem <= nElement; ++CurrElem){//单元形状参数FullMatrixXNode(nNode);FullMatrixYNode(nNode);for (Int32 i = 1; i <= nNode; ++i){Int32 idx=(Int32)LocalNodes(CurrElem, i);XNode(i)=Coord(idx, 1);YNode(i)=Coord(idx, 2);}Floata1=(-XNode(1) + XNode(2) + XNode(3) - XNode(4)) / 4;Floata2=( XNode(1) - XNode(2) + XNode(3) - XNode(4)) / 4;Floata3=(-XNode(1) - XNode(2) + XNode(3) + XNode(4)) / 4;Floatb1=(-YNode(1) + YNode(2) + YNode(3) - YNode(4)) / 4;Floatb2=( YNode(1) - YNode(2) + YNode(3) - YNode(4)) / 4;Floatb3=(-YNode(1) - YNode(2) + YNode(3) + YNode(4)) / 4;FloatJ0=a1 * b3 - a3 * b1;FloatJ1=a1 * b2 - a2 * b1;FloatJ2=a2 * b3 - a3 * b2;//单刚矩阵FullMatrixEStiff(nEDOF, nEDOF, 0);for (Int32 iGauss = 1; iGauss <= nGauss; ++iGauss){Floatx=GaussPoint(iGauss);for (Int32 jGauss = 1; jGauss <= nGauss; ++jGauss){Float y=GaussPoint(jGauss);FullMatrixDeriv4(2, 4);{FullMatrixRight(4, 4);FullMatrixLeft(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;}Floatdjacb = J0 + J1 * x + J2 * y;Floatdvolu = djacb * GaussWeight(iGauss) * GaussWeight(jGauss);FullMatrixxjaci(2, 2);{xjaci.assigner() = b3 + b2 * x, -(b1 + b2 * y), -(a3 + a2 * x), a1 + a2 * y;xjaci /= djacb;}FullMatrixbmatrix(3, 8);{FullMatrixcartd(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; }}//单元载荷FullMatrixEForce(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)){Floatvalue = 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;}}//求解FullMatrixres(GStiff.row(), 1); GaussSolve(res, GStiff, GLoad);cout.precision(16);cout<< res;cout << res(21) << '/t' << res(22) << endl;return 0;}

 

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

mail:[e-mail]