OpenSees开发二源码分析.docx

OpenSees开发二源码分析.docx

《OpenSees开发二源码分析.docx》由会员分享，可在线阅读，更多相关《OpenSees开发二源码分析.docx（7页珍藏版）》请在冰点文库上搜索。

OpenSees开发二源码分析

OpenSees开发

（二）源码分析

这是一个平面桁架静力分析算例，代码位于OpenSees2.3.0\EXAMPLES\Example1\main.cpp这里先给出原始源代码：

[cpp]viewplaincopy//standardC++includes#include<stdlib.h>#include<OPS_Globals.h>#include<StandardStream.h>#include<ArrayOfTaggedObjects.h>//includesforthedomainclasses#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>//includesfortheanalysisclasses#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>//inittheglobalvariableddefinedinOPS_Globals.hStandardStreamsserr;OPS_Stream*opserrPtr=&sserr;doubleops_Dt=0;//Domain*ops_TheActiveDomain=0;Element*ops_TheActiveElement=0;//mainroutineintmain（intargc,char**argv）{////nowcreateadomainandamodelbuilder//andbuildthemodel//Domain*theDomain=newDomain（）;//createthenodesusingconstructor:

//Node（tag,ndof,crd1,crd2）//andthenaddthemtothedomainNode*node1=newNode（1,2,0.0,0.0）;Node*node2=newNode（2,2,144.0,0.0）;Node*node3=newNode（3,2,168.0,0.0）;Node*node4=newNode（4,2,72.0,96.0）;theDomain->addNode（node1）;theDomain->addNode（node2）;theDomain->addNode（node3）;theDomain->addNode（node4）;//createanelasticmaterialusingconstriuctor:

//ElasticMaterialModel（tag,E）UniaxialMaterial*theMaterial=newElasticMaterial（1,3000）;//createthetrusselementsusingconstructor:

//Truss（tag,dim,nd1,nd2,Material&,A）//andthenaddthemtothedomainTruss*truss1=newTruss（1,2,1,4,*theMaterial,10.0）;Truss*truss2=newTruss（2,2,2,4,*theMaterial,5.0）;Truss*truss3=newTruss（3,2,3,4,*theMaterial,5.0）;theDomain->addElement（truss1）;theDomain->addElement（truss2）;theDomain->addElement（truss3）;//createthesingle-pointconstraintobjectsusingconstructor:

//SP_Constraint（tag,nodeTag,dofID,value）//andthenaddthemtothedomainSP_Constraint*sp1=newSP_Constraint（1,1,0,0.0）;SP_Constraint*sp2=newSP_Constraint（2,1,1,0.0）;SP_Constraint*sp3=newSP_Constraint（3,2,0,0.0）;SP_Constraint*sp4=newSP_Constraint（4,2,1,0.0）;SP_Constraint*sp5=newSP_Constraint（5,3,0,0.0）;SP_Constraint*sp6=newSP_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）;//constructalineartimeseriesobjectusingconstructor:

//LinearSeries（）TimeSeries*theSeries=newLinearSeries（）;//constructaloadpattrenusingconstructor:

//LoadPattern（tag）//andthensetit'sTimeSeriesandaddittothedomainLoadPattern*theLoadPattern=newLoadPattern

（1）;theLoadPattern->setTimeSeries（theSeries）;theDomain->addLoadPattern（theLoadPattern）;//constructanodalloadusingconstructor:

//NodalLoad（tag,nodeID,Vector&）//firstconstructaVectorofsize2andsetthevaluesNOTECINDEXING//thenconstructtheloadandaddittothedomainVectortheLoadValues

（2）;theLoadValues（0）=100.0;theLoadValues

（1）=-50.0;NodalLoad*theLoad=newNodalLoad（1,4,theLoadValues）;theDomain->addNodalLoad（theLoad,1）;//createanAnalysisobjecttoperformastaticanalysisofthemodel//-constructs:

