OpenSees開發(fā)(二)源碼分析

OpenSees 開發(fā)（二）源碼分析這是一個平面桁架靜力分析算例，代碼位于OpenSees2.3.0EXAMPLESExample1main.cpp 這里先給 出原始源代碼：cpp view plain copy / standard C+ includes#include<stdlib.h> #include <OPS_Globals.h>#include <StandardStream.h> #include <ArrayOfTaggedObjects.h> / includes for the domain 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 <StaticAnalysis.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 globalvariabled defined in OPS_Globals.h StandardStream sserr; OPS_Stream *opserrPtr = &sserr;doubleops_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/ andbuild the model /Domain *theDomain =new Domain();/ create the nodes usingconstructor:/ Node(tag, ndof, crd1, crd2)/ and then add them to the domainNode*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 Node(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(1, 3000);/create the truss elements using constructor:/Truss(tag, dim, nd1, nd2, Material &,A)/ andthen add them to the domainTruss *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, 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 domainSP_Constraint *sp1 = 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, 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:/ LinearSeries()TimeSeries *theSeries= new LinearSeries();/ construct a loadpattren using constructor:/ LoadPattern(tag)/ and then set its TimeSeries and add it to the domain LoadPattern *theLoadPattern = new LoadPattern(1);theLoadPattern->setTimeSeries(theSeries); theDomain->addLoadPattern(theLoadPattern); / construct a nodal load using constructor:/NodalLoad(tag, nodeID, Vector &)/ firstconstruct a Vector of size 2 and set the values NOTE C INDEXING / then construct the load and add it to the domainVector theLoadValues(2);theLoadValues(0) = 100.0;theLoadValues(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 oftype AnalysisModel,/EquiSolnAlgo of typeLinear/StaticIntegrator of type LoadControl/ ConstraintHandler of type Penalty/DOFNumberer which uses RCM/ LinearSOEof type Band SPD objectAnalysisModel();= new Linear();/ and then the StaticAnalysisAnalysisModelEquiSolnAlgoStaticIntegratornew LoadControl(1.0, 1, 1.0, 1.0);*theModel = new*theSolnAlgo*theIntegrator =ConstraintHandlerRCMDOF_Numberer*theHandler = new PenaltyConstraintHandler(1.0e8,1.0e8);*theRCM = new RCM();*theNumberer = newDOF_Numberer(*theRCM);BandSPDLinSolver*theSolver = new BandSPDLinLapackSolver();LinearSOE*theSOE = newBandSPDLinSOE(*theSolver);StaticAnalysis theAnalysis(*theDomain,*theHandler,*theNumberer,*theModel,*theSolnAlgo,*theSOE,*theIntegrator);/perform the analysis & print out the results for thedomain int numSteps = 1;theAnalysis.analyze(numSteps);opserr <<*theDomain;exit(0); 接下去一步一步解釋代碼：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 = newNode(3, 2, 168.0,0.0); Node *node4 = new Node(4,2, 72.0, 96.0);說明 1 ： Node 構(gòu)造函數(shù)位于 OpenSees2.3.0SRCdomainnodeNode.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);(*Crd)(0) = Crd1;(*Crd)(1) = Crd2;。*參數(shù) tag 為該節(jié)點的標簽，指定給基類 :DomainComponent(tag,NOD_TAG_Node), NOD_TAG_Node 是一個枚舉值，表明該DomainComponent 對象是一個節(jié)點類型；ndof 該節(jié)點的自由度，本例中，節(jié)點都為兩個自由度；Crd1, Crd2 為該 2 維節(jié)點的坐標，賦于成員變量Crd ，這是一個類數(shù)組的數(shù)據(jù)類型，創(chuàng)建了一個含該點坐標信息的數(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, 3000); 說明 2 ：創(chuàng)建材料對象*UniaxialMaterial *theMaterial = newElasticMaterial(1,3000);*UniaxialMaterial 類官方說明：http:/opensees.berkeley.edu/OpenSees/api/doxygen2/html/classElasticMaterial.html其中， ElasticMaterial 為 UniaxialMaterial 派生類意為理想彈性材料http:/opensees.berkeley.edu/wiki/index.php/Elastic_Uniaxial_Material構(gòu)造函數(shù)申明：*ElasticMaterial(int tag, double E, double eta =0.0);*實現(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)*其中，第一個參數(shù)tag 為標簽，傳遞給基類構(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é)點荷載向量Vector theLoadValues(2);theLoadValues(0) = 100.0; theLoadValues(1) = -50.0;/ 第一個參數(shù)tag 標簽，第二個參數(shù)表明施加點荷載的節(jié)點tag ，第三個參數(shù)是一個向量，表明在第一維度施加 100 個單位力，第二維度施加反方向 50 單位力 NodalLoad*theLoad = new NodalLoad(1, 4, theLoadValues);/ 將/ 如果NodalLoad 對象加入模型， 1 表示加入的工況編號theDomain->addNodalLoad(theLoad, 1);new NodalLoad 后的節(jié)點編號未在模型中找到，返回失敗/ 如果 addNodalLoad 第 2 個參數(shù)所表示的工況編號不存在，返回失敗這里為了避免內(nèi)存泄漏，也為了使代碼的封裝性更強，我更改了一部分代碼：cpp view plain copy AnalysisModel*theModel = newAnalysisModel(); EquiSolnAlgo *theSolnAlgo = new Linear(); StaticIntegrator *theIntegrator = newLoadControl(1.0, 1, 1.0, 1.0); ConstraintHandler*theHandler = new PenaltyConstraintHandler(1.0e8,1.0e8);RCM*theRCM = new RCM();DOF_Numberer *theNumberer = newDOF_Numberer(*theRCM);BandSPDLinSolver*theSolver = new BandSPDLinLapackSolver();LinearSOE*theSOE = newBandSPDLinSOE(*theSolver);StaticAnalysistheAnalysis(*theDomain, *theHandler,*theNumberer, *theModel, *theSolnAlgo,*theSOE, *theIntegrator); 改為：cpp view plain copy / 分析對象封裝 structMyStaticAnalysis : public StaticAnalysisConstraintHandler *pConstraintHandler;DOFNumberer*pDOF_Numberer;AnalysisModel*pAnalysisModel;EquiSolnAlgo*pEquiSolnAlgo;LinearSOE*pLinearSOE;StaticIntegrator*pStaticIntegrator;MyStaticAnalysis(Domain *theDomain) :StaticAnalysis(*theDomain,*(pConstraintHandler = newPenaltyConstraintHandler(1.0e8,1.0e8),*(pDOF_Numberer = new DOF_Numberer(*(new RCM(),*(pAnalysisModel = new AnalysisModel(),*(pEquiSolnAlgo = new Linear(),*(pLinearSOE= new BandSPDLinSOE(*(newBandSPDLinLapackSolver(),*(pStaticIntegrator = new LoadControl(1.0, 1, 1.0, 1.0) MyStaticAnalysis()delete pConstraintHandler;delete pLinearSOE; ;delete pDOF_Numberer;delete pAnalysisModel;delete pEquiSolnAlgo;delete pStaticIntegrator;cpp view plain copy / 實例化分析模型對象MyStaticAnalysis &theAnalysis = *(newMyStaticAnalysis(theDomain); / perform the analysis& print out the results for the domain int numSteps =1; theAnalysis.analyze(numSteps); / 釋放分析對象delete &theAnalysis; / 模型信息打印 opserr<< *theDomain;cpp view plain copy / 查看 4 節(jié)點兩個自由度上的位移Vector const &disp4node = theDomain->getNode(4)->getDisp();printf("x4: %lf,z4: %lfn", disp4node0, disp4node1);/ 查看 3 個桁架單元的軸力 Information trussInfo; for(int i=0; i<3; +i) Truss *pTruss = (Truss*)theDomain->getElement(i+1);pTruss->getResponse(2, trussInfo);printf("N%d: %lfn", 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 下載：