nmWTAI-Platform/3rd/meshgen/include/singlePhaseSolver.h

#pragma once
#include <string>
#include <vector>
#include <unordered_set>
#include <algorithm>
namespace nglib {
#include "nglib.h"
}

using namespace nglib;

enum BOUND_TYPE {
    CLOSED,
    CONST_PRESSURE,
    CONST_RATE,
};
enum BOUND_SHAPE {
    CIRCLE,
    RECTANGLE,
    POLYGON
};
enum WELL_TYPE {
    VERTICAL,
    VERTICAL_FRACTURED,
    MULTI_FRACED_HORIZONTAL
};

struct Point {
    double x, y, z;
    std::vector<double> pointData;

    Point(double x0, double y0) {
        x = x0;
        y = y0;
    }

    Point() {};

    // 拷贝构造函数
    Point(const Point& other) {
        //std::cout << "Copy constructor called" << std::endl;
        x = other.x; // 复制数据
        y = other.y;
        z = other.z;
        pointData = other.pointData;
    }

    bool operator==(const Point& other) const {
        const double epsilon = 1e-4; // 误差范围
        return fabs(x - other.x) < epsilon &&
               fabs(y - other.y) < epsilon;
    }
};

struct TimeDis {
    int nLogNum;//每个对数周期的时间点数
    double dTB;//起始时间
    TimeDis() {
        nLogNum = 10; //默认为每个对数周期10个点
        dTB = 1.0e-4; //默认初始时刻为1e-4,如果计算的双对数曲线的井储段未显示，则需要减小该值
    }
};

class Cell {
  public:
    int type;
    std::vector<int> pointIndices;
    std::vector<double> cellData;
    Cell(const Cell& other) {
        type = other.type;
        pointIndices = other.pointIndices;
        cellData = other.cellData;
    }

    Cell() {}

    bool operator==(const Cell&other)const {
        if (type != other.type) {
            return false;
        }

        if (pointIndices.size() != other.pointIndices.size()) {
            return false;
        }

        std::vector<int> a_sorted(pointIndices);
        std::vector<int> b_sorted(other.pointIndices);
        std::sort(a_sorted.begin(), a_sorted.end());
        std::sort(b_sorted.begin(), b_sorted.end());

        return a_sorted == b_sorted;
    }
};

class hash_fun {
  public:
    int operator()(const Cell &A) const {
        return A.type;
    }
};

/*
    class CBoundLine
    {
    //   DECLARE_SERIAL(CBoundLine)

    public:

    CBoundLine() {
        nNodeNum = 8;
        nboundtype = 0;
        ddl = 100;
        bisSel = false;
    };
    virtual ~CBoundLine() {};
    Point pPre, pOri;
    int nboundtype;
    int nNodeNum;
    double ddl;
    int nNodeMin;
    int nNodeMax;
    bool bisSel;

    };*/

class Line {

  public:
    Point p1;
    Point p2;
};

// 裂缝线
class Frac : public Line {
  public:
    Frac(Point p1, Point p2) {
        this->p1 = p1;
        this->p2 = p2;
        dE = 0.01;
        dW = 1;
        dFc = 1000;
        dKr = dFc / dE;
    }

    Frac() {
        dE = 0.01;
        dW = 1;
        dFc = 1000;
        dKr = dFc / dE;
    }

    // 交点
    int nIntersect;                                             // 交点数量
    std::vector<Point> vector_interPoint; // 交点坐标

    //裂缝上的线单元
    std::vector<Cell> fracCells;

    double dE;//裂缝宽度
    double dKr;//裂缝渗透率比,Kf/K
    double dW;//裂缝储能比,(phi*Ct)f/(phi*Ct)
    double dFc;//裂缝导流能力,dFc=dKr*dE;

    //记录netgen的in2d文件中几何点的编号
    //一个裂缝只需要2个节点即可记录
    std::vector<int> ngPointIndex;


};

// 定义外边界的线
class Segment : public Line {
  public:
    Segment(Point p1, Point p2) {
        this->p1 = p1;
        this->p2 = p2;
    }

    Segment() {}

