I am trying to make it work a little c++ code to solve a ElasticBeam3d element but the system cant solve , here i am providing the code and the program output, so is there any posibility that you can see the error on my code and point out so i can proceed?
//code：
Domain* domain = new Domain();
Node* node1 = new Node(1, 6, 0.0, 0.0, 0.0);
Node* node2 = new Node(2, 6, 5000.0, 0.0, 0.0);
domain->addNode(node1);
domain->addNode(node2);
Vector vecInLocXZPlane(3);
vecInLocXZPlane(0) = 5000.0;
vecInLocXZPlane(1) = 0.0;
vecInLocXZPlane(2) = 0.0;
CrdTransf* crdTransf = new LinearCrdTransf3d(1, vecInLocXZPlane);
int tag = 1;
double width = 250;
double height = 500;
double nu = 0.3;
double E = 30000.0;
double G = E / (2 * (1 + nu));
int Nd1 = 1;
int Nd2 = 2;
double A = width * height;
double Jx = (width * pow(height, 3)) / 12.0;
double Iy = (width * pow(height, 3)) / 12.0;
double Iz = (height * pow(width, 3)) / 12.0;
double J = width * height * height + height * width * width;
double rho = 0.0;
int cMass = 0;
int releasez = 0;
int releasey = 0;
Damping* damping = nullptr;
ElasticBeam3d* element = new ElasticBeam3d(tag, A, E, G, Jx, Iy, Iz, 1, 2, *crdTransf, rho, cMass, releasez, releasey, damping);
domain->addElement(element);
SP_Constraint* sp1 = new SP_Constraint(1, 0, 0.0);
SP_Constraint* sp2 = new SP_Constraint(1, 1, 0.0);
SP_Constraint* sp3 = new SP_Constraint(2, 0, 0.0);
SP_Constraint* sp4 = new SP_Constraint(2, 1, 0.0);
SP_Constraint* sp5 = new SP_Constraint(3, 0, 0.0);
SP_Constraint* sp6 = new SP_Constraint(3, 1, 0.0);
SP_Constraint* sp7 = new SP_Constraint(4, 0, 0.0);
SP_Constraint* sp8 = new SP_Constraint(4, 1, 0.0);
domain->addSP_Constraint(sp1);
domain->addSP_Constraint(sp2);
domain->addSP_Constraint(sp3);
domain->addSP_Constraint(sp4);
domain->addSP_Constraint(sp5);
domain->addSP_Constraint(sp6);
domain->addSP_Constraint(sp7);
domain->addSP_Constraint(sp8);
TimeSeries* theSeries = new LinearSeries();
LoadPattern* theLoadPattern = new LoadPattern(1);
theLoadPattern->setTimeSeries(theSeries);
domain->addLoadPattern(theLoadPattern);
Vector theLoadValues(3);
theLoadValues(0) = 0.0;
theLoadValues(1) = -50.0;
NodalLoad* theLoad = new NodalLoad(1, 2, theLoadValues);
domain->addNodalLoad(theLoad, 1);
AnalysisModel* theModel = new AnalysisModel();
EquiSolnAlgo* theSolnAlgo = new Linear();
StaticIntegrator* theIntegrator = new LoadControl(1.0, 1, 1.0, 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 theAnalysis(*domain,
*theHandler,
*theNumberer,
*theModel,
*theSolnAlgo,
*theSOE,
*theIntegrator);
int numSteps = 1;
int info = theAnalysis.analyze(numSteps);
if (info == 0)
{
opserr << *domain;
cout << "求解成功" << endl;
}
else
{
cout << "求解失败，错误代码为：" << info << endl;
}
Vector disp2node = domain->getNode(2)->getDisp();
cout << "节点2的位移为: " << "Ux=" << disp2node[0] << " Uy=" << disp2node[1] << endl;
// 查看单元的轴力
Information trussInfo;
for (int i = 0; i < 1; ++i)
{
Truss* pTruss = (Truss*)domain->getElement(i + 1);
pTruss->getResponse(2, trussInfo);
double AxialLoad = trussInfo.getData()[0];
cout << "单元" << i + 1 << "的轴力为: " << AxialLoad << endl;
}
// Clean up memory before exit
theAnalysis.clearAll();
domain->clearAll();
delete domain;
delete crdTransf;
system("pause");
return 0;
}
////////////////////////////////////////////////////////////////////////////////////////////////
// here is the program output---------------------------------------------------------------------------------------------------
ElasticBeam3d::setDomain tag: 1 -- Element has zero length