when I do some eigenanalysis there's something wrong:
the eigenvalue is numrically undetermined or infinite, but the K is 0 when I output it, I don't know what the problem is? Why the K is 0? Please help.
dear fmk :a problem of eigeanalysis
Moderators: silvia, selimgunay, Moderators
Dear fmk, thank for ur patience
[quote="fmk"]what do you mean the K is 0, if you mean your stifness matrix is all zeros, and you have elements, it means your material properties or section definition is incorrect.
if the matrix is singular, i.e. there are values in it, then you have the problem above or not enough constraints.[/quote]
Dear fmk.I'm a beginner of C++ and opensees. During the course of study I met many problems, since the time is pressing, could you help me to modify this program, it will be an example for future reference, thank you very much and sorry for my English .
this is the code:
#include <stdlib.h>
#include <OPS_Globals.h>
#include <StandardStream.h>
#include <ArrayOfTaggedObjects.h>
#include <math.h>
#include <iostream>
//domain
#include <Domain.h>
#include <Node.h>
#include <ElasticMaterial.h>
#include <Element.h>
#include <Truss.h>
#include <TrilinearBackbone.h>
#include <HystereticBackbone.h>
#include <BackboneMaterial.h>
#include <SP_Constraint.h>
#include <Mp_Constraint.h>
#include <LoadPattern.h>
#include <ElementalLoad.h>
#include <Beam2dUniformLoad.h>
#include <LinearSeries.h>
#include <NodalLoad.h>
#include <StaticAnalysis.h>
#include <AnalysisModel.h>
#include <Linear.h>
#include <PenaltyConstraintHandler.h>
#include <DOF_Numberer.h>
#include <RCM.h>
#include <LoadControl.h>
#include <DisplacementControl.h>
#include <BandSPDLinSOE.h>
#include <BandSPDLinLapackSolver.h>
#include <PathSeries.h>
#include <GroundMotionRecord.h>
#include <UniformExcitation.h>
//analysis
#include <PlainHandler.h>
#include <Newmark.h>
#include <PlainNumberer.h>
#include <NewtonRaphson.h>
#include <CTestNormDispIncr.h>
#include <CTestEnergyIncr.h>
#include <LinearCrdTransf2d.h>
#include <TransformationConstraintHandler.h>
#include <BandGenLinLapackSolver.h>
#include <BandGenLinSOE.h>
#include <BandArpackSolver.h>
#include <BandArpackSOE.h>
#include <StandardEigenAlgo.h>
#include <EigenAlgorithm.h>
#include <EigenAnalysis.h>
#include <EigenIntegrator.h>
#include <EigenSOE.h>
#include <FullGenEigenSolver.h>
#include <EigenSolver.h>
#include <FullGenEigenSOE.h>
#include <SymArpackSolver.h>
#include <FrequencyAlgo.h>
#include <SymArpackSOE.h>
#include <Analysis.h>
#include <StaticIntegrator.h>
#include <TransientIntegrator.h>
#include <TransientAnalysis.h>
#include <DirectIntegrationAnalysis.h>
#include <UmfpackGenLinSolver.h>
#include <UmfpackGenLinSOE.h>
StandardStream sserr;
OPS_Stream *opserrPtr = &sserr;
double ops_Dt = 0;
Domain *ops_TheActiveDomain = 0;
Element *ops_TheActiveElement = 0;
//units
int m=1,N=1,sec=1,sec2=sec*sec;
int KN=1000*N,kg=1,m2=m*m;
// main
int main(int argc, char **argv)
{
Domain *theDomain = new Domain();
//parameters
double g=9.8*m/sec2,q=30*KN/m,h=3*m,Lb=5*m;
double Wb=Lb*q;
//define node
Node *node1=new Node(1,3,0,0);
Node *node2=new Node(2,3,0,h);
Node *node3=new Node(3,3,0,2*h);
Node *node4=new Node(4,3,0,3*h);
Node *node5=new Node(5,3,0,4*h);
Node *node6=new Node(6,3,Lb,0);
Node *node7=new Node(7,3,Lb,h);
Node *node8=new Node(8,3,Lb,2*h);
Node *node9=new Node(9,3,Lb,3*h);
Node *node10=new Node(10,3,Lb,4*h);
Node *node11=new Node(11,3,2*Lb,0);
Node *node12=new Node(12,3,2*Lb,h);
Node *node13=new Node(13,3,2*Lb,2*h);
Node *node14=new Node(14,3,2*Lb,3*h);
Node *node15=new Node(15,3,2*Lb,4*h);
//addnode
theDomain->addNode(node1);
theDomain->addNode(node2);
theDomain->addNode(node3);
theDomain->addNode(node4);
theDomain->addNode(node5);
theDomain->addNode(node6);
theDomain->addNode(node7);
theDomain->addNode(node8);
theDomain->addNode(node9);
theDomain->addNode(node10);
theDomain->addNode(node11);
theDomain->addNode(node12);
theDomain->addNode(node13);
theDomain->addNode(node14);
theDomain->addNode(node15);
//constraint
SP_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, 1, 2, 0.0);
SP_Constraint *sp4 = new SP_Constraint(4, 6, 0, 0.0);
SP_Constraint *sp5 = new SP_Constraint(5, 6, 1, 0.0);
SP_Constraint *sp6 = new SP_Constraint(6, 6, 2, 0.0);
SP_Constraint *sp7 = new SP_Constraint(7, 11, 0, 0.0);
SP_Constraint *sp8 = new SP_Constraint(8, 11, 1, 0.0);
SP_Constraint *sp9 = new SP_Constraint(9, 11, 2, 0.0);
SP_Constraint *sp10= new SP_Constraint(10, 2, 0, 0.0);
SP_Constraint *sp11= new SP_Constraint(11, 2, 1, 0.0);
SP_Constraint *sp12= new SP_Constraint(12, 2, 2, 1.0);
SP_Constraint *sp13= new SP_Constraint(13, 3, 0, 0.0);
SP_Constraint *sp14= new SP_Constraint(14, 3, 1, 0.0);
SP_Constraint *sp15= new SP_Constraint(15, 3, 2, 1.0);
SP_Constraint *sp16= new SP_Constraint(16, 4, 0, 0.0);
SP_Constraint *sp17= new SP_Constraint(17, 4, 1, 0.0);
SP_Constraint *sp18= new SP_Constraint(18, 4, 2, 1.0);
SP_Constraint *sp19= new SP_Constraint(19, 5, 0, 0.0);
SP_Constraint *sp20= new SP_Constraint(20, 5, 1, 0.0);
SP_Constraint *sp21= new SP_Constraint(21, 5, 2, 1.0);
SP_Constraint *sp22= new SP_Constraint(22, 7, 0, 0.0);
SP_Constraint *sp23= new SP_Constraint(23, 7, 1, 0.0);
SP_Constraint *sp24= new SP_Constraint(24, 7, 2, 1.0);
SP_Constraint *sp25= new SP_Constraint(25, 8, 0, 0.0);
SP_Constraint *sp26= new SP_Constraint(26, 8, 1, 0.0);
SP_Constraint *sp27= new SP_Constraint(27, 8, 2, 1.0);
SP_Constraint *sp28= new SP_Constraint(28, 9, 0, 0.0);
SP_Constraint *sp29= new SP_Constraint(29, 9, 1, 0.0);
SP_Constraint *sp30= new SP_Constraint(30, 9, 2, 1.0);
SP_Constraint *sp31= new SP_Constraint(31, 10, 0, 0.0);
SP_Constraint *sp32= new SP_Constraint(32, 10, 1, 0.0);
SP_Constraint *sp33= new SP_Constraint(33, 10, 2, 1.0);
SP_Constraint *sp34= new SP_Constraint(34, 12, 0, 0.0);
SP_Constraint *sp35= new SP_Constraint(35, 12, 1, 0.0);
SP_Constraint *sp36= new SP_Constraint(36, 12, 2, 1.0);
SP_Constraint *sp37= new SP_Constraint(37, 13, 0, 0.0);
SP_Constraint *sp38= new SP_Constraint(38, 13, 1, 0.0);
SP_Constraint *sp39= new SP_Constraint(39, 13, 2, 1.0);
SP_Constraint *sp40= new SP_Constraint(40, 14, 0, 0.0);
SP_Constraint *sp41= new SP_Constraint(41, 14, 1, 0.0);
SP_Constraint *sp42= new SP_Constraint(42, 14, 2, 1.0);
SP_Constraint *sp43= new SP_Constraint(43, 15, 0, 0.0);
SP_Constraint *sp44= new SP_Constraint(44, 15, 1, 0.0);
SP_Constraint *sp45= new SP_Constraint(45, 15, 2, 1.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);
theDomain->addSP_Constraint(sp7);
theDomain->addSP_Constraint(sp8);
theDomain->addSP_Constraint(sp9);
theDomain->addSP_Constraint(sp10);
theDomain->addSP_Constraint(sp11);
theDomain->addSP_Constraint(sp12);
theDomain->addSP_Constraint(sp13);
theDomain->addSP_Constraint(sp14);
theDomain->addSP_Constraint(sp15);
theDomain->addSP_Constraint(sp16);
theDomain->addSP_Constraint(sp17);
theDomain->addSP_Constraint(sp18);
theDomain->addSP_Constraint(sp19);
theDomain->addSP_Constraint(sp20);
theDomain->addSP_Constraint(sp21);
theDomain->addSP_Constraint(sp22);
theDomain->addSP_Constraint(sp23);
theDomain->addSP_Constraint(sp24);
theDomain->addSP_Constraint(sp25);
theDomain->addSP_Constraint(sp26);
theDomain->addSP_Constraint(sp27);
theDomain->addSP_Constraint(sp28);
theDomain->addSP_Constraint(sp29);
theDomain->addSP_Constraint(sp30);
theDomain->addSP_Constraint(sp31);
theDomain->addSP_Constraint(sp32);
theDomain->addSP_Constraint(sp33);
theDomain->addSP_Constraint(sp34);
theDomain->addSP_Constraint(sp35);
theDomain->addSP_Constraint(sp36);
theDomain->addSP_Constraint(sp37);
theDomain->addSP_Constraint(sp38);
theDomain->addSP_Constraint(sp39);
theDomain->addSP_Constraint(sp40);
theDomain->addSP_Constraint(sp41);
theDomain->addSP_Constraint(sp42);
theDomain->addSP_Constraint(sp43);
theDomain->addSP_Constraint(sp44);
theDomain->addSP_Constraint(sp45);
// MP_Constraint(int tag,int nodeRetain,int nodeConstr,Matrix &constrnt,ID &constrainedDOF,ID &retainedDOF);
Matrix constrnt(1,1);
constrnt(0,0)=1;
ID reDOF(1);
reDOF(0)=0;
MP_Constraint *mp1=new MP_Constraint(1,7,2,constrnt, reDOF,reDOF);
MP_Constraint *mp2=new MP_Constraint(2,7,12,constrnt, reDOF,reDOF);
MP_Constraint *mp3=new MP_Constraint(3,8,3,constrnt, reDOF,reDOF);
MP_Constraint *mp4=new MP_Constraint(4,8,13,constrnt, reDOF,reDOF);
MP_Constraint *mp5=new MP_Constraint(5,9,4,constrnt, reDOF,reDOF);
MP_Constraint *mp6=new MP_Constraint(6,9,14,constrnt, reDOF,reDOF);
MP_Constraint *mp7=new MP_Constraint(7,10,5,constrnt, reDOF,reDOF);
MP_Constraint *mp8=new MP_Constraint(8,10,15,constrnt, reDOF,reDOF);
theDomain->addMP_Constraint(mp1);
theDomain->addMP_Constraint(mp2);
theDomain->addMP_Constraint(mp3);
theDomain->addMP_Constraint(mp4);
theDomain->addMP_Constraint(mp5);
theDomain->addMP_Constraint(mp6);
theDomain->addMP_Constraint(mp7);
theDomain->addMP_Constraint(mp8);
//node mass
Matrix mass1(3,3);
mass1(0,0)=0.5*Wb/g;
mass1(0,1)=0;
mass1(0,2)=0;
mass1(1,0)=0;
mass1(1,1)=0.1;
mass1(1,2)=0;
mass1(2,2)=0;
mass1(2,1)=0;
mass1(2,0)=0;
node5->setMass(mass1);
node2->setMass(mass1);
node3->setMass(mass1);
node4->setMass(mass1);
node15->setMass(mass1);
node12->setMass(mass1);
node13->setMass(mass1);
node14->setMass(mass1);
Matrix mass2(3,3);
mass2(0,0)=Wb/g;
mass2(0,1)=0;
mass2(0,2)=0;
mass2(1,0)=0;
mass2(1,1)=0.1;
mass2(1,2)=0;
mass2(2,2)=0;
mass2(2,1)=0;
mass2(2,0)=0;
node10->setMass(mass2);
node7->setMass(mass2);
node8->setMass(mass2);
node9->setMass(mass2);
//total weight
double Weight1=2*q*Lb,Weight3,Weight4;
double Weight2=Weight3=Weight4=Weight1, W=Weight1+Weight2+Weight3+Weight4;
//Fj=W*gi*hi/sum(gi*hi);
double Fi=W*Weight1*h/(Weight1*h*4);
double Load_ce=Fi/3;
//material
HystereticBackbone *theBackBone = new TrilinearBackbone(1,0.009,650,0.03, 800);
UniaxialMaterial *theMaterial=new BackboneMaterial(1, *theBackBone);
//Vector jntOffsetI(2), jntOffsetJ(2);
jntOffsetI.Zero();
jntOffsetJ.Zero();
CrdTransf2d *crdTransf2d;
crdTransf2d = new LinearCrdTransf2d(1, jntOffsetI, jntOffsetJ);
// Truss(int tag,int dimension,int Nd1, int Nd2,UniaxialMaterial &theMaterial,double A, double rho=0.0);
//beam250*500
Element *element1 = new Truss(1, 2, 2, 7, *theMaterial,0.125*m2);
theDomain->addElement(element1);
Element *element2 = new Truss(2, 2, 7, 12,*theMaterial,0.125*m2);
theDomain->addElement(element2);
Element *element3 = new Truss(3, 2, 3, 8, *theMaterial,0.125*m2);
theDomain->addElement(element3);
Element *element4 = new Truss(4,2, 8, 13, *theMaterial,0.125*m2);
theDomain->addElement(element4);
Element *element5 = new Truss(5, 2, 4, 9, *theMaterial,0.125*m2);
theDomain->addElement(element5);
Element *element6 = new Truss(6, 2, 9, 14, *theMaterial,0.125*m2);
theDomain->addElement(element6);
Element *element7 = new Truss(7, 2, 5, 10, *theMaterial,0.125*m2);
theDomain->addElement(element7);
Element *element8 = new Truss(8,2, 10, 15, *theMaterial,0.125*m2);
theDomain->addElement(element8);
//column250*500
Element *c1_element12 = new Truss(12, 2, 1, 2,*theMaterial,0.125*m2);
theDomain->addElement(c1_element12);
Element *c1_element1 = new Truss(13, 2, 2, 3,*theMaterial,0.125*m2);
theDomain->addElement(c1_element1);
Element *c1_element2 = new Truss(14, 2,3, 4,*theMaterial,0.125*m2);
theDomain->addElement(c1_element2);
Element *c1_element3 = new Truss(15, 2, 4, 5,*theMaterial,0.125*m2);
theDomain->addElement(c1_element3);
Element *c2_element4 = new Truss(16, 2, 6, 7,*theMaterial,0.125*m2);
theDomain->addElement(c2_element4);
Element *c2_element5 = new Truss(17, 2, 7, 8,*theMaterial,0.125*m2);
theDomain->addElement(c2_element5);
Element *c3_element6 = new Truss(18, 2, 8, 9,*theMaterial,0.125*m2);
theDomain->addElement(c3_element6);
Element *c3_element7 = new Truss(19, 2, 9, 10,*theMaterial,0.125*m2);
theDomain->addElement(c3_element7);
Element *c1_element8 = new Truss(20, 2, 11, 12,*theMaterial,0.125*m2);
theDomain->addElement(c1_element8);
Element *c1_element9 = new Truss(21, 2, 12, 13,*theMaterial,0.125*m2);
theDomain->addElement(c1_element9);
Element *c1_element10 = new Truss(22, 2, 13,14,*theMaterial,0.125*m2);
theDomain->addElement(c1_element10);
Element *c1_element11 = new Truss(23, 2, 14, 15,*theMaterial,0.125*m2);
theDomain->addElement(c1_element11);
TimeSeries *theSeries = new LinearSeries();
LoadPattern *theLoadPattern = new LoadPattern(1);
theLoadPattern->setTimeSeries(theSeries);
theDomain->addLoadPattern(theLoadPattern);
//nodalload
//NodalLoad(int tag, int node, const Vector &load, bool isLoadConstant = false);
Vector nodalload(3);
nodalload(0)=0;
nodalload(1)=Weight1;
nodalload(2)=0;
NodalLoad *theLoad1=new NodalLoad(1,7,nodalload);
theDomain->addNodalLoad(theLoad1,1);
NodalLoad *theLoad2=new NodalLoad(2,8,nodalload);
theDomain->addNodalLoad(theLoad2,1);
NodalLoad *theLoad3=new NodalLoad(3,9,nodalload);
theDomain->addNodalLoad(theLoad3,1);
NodalLoad *theLoad4=new NodalLoad(4,10,nodalload);
theDomain->addNodalLoad(theLoad4,1);
//gravity and eigenanalysis
//EigenAnalysis(Domain &theDomain,ConstraintHandler &theHandler, DOF_Numberer &theNumberer,AnalysisModel &theModel,
//EigenAlgorithm &theAlgo, EigenSOE &theSOE, EigenIntegrator &theIntegrator);
AnalysisModel *theEigenModel = new AnalysisModel();
RCM *theEigenRCM = new RCM();
DOF_Numberer *theEigenNumberer = new DOF_Numberer(*theEigenRCM);
ConstraintHandler *theEigenHandler = new PenaltyConstraintHandler(1e8,1e8);
SymArpackSolver *eigen_theSolver = new SymArpackSolver();
EigenAlgorithm *eigen_theSolnAlgo = new FrequencyAlgo();
SymArpackSOE *eigen_theSOE=new SymArpackSOE(*eigen_theSolver,*theEigenModel);
EigenIntegrator *eigen_theIntegrator = new EigenIntegrator();
EigenAnalysis eigen_Analysis(*theDomain,*theEigenHandler,*theEigenNumberer, *theEigenModel,*eigen_theSolnAlgo,*eigen_theSOE,*eigen_theIntegrator);
double lambda=eigen_Analysis.analyze(1);
double omega=sqrt(double(lambda));
double k=lambda*W/g;//intial stiff
double T=2*PI/omega;//period
opserr<<"K is "<<k<<"\n";
opserr<<"T is:"<<T<<"\n";
AnalysisModel *theModel = new AnalysisModel();
ConstraintHandler *theHandler = new PenaltyConstraintHandler(1e8, 1e8);
RCM *theRCM = new RCM();
DOF_Numberer *theNumberer = new DOF_Numberer(*theRCM);
CTestNormDispIncr *theTest = new CTestNormDispIncr(1e-8,6,3);
EquiSolnAlgo *theSolnAlgo = new NewtonRaphson(*theTest,0);
StaticIntegrator *theIntegrator = new LoadControl(0.1, 1, 0.1, 0.1);
BandGenLinSolver *theSolver = new BandGenLinLapackSolver();
LinearSOE *theSOE = new BandGenLinSOE(*theSolver);
StaticAnalysis theAnalysis(*theDomain,
*theHandler,
*theNumberer,
*theModel,
*theSolnAlgo,
*theSOE,
*theIntegrator,theTest);
int numSteps = 10;
theAnalysis.analyze(numSteps);
opserr << *theDomain;
opserr<<node5->getDisp();
//node5->getDisp();
theDomain->setLoadConstant();
//groudmotion
double dt=0.01, Tmax=10;
double lamda1=0.5;
const char *filename="BM68elc.acc";
TimeSeries *theaccelSeries = new PathSeries(filename,0.01,1);
GroundMotion *theMotion=new GroundMotionRecord(filename, 1.0,1.0,0.01);
//GroundMotion(0,0,theaccelSeries);
//UniformExcitation(GroundMotion &theMotion,int dof, int tag, double vel0 = 0.0);
LoadPattern *theLoadPattern1 = new UniformExcitation(*theMotion,1, 400,0.0);
theLoadPattern1->setTimeSeries(theaccelSeries);
theDomain->addLoadPattern(theLoadPattern1);
/*************How can I define and get the output file??????????
const ID theDofs[3]={0,1,2};
const ID theNodes[3]={1,2,3};
const char *dataToStore="disp";
const char *path="D:/";
const char *fileName="node.out";
File *file=new File(fileName, dataToStore,1);
file->addFile(fileName,path,dataToStore);
OPS_Stream *theOutputHandler=new DataFileStream(fileName);
Recorder *therecord=new NodeRecorder(*theDofs,theNodes,0,dataToStore,*theDomain,*theOutputHandler);
theDomain->addRecorder(*therecord);
**************/
//POA
//DisplacementControl::DisplacementControl(int node, int dof, double increment, Domain *domain,int numIncr,double min, double max)
AnalysisModel *theModel1 = new AnalysisModel();
ConvergenceTest *test1=new CTestEnergyIncr(1e-8,10,0);
EquiSolnAlgo *theSolnAlgo1 = new NewtonRaphson(*test1);
ConstraintHandler *theHandler1 = new PenaltyConstraintHandler(1e8,1e8);
RCM *theRCM1= new RCM();
DOF_Numberer *theNumberer1 = new DOF_Numberer(*theRCM1);
BandGenLinSolver *theSolver1 = new BandGenLinLapackSolver();
BandGenLinSOE *theSOE1 = new BandGenLinSOE(*theSolver1);
//Newmark(double gamma, double beta, double alphaM, double betaK, double betaKi, double betaKc, bool disp = true);
//the parameters of Newmark would be wrong.
TransientIntegrator *theIntegrator1 = new Newmark(0.5, 0.25, 0.0, 0.0,0.0, 2*0.02/sqrt(lamda1), 1);
DirectIntegrationAnalysis theAnalysis1(*theDomain, *theHandler1,*theNumberer1, *theModel1,*theSolnAlgo1,*theSOE1,*theIntegrator1,test1);
int numSteps2 =100;
theAnalysis1.analyze(numSteps2,0.01);
opserr<<*theDomain;
opserr<<node15->getDisp();
exit (0);
}
if the matrix is singular, i.e. there are values in it, then you have the problem above or not enough constraints.[/quote]
Dear fmk.I'm a beginner of C++ and opensees. During the course of study I met many problems, since the time is pressing, could you help me to modify this program, it will be an example for future reference, thank you very much and sorry for my English .
this is the code:
#include <stdlib.h>
#include <OPS_Globals.h>
#include <StandardStream.h>
#include <ArrayOfTaggedObjects.h>
#include <math.h>
#include <iostream>
//domain
#include <Domain.h>
#include <Node.h>
#include <ElasticMaterial.h>
#include <Element.h>
#include <Truss.h>
#include <TrilinearBackbone.h>
#include <HystereticBackbone.h>
#include <BackboneMaterial.h>
#include <SP_Constraint.h>
#include <Mp_Constraint.h>
#include <LoadPattern.h>
#include <ElementalLoad.h>
#include <Beam2dUniformLoad.h>
#include <LinearSeries.h>
#include <NodalLoad.h>
#include <StaticAnalysis.h>
#include <AnalysisModel.h>
#include <Linear.h>
#include <PenaltyConstraintHandler.h>
#include <DOF_Numberer.h>
#include <RCM.h>
#include <LoadControl.h>
#include <DisplacementControl.h>
#include <BandSPDLinSOE.h>
#include <BandSPDLinLapackSolver.h>
#include <PathSeries.h>
#include <GroundMotionRecord.h>
#include <UniformExcitation.h>
//analysis
#include <PlainHandler.h>
#include <Newmark.h>
#include <PlainNumberer.h>
#include <NewtonRaphson.h>
#include <CTestNormDispIncr.h>
#include <CTestEnergyIncr.h>
#include <LinearCrdTransf2d.h>
#include <TransformationConstraintHandler.h>
#include <BandGenLinLapackSolver.h>
#include <BandGenLinSOE.h>
#include <BandArpackSolver.h>
#include <BandArpackSOE.h>
#include <StandardEigenAlgo.h>
#include <EigenAlgorithm.h>
#include <EigenAnalysis.h>
#include <EigenIntegrator.h>
#include <EigenSOE.h>
#include <FullGenEigenSolver.h>
#include <EigenSolver.h>
#include <FullGenEigenSOE.h>
#include <SymArpackSolver.h>
#include <FrequencyAlgo.h>
#include <SymArpackSOE.h>
#include <Analysis.h>
#include <StaticIntegrator.h>
#include <TransientIntegrator.h>
#include <TransientAnalysis.h>
#include <DirectIntegrationAnalysis.h>
#include <UmfpackGenLinSolver.h>
#include <UmfpackGenLinSOE.h>
StandardStream sserr;
OPS_Stream *opserrPtr = &sserr;
double ops_Dt = 0;
Domain *ops_TheActiveDomain = 0;
Element *ops_TheActiveElement = 0;
//units
int m=1,N=1,sec=1,sec2=sec*sec;
int KN=1000*N,kg=1,m2=m*m;
// main
int main(int argc, char **argv)
{
Domain *theDomain = new Domain();
//parameters
double g=9.8*m/sec2,q=30*KN/m,h=3*m,Lb=5*m;
double Wb=Lb*q;
//define node
Node *node1=new Node(1,3,0,0);
Node *node2=new Node(2,3,0,h);
Node *node3=new Node(3,3,0,2*h);
Node *node4=new Node(4,3,0,3*h);
Node *node5=new Node(5,3,0,4*h);
Node *node6=new Node(6,3,Lb,0);
Node *node7=new Node(7,3,Lb,h);
Node *node8=new Node(8,3,Lb,2*h);
Node *node9=new Node(9,3,Lb,3*h);
Node *node10=new Node(10,3,Lb,4*h);
Node *node11=new Node(11,3,2*Lb,0);
Node *node12=new Node(12,3,2*Lb,h);
Node *node13=new Node(13,3,2*Lb,2*h);
Node *node14=new Node(14,3,2*Lb,3*h);
Node *node15=new Node(15,3,2*Lb,4*h);
//addnode
theDomain->addNode(node1);
theDomain->addNode(node2);
theDomain->addNode(node3);
theDomain->addNode(node4);
theDomain->addNode(node5);
theDomain->addNode(node6);
theDomain->addNode(node7);
theDomain->addNode(node8);
theDomain->addNode(node9);
theDomain->addNode(node10);
theDomain->addNode(node11);
theDomain->addNode(node12);
theDomain->addNode(node13);
theDomain->addNode(node14);
theDomain->addNode(node15);
//constraint
SP_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, 1, 2, 0.0);
SP_Constraint *sp4 = new SP_Constraint(4, 6, 0, 0.0);
SP_Constraint *sp5 = new SP_Constraint(5, 6, 1, 0.0);
SP_Constraint *sp6 = new SP_Constraint(6, 6, 2, 0.0);
SP_Constraint *sp7 = new SP_Constraint(7, 11, 0, 0.0);
SP_Constraint *sp8 = new SP_Constraint(8, 11, 1, 0.0);
SP_Constraint *sp9 = new SP_Constraint(9, 11, 2, 0.0);
SP_Constraint *sp10= new SP_Constraint(10, 2, 0, 0.0);
SP_Constraint *sp11= new SP_Constraint(11, 2, 1, 0.0);
SP_Constraint *sp12= new SP_Constraint(12, 2, 2, 1.0);
SP_Constraint *sp13= new SP_Constraint(13, 3, 0, 0.0);
SP_Constraint *sp14= new SP_Constraint(14, 3, 1, 0.0);
SP_Constraint *sp15= new SP_Constraint(15, 3, 2, 1.0);
SP_Constraint *sp16= new SP_Constraint(16, 4, 0, 0.0);
SP_Constraint *sp17= new SP_Constraint(17, 4, 1, 0.0);
SP_Constraint *sp18= new SP_Constraint(18, 4, 2, 1.0);
SP_Constraint *sp19= new SP_Constraint(19, 5, 0, 0.0);
SP_Constraint *sp20= new SP_Constraint(20, 5, 1, 0.0);
SP_Constraint *sp21= new SP_Constraint(21, 5, 2, 1.0);
SP_Constraint *sp22= new SP_Constraint(22, 7, 0, 0.0);
SP_Constraint *sp23= new SP_Constraint(23, 7, 1, 0.0);
SP_Constraint *sp24= new SP_Constraint(24, 7, 2, 1.0);
SP_Constraint *sp25= new SP_Constraint(25, 8, 0, 0.0);
SP_Constraint *sp26= new SP_Constraint(26, 8, 1, 0.0);
SP_Constraint *sp27= new SP_Constraint(27, 8, 2, 1.0);
SP_Constraint *sp28= new SP_Constraint(28, 9, 0, 0.0);
SP_Constraint *sp29= new SP_Constraint(29, 9, 1, 0.0);
SP_Constraint *sp30= new SP_Constraint(30, 9, 2, 1.0);
SP_Constraint *sp31= new SP_Constraint(31, 10, 0, 0.0);
SP_Constraint *sp32= new SP_Constraint(32, 10, 1, 0.0);
SP_Constraint *sp33= new SP_Constraint(33, 10, 2, 1.0);
SP_Constraint *sp34= new SP_Constraint(34, 12, 0, 0.0);
SP_Constraint *sp35= new SP_Constraint(35, 12, 1, 0.0);
SP_Constraint *sp36= new SP_Constraint(36, 12, 2, 1.0);
SP_Constraint *sp37= new SP_Constraint(37, 13, 0, 0.0);
SP_Constraint *sp38= new SP_Constraint(38, 13, 1, 0.0);
SP_Constraint *sp39= new SP_Constraint(39, 13, 2, 1.0);
SP_Constraint *sp40= new SP_Constraint(40, 14, 0, 0.0);
SP_Constraint *sp41= new SP_Constraint(41, 14, 1, 0.0);
SP_Constraint *sp42= new SP_Constraint(42, 14, 2, 1.0);
SP_Constraint *sp43= new SP_Constraint(43, 15, 0, 0.0);
SP_Constraint *sp44= new SP_Constraint(44, 15, 1, 0.0);
SP_Constraint *sp45= new SP_Constraint(45, 15, 2, 1.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);
theDomain->addSP_Constraint(sp7);
theDomain->addSP_Constraint(sp8);
theDomain->addSP_Constraint(sp9);
theDomain->addSP_Constraint(sp10);
theDomain->addSP_Constraint(sp11);
theDomain->addSP_Constraint(sp12);
theDomain->addSP_Constraint(sp13);
theDomain->addSP_Constraint(sp14);
theDomain->addSP_Constraint(sp15);
theDomain->addSP_Constraint(sp16);
theDomain->addSP_Constraint(sp17);
theDomain->addSP_Constraint(sp18);
theDomain->addSP_Constraint(sp19);
theDomain->addSP_Constraint(sp20);
theDomain->addSP_Constraint(sp21);
theDomain->addSP_Constraint(sp22);
theDomain->addSP_Constraint(sp23);
theDomain->addSP_Constraint(sp24);
theDomain->addSP_Constraint(sp25);
theDomain->addSP_Constraint(sp26);
theDomain->addSP_Constraint(sp27);
theDomain->addSP_Constraint(sp28);
theDomain->addSP_Constraint(sp29);
theDomain->addSP_Constraint(sp30);
theDomain->addSP_Constraint(sp31);
theDomain->addSP_Constraint(sp32);
theDomain->addSP_Constraint(sp33);
theDomain->addSP_Constraint(sp34);
theDomain->addSP_Constraint(sp35);
theDomain->addSP_Constraint(sp36);
theDomain->addSP_Constraint(sp37);
theDomain->addSP_Constraint(sp38);
theDomain->addSP_Constraint(sp39);
theDomain->addSP_Constraint(sp40);
theDomain->addSP_Constraint(sp41);
theDomain->addSP_Constraint(sp42);
theDomain->addSP_Constraint(sp43);
theDomain->addSP_Constraint(sp44);
theDomain->addSP_Constraint(sp45);
// MP_Constraint(int tag,int nodeRetain,int nodeConstr,Matrix &constrnt,ID &constrainedDOF,ID &retainedDOF);
Matrix constrnt(1,1);
constrnt(0,0)=1;
ID reDOF(1);
reDOF(0)=0;
MP_Constraint *mp1=new MP_Constraint(1,7,2,constrnt, reDOF,reDOF);
MP_Constraint *mp2=new MP_Constraint(2,7,12,constrnt, reDOF,reDOF);
MP_Constraint *mp3=new MP_Constraint(3,8,3,constrnt, reDOF,reDOF);
MP_Constraint *mp4=new MP_Constraint(4,8,13,constrnt, reDOF,reDOF);
MP_Constraint *mp5=new MP_Constraint(5,9,4,constrnt, reDOF,reDOF);
MP_Constraint *mp6=new MP_Constraint(6,9,14,constrnt, reDOF,reDOF);
MP_Constraint *mp7=new MP_Constraint(7,10,5,constrnt, reDOF,reDOF);
MP_Constraint *mp8=new MP_Constraint(8,10,15,constrnt, reDOF,reDOF);
theDomain->addMP_Constraint(mp1);
theDomain->addMP_Constraint(mp2);
theDomain->addMP_Constraint(mp3);
theDomain->addMP_Constraint(mp4);
theDomain->addMP_Constraint(mp5);
theDomain->addMP_Constraint(mp6);
theDomain->addMP_Constraint(mp7);
theDomain->addMP_Constraint(mp8);
//node mass
Matrix mass1(3,3);
mass1(0,0)=0.5*Wb/g;
mass1(0,1)=0;
mass1(0,2)=0;
mass1(1,0)=0;
mass1(1,1)=0.1;
mass1(1,2)=0;
mass1(2,2)=0;
mass1(2,1)=0;
mass1(2,0)=0;
node5->setMass(mass1);
node2->setMass(mass1);
node3->setMass(mass1);
node4->setMass(mass1);
node15->setMass(mass1);
node12->setMass(mass1);
node13->setMass(mass1);
node14->setMass(mass1);
Matrix mass2(3,3);
mass2(0,0)=Wb/g;
mass2(0,1)=0;
mass2(0,2)=0;
mass2(1,0)=0;
mass2(1,1)=0.1;
mass2(1,2)=0;
mass2(2,2)=0;
mass2(2,1)=0;
mass2(2,0)=0;
node10->setMass(mass2);
node7->setMass(mass2);
node8->setMass(mass2);
node9->setMass(mass2);
//total weight
double Weight1=2*q*Lb,Weight3,Weight4;
double Weight2=Weight3=Weight4=Weight1, W=Weight1+Weight2+Weight3+Weight4;
//Fj=W*gi*hi/sum(gi*hi);
double Fi=W*Weight1*h/(Weight1*h*4);
double Load_ce=Fi/3;
//material
HystereticBackbone *theBackBone = new TrilinearBackbone(1,0.009,650,0.03, 800);
UniaxialMaterial *theMaterial=new BackboneMaterial(1, *theBackBone);
//Vector jntOffsetI(2), jntOffsetJ(2);
jntOffsetI.Zero();
jntOffsetJ.Zero();
CrdTransf2d *crdTransf2d;
crdTransf2d = new LinearCrdTransf2d(1, jntOffsetI, jntOffsetJ);
// Truss(int tag,int dimension,int Nd1, int Nd2,UniaxialMaterial &theMaterial,double A, double rho=0.0);
//beam250*500
Element *element1 = new Truss(1, 2, 2, 7, *theMaterial,0.125*m2);
theDomain->addElement(element1);
Element *element2 = new Truss(2, 2, 7, 12,*theMaterial,0.125*m2);
theDomain->addElement(element2);
Element *element3 = new Truss(3, 2, 3, 8, *theMaterial,0.125*m2);
theDomain->addElement(element3);
Element *element4 = new Truss(4,2, 8, 13, *theMaterial,0.125*m2);
theDomain->addElement(element4);
Element *element5 = new Truss(5, 2, 4, 9, *theMaterial,0.125*m2);
theDomain->addElement(element5);
Element *element6 = new Truss(6, 2, 9, 14, *theMaterial,0.125*m2);
theDomain->addElement(element6);
Element *element7 = new Truss(7, 2, 5, 10, *theMaterial,0.125*m2);
theDomain->addElement(element7);
Element *element8 = new Truss(8,2, 10, 15, *theMaterial,0.125*m2);
theDomain->addElement(element8);
//column250*500
Element *c1_element12 = new Truss(12, 2, 1, 2,*theMaterial,0.125*m2);
theDomain->addElement(c1_element12);
Element *c1_element1 = new Truss(13, 2, 2, 3,*theMaterial,0.125*m2);
theDomain->addElement(c1_element1);
Element *c1_element2 = new Truss(14, 2,3, 4,*theMaterial,0.125*m2);
theDomain->addElement(c1_element2);
Element *c1_element3 = new Truss(15, 2, 4, 5,*theMaterial,0.125*m2);
theDomain->addElement(c1_element3);
Element *c2_element4 = new Truss(16, 2, 6, 7,*theMaterial,0.125*m2);
theDomain->addElement(c2_element4);
Element *c2_element5 = new Truss(17, 2, 7, 8,*theMaterial,0.125*m2);
theDomain->addElement(c2_element5);
Element *c3_element6 = new Truss(18, 2, 8, 9,*theMaterial,0.125*m2);
theDomain->addElement(c3_element6);
Element *c3_element7 = new Truss(19, 2, 9, 10,*theMaterial,0.125*m2);
theDomain->addElement(c3_element7);
Element *c1_element8 = new Truss(20, 2, 11, 12,*theMaterial,0.125*m2);
theDomain->addElement(c1_element8);
Element *c1_element9 = new Truss(21, 2, 12, 13,*theMaterial,0.125*m2);
theDomain->addElement(c1_element9);
Element *c1_element10 = new Truss(22, 2, 13,14,*theMaterial,0.125*m2);
theDomain->addElement(c1_element10);
Element *c1_element11 = new Truss(23, 2, 14, 15,*theMaterial,0.125*m2);
theDomain->addElement(c1_element11);
TimeSeries *theSeries = new LinearSeries();
LoadPattern *theLoadPattern = new LoadPattern(1);
theLoadPattern->setTimeSeries(theSeries);
theDomain->addLoadPattern(theLoadPattern);
//nodalload
//NodalLoad(int tag, int node, const Vector &load, bool isLoadConstant = false);
Vector nodalload(3);
nodalload(0)=0;
nodalload(1)=Weight1;
nodalload(2)=0;
NodalLoad *theLoad1=new NodalLoad(1,7,nodalload);
theDomain->addNodalLoad(theLoad1,1);
NodalLoad *theLoad2=new NodalLoad(2,8,nodalload);
theDomain->addNodalLoad(theLoad2,1);
NodalLoad *theLoad3=new NodalLoad(3,9,nodalload);
theDomain->addNodalLoad(theLoad3,1);
NodalLoad *theLoad4=new NodalLoad(4,10,nodalload);
theDomain->addNodalLoad(theLoad4,1);
//gravity and eigenanalysis
//EigenAnalysis(Domain &theDomain,ConstraintHandler &theHandler, DOF_Numberer &theNumberer,AnalysisModel &theModel,
//EigenAlgorithm &theAlgo, EigenSOE &theSOE, EigenIntegrator &theIntegrator);
AnalysisModel *theEigenModel = new AnalysisModel();
RCM *theEigenRCM = new RCM();
DOF_Numberer *theEigenNumberer = new DOF_Numberer(*theEigenRCM);
ConstraintHandler *theEigenHandler = new PenaltyConstraintHandler(1e8,1e8);
SymArpackSolver *eigen_theSolver = new SymArpackSolver();
EigenAlgorithm *eigen_theSolnAlgo = new FrequencyAlgo();
SymArpackSOE *eigen_theSOE=new SymArpackSOE(*eigen_theSolver,*theEigenModel);
EigenIntegrator *eigen_theIntegrator = new EigenIntegrator();
EigenAnalysis eigen_Analysis(*theDomain,*theEigenHandler,*theEigenNumberer, *theEigenModel,*eigen_theSolnAlgo,*eigen_theSOE,*eigen_theIntegrator);
double lambda=eigen_Analysis.analyze(1);
double omega=sqrt(double(lambda));
double k=lambda*W/g;//intial stiff
double T=2*PI/omega;//period
opserr<<"K is "<<k<<"\n";
opserr<<"T is:"<<T<<"\n";
AnalysisModel *theModel = new AnalysisModel();
ConstraintHandler *theHandler = new PenaltyConstraintHandler(1e8, 1e8);
RCM *theRCM = new RCM();
DOF_Numberer *theNumberer = new DOF_Numberer(*theRCM);
CTestNormDispIncr *theTest = new CTestNormDispIncr(1e-8,6,3);
EquiSolnAlgo *theSolnAlgo = new NewtonRaphson(*theTest,0);
StaticIntegrator *theIntegrator = new LoadControl(0.1, 1, 0.1, 0.1);
BandGenLinSolver *theSolver = new BandGenLinLapackSolver();
LinearSOE *theSOE = new BandGenLinSOE(*theSolver);
StaticAnalysis theAnalysis(*theDomain,
*theHandler,
*theNumberer,
*theModel,
*theSolnAlgo,
*theSOE,
*theIntegrator,theTest);
int numSteps = 10;
theAnalysis.analyze(numSteps);
opserr << *theDomain;
opserr<<node5->getDisp();
//node5->getDisp();
theDomain->setLoadConstant();
//groudmotion
double dt=0.01, Tmax=10;
double lamda1=0.5;
const char *filename="BM68elc.acc";
TimeSeries *theaccelSeries = new PathSeries(filename,0.01,1);
GroundMotion *theMotion=new GroundMotionRecord(filename, 1.0,1.0,0.01);
//GroundMotion(0,0,theaccelSeries);
//UniformExcitation(GroundMotion &theMotion,int dof, int tag, double vel0 = 0.0);
LoadPattern *theLoadPattern1 = new UniformExcitation(*theMotion,1, 400,0.0);
theLoadPattern1->setTimeSeries(theaccelSeries);
theDomain->addLoadPattern(theLoadPattern1);
/*************How can I define and get the output file??????????
const ID theDofs[3]={0,1,2};
const ID theNodes[3]={1,2,3};
const char *dataToStore="disp";
const char *path="D:/";
const char *fileName="node.out";
File *file=new File(fileName, dataToStore,1);
file->addFile(fileName,path,dataToStore);
OPS_Stream *theOutputHandler=new DataFileStream(fileName);
Recorder *therecord=new NodeRecorder(*theDofs,theNodes,0,dataToStore,*theDomain,*theOutputHandler);
theDomain->addRecorder(*therecord);
**************/
//POA
//DisplacementControl::DisplacementControl(int node, int dof, double increment, Domain *domain,int numIncr,double min, double max)
AnalysisModel *theModel1 = new AnalysisModel();
ConvergenceTest *test1=new CTestEnergyIncr(1e-8,10,0);
EquiSolnAlgo *theSolnAlgo1 = new NewtonRaphson(*test1);
ConstraintHandler *theHandler1 = new PenaltyConstraintHandler(1e8,1e8);
RCM *theRCM1= new RCM();
DOF_Numberer *theNumberer1 = new DOF_Numberer(*theRCM1);
BandGenLinSolver *theSolver1 = new BandGenLinLapackSolver();
BandGenLinSOE *theSOE1 = new BandGenLinSOE(*theSolver1);
//Newmark(double gamma, double beta, double alphaM, double betaK, double betaKi, double betaKc, bool disp = true);
//the parameters of Newmark would be wrong.
TransientIntegrator *theIntegrator1 = new Newmark(0.5, 0.25, 0.0, 0.0,0.0, 2*0.02/sqrt(lamda1), 1);
DirectIntegrationAnalysis theAnalysis1(*theDomain, *theHandler1,*theNumberer1, *theModel1,*theSolnAlgo1,*theSOE1,*theIntegrator1,test1);
int numSteps2 =100;
theAnalysis1.analyze(numSteps2,0.01);
opserr<<*theDomain;
opserr<<node15->getDisp();
exit (0);
}