《OpenSees開發(fā)(二)源碼分析》由會員分享，可在線閱讀，更多相關(guān)《OpenSees開發(fā)(二)源碼分析（15頁珍藏版）》請在裝配圖網(wǎng)上搜索。

1、OpenSees 開發(fā)（二）源碼分析 這是一個平面桁架靜力分析算例，代碼位于 OpenSees2.3.0\EXAMPLES\Example1\main.cpp 這里先給 出原始源代碼： [cpp] view plain copy // standard C++ includes #include <stdlib.h> #include <OPS_Globals.h> #include <StandardStream.h> #include <ArrayOfTaggedObjects.h> // includes for the doma

2、in classes #include <Domain.h> #include <Node.h> #include <Truss.h> #include <ElasticMaterial.h> #include <SP_Constraint.h> #include <LoadPattern.h> #include <LinearSeries.h> #include <NodalLoad.h> // includes for the analysis classes #include <St

3、aticAnalysis.h> #include <AnalysisModel.h> #include <Linear.h> #include <PenaltyConstraintHandler.h> #include <DOF_Numberer.h> #include <RCM.h> #include <LoadControl.h> #include <BandSPDLinSOE.h> #include <BandSPDLinLapackSolver.h> // init the

4、global variabled defined in OPS_Globals.h StandardStream sserr; OPS_Stream *opserrPtr = &sserr; double ops_Dt = 0; // Domain *ops_TheActiveDomain = 0; Element *ops_TheActiveElement = 0; // main routine int main(int argc, char **argv) { // // now create a domain and a modelbuilder // and bu

5、ild the model // Domain *theDomain = new Domain(); // create the nodes using constructor: // Node(tag, ndof, crd1, crd2) // and then add them to the domain Node *node1 = new Node(1, 2, 0.0, 0.0); Node *node2 = new Node(2, 2, 144.0, 0.0); Node *node3 = new Node(3, 2, 168.0, 0.0); new Nod

6、e(4, 2, 72.0, 96.0); theDomain->addNode(node1); theDomain->addNode(node2); theDomain->addNode(node3); theDomain->addNode(node4); elastic material using constriuctor: Node *node4 = // create an // ElasticMaterialModel(tag, E) UniaxialMaterial *theMaterial = new ElasticMaterial(

7、1, 3000); //  create the truss elements using constructor: // Truss(tag, dim, nd1, nd2, Material &,A) // and then add them to the domain Truss *truss1 = new Truss(1, 2, 1, 4, *theMaterial, 10.0); Truss *truss2 = new Truss(2, 2, 2, 4, *theMaterial, 5.0); Truss *truss3 = new Truss(3, 2, 3

8、, 4, *theMaterial, 5.0); theDomain->addElement(truss1); theDomain->addElement(truss2); theDomain->addElement(truss3); // create the single-point constraint objects using constructor: // // SP_Constraint(tag, nodeTag, dofID, value) and then add them to the domain SP_Constraint *sp1 =

9、 new SP_Constraint(1, 1, 0, 0.0); SP_Constraint *sp2 = new SP_Constraint(2, 1, 1, 0.0); SP_Constraint *sp3 = new SP_Constraint(3, 2, 0, 0.0); SP_Constraint *sp4 = new SP_Constraint(4, 2, 1, 0.0); SP_Constraint *sp5 = new SP_Constraint(5, 3, 0, 0.0); SP_Constraint *sp6 = new SP_Constraint(6, 3,

10、1, 0.0); theDomain->addSP_Constraint(sp1); theDomain->addSP_Constraint(sp2); theDomain->addSP_Constraint(sp3); theDomain->addSP_Constraint(sp4); // theDomain->addSP_Constraint(sp5); theDomain->addSP_Constraint(sp6); construct a linear time series object using constructor:

11、 // LinearSeries() TimeSeries *theSeries = new LinearSeries(); // construct a load pattren using constructor: // LoadPattern(tag) // and then set its TimeSeries and add it to the domain LoadPattern *theLoadPattern = new LoadPattern(1); theLoadPattern->setTimeSeries(theSeries); theDomain->a

12、ddLoadPattern(theLoadPattern); // construct a nodal load using constructor: // NodalLoad(tag, nodeID, Vector &) // first construct a Vector of size 2 and set the values NOTE C INDEXING // then construct the load and add it to the domain Vector theLoadValues(2); theLoadValues(0) = 100.0; theLo

13、adValues(1) = -50.0; NodalLoad *theLoad = new NodalLoad(1, 4, theLoadValues); theDomain->addNodalLoad(theLoad, 1); // create an Analysis object to perform a static analysis of the model // - constructs: // AnalysisModel of type AnalysisModel, // EquiSolnAlgo of type  Linear // StaticIn

