封装 C++ 写的 OpenSees 滞回曲线计算程序为 Python 库

#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 <twoNodeLink\TwoNodeLink.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 <TransformationConstraintHandler.h>
#include <DOF_Numberer.h>
#include <RCM.h>
#include <LoadControl.h>
#include <DisplacementControl.h>
#include <BandSPDLinSOE.h>
#include <BandSPDLinLapackSolver.h>

void test_cyclic() {
    // Create the domain
    Domain* theDomain = new Domain();
    Node* node0 = new Node(1, 1, 0.0);
    Node* node1 = new Node(2, 1, 0.0);
    theDomain->addNode(node0);
    theDomain->addNode(node1);
    SP_Constraint* sp1 = new SP_Constraint(1, 0, 0.0);
    theDomain->addSP_Constraint(sp1);
    UniaxialMaterial* theMaterial = new ElasticMaterial(1, 3000);
    int directions[1] = { 0 };
    const ID* id_dir = new ID(directions, 1);
    TwoNodeLink* ele = new TwoNodeLink(1, 1, 1, 2, *id_dir, &theMaterial);
    theDomain->addElement(ele);
    // Load
    TimeSeries* theSeries = new LinearSeries();
    LoadPattern* theLoadPattern = new LoadPattern(1);
    theLoadPattern->setTimeSeries(theSeries);
    theDomain->addLoadPattern(theLoadPattern);
    Vector theLoadValues(1);
    theLoadValues(0) = 7.0;
    NodalLoad* theLoad = new NodalLoad(1, 2, theLoadValues);
    theDomain->addNodalLoad(theLoad, 1);
    // Analysis
    AnalysisModel* theModel = new AnalysisModel();
    EquiSolnAlgo* theSolnAlgo = new Linear();
    StaticIntegrator* theIntegrator = new DisplacementControl(2, 0, 0.1, theDomain, 1, 0.1, 0.1);
    ConstraintHandler* theHandler = new TransformationConstraintHandler();
    RCM* theRCM = new RCM();
    DOF_Numberer* theNumberer = new DOF_Numberer(*theRCM);
    BandSPDLinSolver* theSolver = new BandSPDLinLapackSolver();
    LinearSOE* theSOE = new BandSPDLinSOE(*theSolver);

    StaticAnalysis    theAnalysis(*theDomain,
        *theHandler,
        *theNumberer,
        *theModel,
        *theSolnAlgo,
        *theSOE,
        *theIntegrator);
    // Analyze
    int numSteps = 1;
    theAnalysis.analyze(numSteps);
    // Output
    const char* response = "forces";
    const Vector* theForces = theDomain->getElementResponse(1, &response, 1);
    opserr << *theDomain;
    opserr << theForces[0];

    // Clean up memory before exit
    theAnalysis.clearAll();
    theDomain->clearAll();
    delete theDomain;
    delete theMaterial;
}

ID();
ID(int);
ID(int size, int arraySize);
ID(int *data, int size, bool cleanIt = false);
ID(const ID &);    

int main() {
    cyclic();
}

std::vector<double> cyclic(std::vector<double> protocol, std::vector<double> params)

#include <vector>

UniaxialMaterial* theMaterial = new Steel01(1, params[0], params[1], params[2]);

// Analyze
int protocol_size = protocol.size();
std::vector<double> result(protocol_size);
double disp;
int ok;
bool success = true;
const char* response = "forces";
for (int i = 0; i < protocol_size; i++) {
    disp = (i == 0) ? protocol[0] : protocol[i] - protocol[i - 1];
    StaticIntegrator* theIntegrator = new DisplacementControl(2, 0, disp, theDomain, 1, disp, disp);
    theAnalysis.setIntegrator(* theIntegrator);
    ok = theAnalysis.analyze(1);
    if (ok == 0) {
        const Vector* theForces = theDomain->getElementResponse(1, &response, 1);
        result[i] = theForces[0][1];
    }
    else {
        result[i] = 0.0;
        success = false;
    }
}

int main() {
    double protocol[12] = { 1., 2., 3., 2., 1., 0., -1., -2., -3., -2., -1., 0. };
    std::vector<double> vec_protocol(protocol, protocol + 12);
    double params[3] = { 233., 123., 0.1 };
    std::vector<double> vec_params(params, params + 3);
    std::vector<double> res = cyclic(vec_protocol, vec_params);
    for (int i = 0; i < res.size(); i++) {
        opserr << res[i] << "\n";
    }
}

std::vector<double> cyclic_batch(std::vector<double> protocol, std::vector<double> params)

const int NUM_MAT = 3;
int mat_size = params.size() / NUM_MAT;
for (int i = 0; i < mat_size; i++) {
    UniaxialMaterial* theMaterial = new Steel01(1, params[i * NUM_MAT], params[i * NUM_MAT + 1], params[i * NUM_MAT + 2]);
    TwoNodeLink* ele = new TwoNodeLink(1, 1, 1, 2, *id_dir, &theMaterial);
    theDomain->addElement(ele);
}

std::vector<double> result(protocol_size * mat_size);
if (ok == 0) {
    for (int j = 0; j < mat_size; j++) {
        const Vector* theForces = theDomain->getElementResponse(j + 1, &response, 1);
        result[j * protocol_size + i] = theForces[0][1];
    }
}
else {
    for (int j = 0; j < mat_size; j++) {
        result[j * protocol_size + i] = 0.0;
    }
    success = false;
}