//AnalysisModeloftypeAnalysisModel,//EquiSolnAlgooftypeLinear//StaticIntegratoroftypeLoadControl//ConstraintHandleroftypePenalty//DOF_NumbererwhichusesRCM//LinearSOEoftypeBandSPD//andthentheStaticAnalysisobjectAnalysisModel*theModel=newAnalysisModel（）;EquiSolnAlgo*theSolnAlgo=newLinear（）;StaticIntegrator*theIntegrator=newLoadControl（1.0,1,1.0,1.0）;ConstraintHandler*theHandler=newPenaltyConstraintHandler（1.0e8,1.0e8）;RCM*theRCM=newRCM（）;DOF_Numberer*theNumberer=newDOF_Numberer（*theRCM）;BandSPDLinSolver*theSolver=newBandSPDLinLapackSolver（）;LinearSOE*theSOE=newBandSPDLinSOE（*theSolver）;StaticAnalysistheAnalysis（*theDomain,*theHandler,*theNumberer,*theModel,*theSolnAlgo,*theSOE,*theIntegrator）;//performtheanalysis&printouttheresultsforthedomainintnumSteps=1;theAnalysis.analyze（numSteps）;opserr<<*theDomain;exit（0）;}接下去一步一步解释代码：

[cpp]viewplaincopy//创建一个有限元模型Domain*theDomain=newDomain（）;

[cpp]viewplaincopy//创建4个节点，详细见说明1Node*node1=newNode（1,2,0.0,0.0）;Node*node2=newNode（2,2,144.0,0.0）;Node*node3=newNode（3,2,168.0,0.0）;Node*node4=newNode（4,2,72.0,96.0）;

说明1：

Node构造函数

位于OpenSees2.3.0\SRC\domain\node\Node.cpp

源码定义如下：

*****************************************************

Node:

:

Node（inttag,intndof,doubleCrd1,doubleCrd2）

:

DomainComponent（tag,NOD_TAG_Node）,

numberDOF（ndof）,theDOF_GroupPtr（0）,

Crd（0）,。

。

。

。

。

。

。