14、tegrator of type LoadControl // ConstraintHandler of type Penalty // DOFNumberer which uses RCM // LinearSOE of type Band SPD object AnalysisModel(); = new Linear(); // and then the StaticAnalysis AnalysisModel EquiSolnAlgo StaticIntegrator new LoadControl(1.0, 1, 1.0, 1.0); *theModel =

15、 new *theSolnAlgo *theIntegrator = ConstraintHandler RCM DOF_Numberer *theHandler = new PenaltyConstraintHandler(1.0e8,1.0e8); *theRCM = new RCM(); *theNumberer = new DOF_Numberer(*theRCM); BandSPDLinSolver *theSolver = new BandSPDLinLapackSolver(); LinearSOE *theSOE = new BandSPDL

16、inSOE(*theSolver); StaticAnalysis theAnalysis(*theDomain, *theHandler, *theNumberer, *theModel, *theSolnAlgo, *theSOE, *theIntegrator); // perform the analysis & print out the results for the domain int numSteps = 1; theAnalysis.analyze(numSteps); opserr << *theDomain; exit(0); }

17、接下去一步一步解釋 代碼： [cpp] view plain copy // 創(chuàng)建一個有限元模型 Domain *theDomain = new Domain(); [cpp] view plain copy // 創(chuàng)建 4 個節(jié)點， 詳細見說明 1 Node *node1 = new Node(1, 2, 0.0, 0.0); Node *node2 = new Node(2, 2, 144.0, 0.0); Node *node3 = new Node(3, 2, 168.0, 0.0); Node *node4 = new Node(4, 2, 72.0, 96.0); 說

18、明 1 ： Node 構(gòu)造函數(shù) 位于 OpenSees2.3.0\SRC\domain\node\Node.cpp 源碼定義如下： ***************************************************** Node::Node(int tag, int ndof, double Crd1, double Crd2) :DomainComponent(tag,NOD_TAG_Node), numberDOF(ndof), theDOF_GroupPtr(0), Crd(0), 。。。。。。 Crd = new Vector(2); (*C

19、rd)(0) = Crd1; (*Crd)(1) = Crd2; 。。。。。。 ***************************************************** 參數(shù) tag 為該節(jié)點的標(biāo)簽，指定給基類 :DomainComponent(tag,NOD_TAG_Node), NOD_TAG_Node 是一個枚舉值，表明該 DomainComponent 對象是一個節(jié)點類型； ndof 該節(jié)點的自由度，本例中，節(jié)點都為兩個自由度； Crd1, Crd2 為該 2 維節(jié)點的坐標(biāo)，賦于成員變量 Crd ，這是 一個類數(shù)組的數(shù)據(jù)類型，創(chuàng)建了一個含該點坐標(biāo)信息的

20、數(shù)組。 [cpp] view plain copy // 將 4 個節(jié)點對象加入有限元模型中 // 如果兩個 node 對象 tag 相同，則會返回失敗 theDomain->addNode(node1); theDomain->addNode(node2); theDomain->addNode(node3); theDomain->addNode(node4); [cpp] view plain copy // 創(chuàng)建一個彈性材料， 見說明 2 UniaxialMaterial *theMaterial = new ElasticMaterial(1, 3

21、000); 說明 2 ：創(chuàng)建材料對象 ***************************************************** UniaxialMaterial *theMaterial = new ElasticMaterial(1,3000); ***************************************************** UniaxialMaterial 類官方說明： http://opensees.berkeley.edu/OpenSees/api/doxygen2/ht ml/classElasticMaterial.htm

22、l 其中， ElasticMaterial 為 UniaxialMaterial 派生類 意為理想彈性材料 http://opensees.berkeley.edu/wiki/index.php/Elastic_Uniaxi al_Material 構(gòu)造函數(shù) 申明： ***************************************************** ElasticMaterial(int tag, double E, double eta =0.0); **************************************************

23、*** 實現(xiàn)： ***************************************************** ElasticMaterial::ElasticMaterial(int tag, double e, doubleet) :UniaxialMaterial(tag,MAT_TAG_ElasticMaterial), trialStrain(0.0), trialStrainRate(0.0), E(e), eta(et), parameterID(0) { } ********************************************

24、********* 其中，第一個參數(shù) tag 為標(biāo)簽，傳遞給基類構(gòu)造函數(shù)， e 為 彈性模型， et 為材料阻尼比， 默認為 0.[cpp] view plain copy // 創(chuàng)建一個工況，編號為 1 ，暫時未知 LoadPattern *theLoadPattern = new LoadPattern(1); theDomain->addLoadPattern(theLoadPattern); // 暫 時未知這句話什么意思 theLoadPattern->setTimeSeries(new LinearSeries()); // 創(chuàng)建一個節(jié)點荷載向量 Vec