int main() {
    double protocol[12] = { 1, 2, 3, 2, 1, 0, -1, -2, -3, -2, -1, 0 };
    Vector vec_protocol;
    vec_protocol.setData(protocol, 12);
    double params[6] = { 233, 123, 0.1, 12, 13, 0.2 };
    Vector vec_params;
    vec_params.setData(params, 6);
    Vector c = cyclic_batch(vec_protocol, vec_params);
    for (int i = 0; i < c.Size(); i++) {
        opserr << c(i) << "\n";
    }
}

PYBIND11_MODULE(opsutils, m) {
    m.doc() = "Utility OpenSees functions.";
    m.def("cyclic", &cyclic, "Cyclic tests for a single material");
    m.def("cyclic_batch", &cyclic_batch, "Cyclic tests for a batch of materials");
}

#include <fstream>
#include <iostream>
#include <iomanip>

ofstream fout;
fout.open("out.txt");
for (int i = 0; i < res.size(); i++) {
    fout << std::fixed << std::setprecision(14) << res[i] << "\n";
}
fout.close();

py::array_t<double> cyclic2(py::array_t<double> protocol, py::array_t<double> mats) {
    
    if (protocol.ndim() != 1) {
        throw std::runtime_error("Protocol must in dimension 1.");
    }
    auto p = protocol.mutable_unchecked<1>();
    py::ssize_t len_p = p.shape(0);

    py::buffer_info ptc = protocol.request();

    if (mats.ndim() != 2) {
        throw std::runtime_error("Material parameters must in dimension 2.");
    }
    auto m = mats.mutable_unchecked<2>();
    py::ssize_t len_m = m.shape(0);
    py::ssize_t len_param = m.shape(1);
    if (len_param != 3) {
        throw std::runtime_error("Number of parameters must be 3.");
    }
    // Create the domain
    Domain* theDomain = new Domain();
    Node* node0 = new Node(1, 1, 0.0);
    Node* node1 = new Node(2, 1, 0.0);
    theDomain->addNode(node0);
    theDomain->addNode(node1);
    SP_Constraint* sp1 = new SP_Constraint(1, 0, 0.0);
    theDomain->addSP_Constraint(sp1);
    int directions[1] = { 0 };
    const ID* id_dir = new ID(directions, 1);

    std::vector<UniaxialMaterial>* mats = new std::vector<UniaxialMaterial>(len_m);
    for (int i = 0; i < len_m; i++) {
        UniaxialMaterial* theMaterial = new Steel01(i + 1, m(i, 0), m(i, 1), m(i, 2));
        TwoNodeLink* ele = new TwoNodeLink(i + 1, 1, 1, 2, *id_dir, &theMaterial);
        theDomain->addElement(ele);
    }
    // Load
    TimeSeries* theSeries = new LinearSeries();
    LoadPattern* theLoadPattern = new LoadPattern(1);
    theLoadPattern->setTimeSeries(theSeries);
    theDomain->addLoadPattern(theLoadPattern);
    Vector theLoadValues(1);
    theLoadValues(0) = 1.0;
    NodalLoad* theLoad = new NodalLoad(1, 2, theLoadValues);
    theDomain->addNodalLoad(theLoad, 1);
    // Analysis
    AnalysisModel* theModel = new AnalysisModel();
    EquiSolnAlgo* theSolnAlgo = new Linear();
    StaticIntegrator* theIntegrator = new DisplacementControl(2, 0, 1.0, theDomain, 1, 1.0, 1.0); // fake
    ConstraintHandler* theHandler = new TransformationConstraintHandler();
    RCM* theRCM = new RCM();
    DOF_Numberer* theNumberer = new DOF_Numberer(*theRCM);
    BandSPDLinSolver* theSolver = new BandSPDLinLapackSolver();
    LinearSOE* theSOE = new BandSPDLinSOE(*theSolver);

    StaticAnalysis    theAnalysis(*theDomain,
        *theHandler,
        *theNumberer,
        *theModel,
        *theSolnAlgo,
        *theSOE,
        *theIntegrator);

    // Analyze
    auto result = py::array_t<double>({ len_m, len_p });
    auto r = result.mutable_unchecked<2>();
    double disp;
    int ok;
    bool success = true;
    const char* response = "forces";
    for (int i = 0; i < len_p; i++) {
        disp = (i == 0) ? p(i) : p(i) - p(i-1);
        StaticIntegrator* theIntegrator = new DisplacementControl(2, 0, disp, theDomain, 1, disp, disp);
        theAnalysis.setIntegrator(*theIntegrator);
        ok = theAnalysis.analyze(1);
        if (ok == 0) {
            for (int j = 0; j < len_m; j++) {
                const Vector* theForces = theDomain->getElementResponse(j + 1, &response, 1);
                r(j, i) = theForces[0][1];
            }
        }
        else {
            for (int j = 0; j < len_m; j++) {
                r(j, i) = 0.0;
            }
            success = false;
        }
    }

    // Clean up memory before exit
    theAnalysis.clearAll();
    theDomain->clearAll();
    delete theDomain;
    // delete theMaterial;
    return result;
}