    BOUND_TYPE boundType; // 边界线的边界条件类型
    double dConstPressure;//定义边界时的压力值
    int nNumNodes;      // 边界线上的网格点数
    std::vector<int> vecNodes;//边界线上的网格节点编号
};

// 外边界的基类
class CBoundShape {
  public:
    BOUND_SHAPE boundShape; // 定义外边界类型，0为圆形，1为四边形，2为多边形

    // 表示多边形（包括四边形）的参数
    int nNumSegs;                                           // 边的数量
    std::vector<Segment> vecSegments; // 边

    // 表示圆形
    Point cCenter;
    double dRadius;
    int nNodeNum;
    BOUND_TYPE boundType;
    double dConstPressure;//圆形定压边界时的压力值
    std::vector<int> vecNodes;//圆形边界上的网格点编号

    //记录netgen的in2d文件中几何点的编号
    std::vector<int> ngPointIndex;
};

// 多边形,内部的复合区域
class CBoundPolygon : public CBoundShape {
  public:
    CBoundPolygon() {
        bAnchor = false;
    }

    CBoundPolygon(std::vector<Point> points) { //按拟时针给定的点来构建多边形
        for (int i = 0; i < points.size() - 1; ++i) {
            vecSegments.push_back(Segment(points[i], points[i + 1]));
        }

        vecSegments.push_back(Segment(points[points.size() - 1], points[0]));
        nNumSegs = points.size();
        bAnchor = false;
    }

    int nNumSegs;                                           // 边的数量
    std::vector<Segment> vecSegments; // 边

    double dComKr;//复合区的渗透率比
    double dComW;//复合区的储能比
    bool bAnchor;//是否作为复合区进行材料属性的分区赋值，false表示不进行复合区的材料编号赋值，即其材料编号与背景值相同

    //记录netgen的in2d文件中几何点的编号
    std::vector<int> ngPointIndex;

};

class CBoundWell {
    //    DECLARE_SERIAL(CWell)
  public:
    CBoundWell() {
        nNodeNum = 12;
        //      dRc = 0.1;
        nFlowType = 0;
        nTimeNumQ = 0;
        nTimeNumP = 0;
        dSkin = 0;
        dC = 0.01;
        dMinPress = 0.101325;
        bisSel = false;
        nDenseNum = 10;
        dDenseRadius[0] = 1;
        nDenseNodeNUm[0] = 18;
        dDenseRadius[1] = 10;
        nDenseNodeNUm[1] = 23;
        dDenseRadius[2] = 50;
        nDenseNodeNUm[2] = 22;
        dDenseRadius[3] = 100;
        nDenseNodeNUm[3] = 25;
        dDenseRadius[4] = 500;
        nDenseNodeNUm[4] = 26;
        dDenseRadius[5] = 1000;
        nDenseNodeNUm[5] = 27;
        dDenseRadius[6] = 2000;
        nDenseNodeNUm[6] = 32;
        dDenseRadius[7] = 3000;
        nDenseNodeNUm[7] = 34;
        dDenseRadius[8] = 4000;
        nDenseNodeNUm[8] = 35;
        dDenseRadius[9] = 5000;
        nDenseNodeNUm[9] = 47;
        dDenseRadius[10] = 10000;
        nDenseNodeNUm[10] = 48;
        dDenseRadius[11] = 20000;
        nDenseNodeNUm[11] = 48;
        bisAutoChoose = false;
        bisMultFlow = true;
        nInComNote = 0;
        wellType = VERTICAL;
        dWidth = 0.1;
    };

    virtual ~CBoundWell() {};

    void SetDenseValue();
  public:
    //void Serialize(CArchive& ar);
    std::string sWellname;
    Point pCenter;
    int nFlowType;
    int nNodeNum;
    //double dRc;
    double dSkin;
    double dC;
    double dMinPress;
    int nTimeNumQ;
    int nTimeNumP;
    double pdTimeQ[100];
    double pdTimeP[100];
    double pdQ[100];
    double pdP[100];
    int nNodeMin;
    int nNodeMax;
    bool bisSel;
    int nDenseNum;
    double dDenseRadius[20];
    int nDenseNodeNUm[20];
    double dS;
    bool bisAutoChoose;
    bool bisMultFlow;
    int nInComNote;
    std::vector<int> m_NodeIndex;//all the nodes index of the wellbore
    std::vector<Point> m_PressureData;//计算的井筒压力数据