25、tor theLoadValues(2); theLoadValues(0) = 100.0; theLoadValues(1) = -50.0; // 第一個參數(shù) tag 標(biāo)簽，第二個參數(shù)表明施加點荷載的節(jié)點 tag ，第三個參數(shù)是一個向量，表明在第一維度施加 100 個 單位力，第二維度施加反方向 50 單位力 NodalLoad *theLoad = new NodalLoad(1, 4, theLoadValues); // 將 // 如果 NodalLoad 對象加入模型， 1 表示加入的工況編號 theDomain->addNodalLoad(theLoad,

26、 1); new NodalLoad 后的節(jié)點編號未在模型中找到，返回失敗 // 如果 addNodalLoad 第 2 個參數(shù)所表示的工況編號不存在， 返回失敗 這里為了避免內(nèi)存泄漏，也為了使代碼的封裝性更強，我更 改了一部分代碼： [cpp] view plain copy AnalysisModel *theModel = new AnalysisModel(); EquiSolnAlgo *theSolnAlgo = new Linear(); StaticIntegrator *theIntegrator = new LoadControl(1.0, 1, 1.0,

27、1.0); ConstraintHandler *theHandler = new PenaltyConstraintHandler(1.0e8,1.0e8); RCM *theRCM = new RCM(); DOF_Numberer *theNumberer = new DOF_Numberer(*theRCM); BandSPDLinSolver *theSolver = new BandSPDLinLapackSolver(); LinearSOE *theSOE = new BandSPDLinSOE(*theSolver); StaticAnalysis theAn

28、alysis(*theDomain, *theHandler, *theNumberer, *theModel, *theSolnAlgo, *theSOE, *theIntegrator); 改為： [cpp] view plain copy // 分析對象封裝 struct MyStaticAnalysis : public StaticAnalysis ConstraintHandler *pConstraintHandler; DOFNumberer *pDOF_Numberer; AnalysisModel *pAnalysisModel; EquiSolnAl

29、go *pEquiSolnAlgo; LinearSOE *pLinearSOE; StaticIntegrator *pStaticIntegrator; MyStaticAnalysis(Domain *theDomain) : StaticAnalysis(*theDomain, *(pConstraintHandler = new PenaltyConstraintHandler(1.0e8,1.0e8)), *(pDOF_Numberer = new DOF_Numberer(*(new RCM()))), *(pAnalysisModel = ne

30、w AnalysisModel()), *(pEquiSolnAlgo = new Linear()), *(pLinearSOE = new BandSPDLinSOE(*(new BandSPDLinLapackSolver()))), *(pStaticIntegrator = new LoadControl(1.0, 1, 1.0, 1.0))) {} ~MyStaticAnalysis(){ delete pConstraintHandler; delete pLinearSOE; } }; delete pDOF_Numberer; delete pAnalysis

31、Model; delete pEquiSolnAlgo; delete pStaticIntegrator; [cpp] view plain copy // 實例化分析模型對象 MyStaticAnalysis &theAnalysis = *(new MyStaticAnalysis(theDomain)); // perform the analysis & print out the results for the domain int numSteps = 1; theAnalysis.analyze(numSteps); // 釋放分析對象 dele

32、te &theAnalysis; // 模型信息打印 opserr << *theDomain; [cpp] view plain copy // 查看 4 節(jié)點兩個自由度上的位移 Vector const &disp4node = theDomain->getNode(4)->getDisp(); printf("x4: %lf, z4: %lf\n", disp4node[0], disp4node[1]); // 查看 3 個桁 架單元的軸力 Information trussInfo; for(int i=0; i<3; ++i)

33、 { Truss *pTruss = (Truss *)theDomain->getElement(i+1); pTruss->getResponse(2, trussInfo); printf("N%d: %lf\n", i+1, trussInfo.theDouble); } // Domain 類的析構(gòu)會釋放加入 domain 的所有元素， 所以 node 之類的對象不用單獨析構(gòu) delete theDomain; 編譯—— 運行——屏幕輸出：第一自由度位移 0.530094 ，第二自由 度位移 -0.177894 桿件 1 軸力： 43.94 桿件 2 軸力： -57.55 桿件 3 軸力： -55.31 與 sap2000 計算結(jié)果一致： sap2000 模型文件*SDB(v14)和* s2k文件，及修改后的源文件 first.cpp 下載：