{

Crd=newVector

（2）;

（*Crd）（0）=Crd1;

（*Crd）

（1）=Crd2;

。

。

。

。

。

。

*****************************************************

参数tag为该节点的标签，指定给基类

:

DomainComponent（tag,NOD_TAG_Node）,NOD_TAG_Node是一个枚举值，表明该DomainComponent对象是一个节点类型；

ndof该节点的自由度，本例中，节点都为两个自由度；

Crd1,Crd2为该2维节点的坐标，赋于成员变量Crd，这是一个类数组的数据类型，创建了一个含该点坐标信息的数组。

[cpp]viewplaincopy//将4个节点对象加入有限元模型中//如果两个node对象tag相同，则会返回失败theDomain->addNode（node1）;theDomain->addNode（node2）;theDomain->addNode（node3）;theDomain->addNode（node4）;[cpp]viewplaincopy//创建一个弹性材料，见说明2UniaxialMaterial*theMaterial=newElasticMaterial（1,3000）;说明2：

创建材料对象

*****************************************************

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

构造函数

申明：

*****************************************************

ElasticMaterial（inttag,doubleE,doubleeta=0.0）;

*****************************************************

实现：

*****************************************************

ElasticMaterial:

:

ElasticMaterial（inttag,doublee,doubleet）

:

UniaxialMaterial（tag,MAT_TAG_ElasticMaterial）,

trialStrain（0.0）,trialStrainRate（0.0）,

E（e）,eta（et）,parameterID（0）

{

}

*****************************************************

其中，第一个参数tag为标签，传递给基类构造函数，e为弹性模型，et为材料阻尼比，默认为0.[cpp]viewplaincopy//创建一个工况，编号为1，暂时未知LoadPattern*theLoadPattern=newLoadPattern

（1）;theDomain->addLoadPattern（theLoadPattern）;//暂时未知这句话什么意思theLoadPattern->setTimeSeries（newLinearSeries（））;//创建一个节点荷载向量VectortheLoadValues

（2）;theLoadValues（0）=100.0;theLoadValues

（1）=-50.0;//第一个参数tag标签，第二个参数表明施加点荷载的节点tag，第三个参数是一个向量，表明在第一维度施加100个单位力，第二维度施加反方向50单位力NodalLoad*theLoad=newNodalLoad（1,4,theLoadValues）;//将NodalLoad对象加入模型，1表示加入的工况编号theDomain->addNodalLoad（theLoad,1）;//如果newNodalLoad后的节点编号未在模型中找到，返回失败//如果addNodalLoad第2个参数所表示的工况编号不存在，返回失败

这里为了避免内存泄漏，也为了使代码的封装性更强，我更改了一部分代码：

[cpp]viewplaincopyAnalysisModel*theModel=newAnalysisModel（）;EquiSolnAlgo*theSolnAlgo=newLinear（）;StaticIntegrator*theIntegrator=newLoadControl（1.0,1,1.0,1.0）;ConstraintHandler*theHandler=newPenaltyConstraintHandler（1.0e8,1.0e8）;RCM*theRCM=newRCM（）;DOF_Numberer*theNumberer=newDOF_Numberer（*theRCM）;BandSPDLinSolver*theSolver=newBandSPDLinLapackSolver（）;LinearSOE*theSOE=newBandSPDLinSOE（*theSolver）;StaticAnalysistheAnalysis（*theDomain,*theHandler,*theNumberer,*theModel,*theSolnAlgo,*theSOE,*theIntegrator）;改为：

[cpp]viewplaincopy//分析对象封装structMyStaticAnalysis:

publicStaticAnalysis{ConstraintHandler*pConstraintHandler;DOF_Numberer*pDOF_Numberer;AnalysisModel*pAnalysisModel;EquiSolnAlgo*pEquiSolnAlgo;LinearSOE*pLinearSOE;StaticIntegrator*pStaticIntegrator;MyStaticAnalysis（Domain*theDomain）:

StaticAnalysis（*theDomain,*（pConstraintHandler=newPenaltyConstraintHandler（1.0e8,1.0e8））,*（pDOF_Numberer=newDOF_Numberer（*（newRCM（））））,*（pAnalysisModel=newAnalysisModel（））,*（pEquiSolnAlgo=newLinear（））,*（pLinearSOE=newBandSPDLinSOE（*（newBandSPDLinLapackSolver（））））,*（pStaticIntegrator=newLoadControl（1.0,1,1.0,1.0）））{}~MyStaticAnalysis（）{deletepConstraintHandler;deletepDOF_Numberer;deletepAnalysisModel;deletepEquiSolnAlgo;deletepLinearSOE;deletepStaticIntegrator;}};

[cpp]viewplaincopy//实例化分析模型对象MyStaticAnalysis&theAnalysis=*（newMyStaticAnalysis（theDomain））;//performtheanalysis&printouttheresultsforthedomainintnumSteps=1;theAnalysis.analyze（numSteps）;//释放分析对象delete&theAnalysis;//模型信息打印opserr<<*theDomain;

[cpp]viewplaincopy//查看4节点两个自由度上的位移Vectorconst&disp4node=theDomain->getNode（4）->getDisp（）;printf（"x4:

%lf,z4:

%lf\n",disp4node[0],disp4node[1]）;//查看3个桁架单元的轴力InformationtrussInfo;for（inti=0;i<3;++i）{Truss*pTruss=（Truss*）theDomain->getElement（i+1）;pTruss->getResponse（2,trussInfo）;printf（"N%d:

%lf\n",i+1,trussInfo.theDouble）;}//Domain类的析构会释放加入domain的所有元素，所以node之类的对象不用单独析构deletetheDomain;编译——运行——屏幕输出：

第一自由度位移0.530094，第二自由度位移-0.177894

杆件1轴力：

43.94

杆件2轴力：

-57.55

杆件3轴力：

-55.31与sap2000计算结果一致：

sap2000模型文件*.SDB（v14）和*.s2k文件，及修改后的源文件first.cpp下载：