    WELL_TYPE wellType;
    //直井的几何参数，wellType=VERTICAL
    Point cCenter;
    double dRadius;
    //int nNodeNum;
    BOUND_TYPE boundType;

    //直井压裂井独有的参数，wellType=VERTICAL_FRACTURED
    //直井的中心和半径仍然使用cCenter和dRaduis表示
    double dHf;//裂缝半长
    double dTheata;//裂缝角度(相对于x轴)
    Frac wFrac[2];//通过dHf和dTheata计算的裂缝存储在该变量中，用于后续的计算

    //多段压裂水平井的参数，cCenter表示水平井中心点坐标，wellType=MULTI_FRACED_HORIZONTAL
    double dLength;//长度
    double dBeta;//角度
    int nFracNum;//压裂裂缝的数量
    std::vector<Frac> vecFracs;//压裂裂缝
    double dWidth;//2d模型中水平井宽度
    CBoundPolygon wellPoly;//水平井的边界形状

    //记录netgen的in2d文件中几何点的编号
    std::vector<int> ngPointIndex;
};

/*
    class CBoundCir
    {
    //  DECLARE_SERIAL(CBoundCir)
    public:
    CBoundCir() {
        nNodeNum = 24;
        nboundtype = 0;
        bisSel = false;
        dComKr = 10;
        dComW = 3;
    };
    virtual ~CBoundCir() {};

    public:
    Point pCenter;
    double dRc;
    int nboundtype;
    int nNodeNum;
    bool bisSel;
    int nNodeMin;
    int nNodeMax;
    int nMeshMin;
    int nMeshMax;
    double dComKr;
    double dComW;
    };*/


#ifdef __cplusplus
extern "C" {
#endif

//mesh generation API
__declspec(dllexport) Ng_Mesh* NgMeshGen2D(std::string strFileIn2D);
__declspec(dllexport) Ng_Mesh* NgMeshGen2DByIn2d(std::vector<CBoundWell>& wells, std::vector<CBoundPolygon> limits,
        std::vector<Frac>& faults, std::vector<Frac>& fracs, CBoundShape &shape);
__declspec(dllexport) void Ng2VTK(Ng_Mesh *mesh);


//singlePhase solver API
__declspec(dllexport) int singlePhaseSolver(double* pBaseData, std::vector<Point> points, std::vector<Cell> cells,
        int m_nWellNum, CBoundWell* pwells, std::vector<std::vector<Point>>& vecBP,
        std::vector<std::pair<double, std::vector<double>>>& vecNodePre);
__declspec(dllexport) int singlePhaseSolverNgmesh(const TimeDis time, double* m_pBaseData, nglib::Ng_Mesh* mesh,
        std::vector<CBoundWell>& m_wWell, const std::vector<CBoundPolygon>& comPolygon, const std::vector<Frac>& faultPolygon,
        std::vector<Frac>& fracs, CBoundShape &outBound, std::vector<std::pair<double, std::vector<double>>>& vecNodePre,
        int& timeStep, int& totalSteps);
__declspec(dllexport) int singlePhaseGasSolverNgmesh(double* component, double temperature, const TimeDis time,
        double* m_pBaseData, nglib::Ng_Mesh* mesh, std::vector<CBoundWell>& m_wWell,
        const std::vector<CBoundPolygon>& comPolygon, const std::vector<Frac>& faultPolygon, std::vector<Frac>& fracs,
        CBoundShape &outBound, std::vector<std::pair<double, std::vector<double>>>& vecNodePre, int& timeStep, int& totalSteps);
__declspec(dllexport) bool logLogPre(const std::vector<Point>& pressure, const int& nSection, double* pdTimeQ,
                                     double* pdQ, int nTimeQ, std::vector<Point>& logPre);
#ifdef __cplusplus
}

#endif