《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 下載：
OpenSees開發(fā)(二)源碼分析
