nmWTAI-Platform/3rd/meshgen/example/solver_example-3-test.cpp

// netgen-2010.cpp : 定义控制台应用程序的入口点。
//

#include <stdafx.h>
#include "../include/singlePhaseSolver.h"
#include "windows.h"
#include <iostream>
#include <fstream>
#include <vector>
#include <string>
#include <sstream>
#include <cstdlib>
#include <algorithm>
#include <stdexcept>
#include <iterator>
using namespace std;

// 添加全局变量用于进度跟踪
int g_timeStep = 0;
int g_totalSteps = 100;
volatile bool g_calculationRunning = true;
UINT_PTR g_timerID = 0;

// 添加一个进度条显示函数
void ShowProgressBar(int timeStep, int totalSteps, int barWidth = 50) {
    float progress = (float)timeStep / totalSteps;
    int pos = barWidth * progress;

    std::cout << "[";

    for (int i = 0; i < barWidth; ++i) {
        if (i < pos) {
            std::cout << "=";
        } else if (i == pos) {
            std::cout << ">";
        } else {
            std::cout << " ";
        }
    }

    std::cout << "] " << int(progress * 100.0) << "% (" << timeStep << "/" << totalSteps << ")\r";
    std::cout.flush();
}

// 添加定时器回调函数
VOID CALLBACK ProgressTimerProc(HWND hwnd, UINT uMsg, UINT_PTR idEvent, DWORD dwTime) {
    if (g_calculationRunning && g_totalSteps > 0) {
        ShowProgressBar(g_timeStep, g_totalSteps);
    }
}

// 添加油藏类型枚举
enum RESERVOIR_TYPE {
    OIL_RESERVOIR = 0, // 油藏
    GAS_RESERVOIR = 1  // 气藏
};

void Ng2VTKWithData(Ng_Mesh *mesh, std::vector<double> vecField, std::string filename) {
    std::cout << "Convert Ng_Mesh to vtk legacy file version of 4.2" << std::endl;
    int i, nseg, ne, np, matnum;
    int nodes[3];
    double point[2];
    //    int *EleNodeNum;
    //    double *NodeCoors;

    np = Ng_GetNP_2D(mesh);
    ne = Ng_GetNE_2D(mesh);
    nseg = Ng_GetNSeg_2D(mesh);

    std::ofstream vtkFile;
    vtkFile.open(filename);

    if (!vtkFile.is_open()) {
        std::cerr << "Error opening file: " << "mesh.vtk" << std::endl;
        return;
    }

    // Write VTK file header
    vtkFile << "# vtk DataFile Version 4.1" << std::endl;
    vtkFile << "Mesh data" << std::endl;
    vtkFile << "ASCII" << std::endl;
    vtkFile << "DATASET UNSTRUCTURED_GRID" << std::endl;

    // Write points
    vtkFile << "POINTS " << np << " float" << std::endl;

    for (i = 1; i <= np; i++) {
        Ng_GetPoint_2D(mesh, i, point);
        vtkFile << point[0] << " " << point[1] << " " << 0 << std::endl;
    }

    // Write cells
    vtkFile << "CELLS " << ne << " " << ne * 4 << std::endl;
    std::vector<int> vecMater;

    for (int i = 1; i <= ne; i++) {
        Ng_GetElement_2D(mesh, i, nodes, &matnum);
        vecMater.emplace_back(matnum);
        vtkFile << 3 << " " << nodes[0] - 1 << " " << nodes[1] - 1 << " " << nodes[2] - 1 << std::endl;
    }

    // Write cell types
    vtkFile << "CELL_TYPES " << ne << std::endl;

    for (int i = 1; i <= ne; i++) {
        vtkFile << 5 << std::endl; // Assuming getCellType is a method that returns the VTK cell type
    }

    // Write cell data
    vtkFile << "CELL_DATA " << ne << std::endl;
    vtkFile << "SCALARS Material int 1" << std::endl;
    vtkFile << "LOOKUP_TABLE default" << std::endl;

    for (int i = 0; i < vecMater.size(); i++) {
        vtkFile << vecMater[i] << std::endl;
    }

    // Write point data
    vtkFile << "POINT_DATA " << np << std::endl;
    vtkFile << "SCALARS Pressure double 1" << std::endl;
    vtkFile << "LOOKUP_TABLE default" << std::endl;

    for (int i = 0; i < vecField.size(); i++) {
        vtkFile << vecField[i] << std::endl;
    }

    vtkFile.close();
}

// 辅助函数
string trim(const string &s) {
    size_t start = s.find_first_not_of(" \t");
    size_t end = s.find_last_not_of(" \t");
    return (start == string::npos) ? "" : s.substr(start, end - start + 1);
}

vector<string> split(const string &s, char delimiter) {
    vector<string> tokens;
    string token;
    stringstream iss(s);

    while (getline(iss, token, delimiter)) {
        if (!token.empty()) {
            tokens.push_back(token);
        }
    }

    return tokens;
}

// 边界类型字符串转枚举值的辅助函数
BOUND_TYPE StringToBoundType(const string &s) {
    string str = trim(s);
    transform(str.begin(), str.end(), str.begin(), ::toupper);

    if (str == "CONST_PRESSURE") {
        return CONST_PRESSURE;
    }

    if (str == "CONST_RATE") {
        return CONST_RATE;
    }

    if (str == "CLOSED") {
        return CLOSED;
    }

    throw runtime_error("未知边界条件类型: " + str);
}

// 井类型字符串转枚举值的辅助函数
WELL_TYPE StringToWellType(const string &s) {
    string str = trim(s);
    transform(str.begin(), str.end(), str.begin(), ::toupper);

    if (str == "VERTICAL") {
        return VERTICAL;
    }

    if (str == "VERTICAL_FRACTURED") {
        return VERTICAL_FRACTURED;
    }

    if (str == "MULTI_FRACED_HORIZONTAL") {
        return MULTI_FRACED_HORIZONTAL;
    }

    throw runtime_error("未知井类型: " + str);
}

void ReadParameters(const string &filename,
                    CBoundShape &shape,
                    vector<CBoundWell>& wells,
                    vector<CBoundPolygon>& limits,
                    vector<Frac>& faults,
                    vector<Frac>& fracs,
                    double pBaseData[8],
                    TimeDis &timeDis,
                    RESERVOIR_TYPE &reservoirType,
                    double& reservoirTemp,
                    double gasComponents[32]) {
    ifstream fin(filename.c_str());

    if (!fin) {
        throw runtime_error("无法打开文件: " + filename);
    }

    string line, section;
    int currentWell = -1;
    vector<CBoundWell> tempWells;
    vector<CBoundPolygon> tempLimits;
    CBoundPolygon *currentLimit = nullptr;
    vector<Frac> tempFaults, tempFracs;

    // 设置默认值
    reservoirType = OIL_RESERVOIR; // 默认为油藏
    reservoirTemp = 43;            // 默认43℃

    for (int i = 0; i < 32; i++) {
        gasComponents[i] = 0.0; // 默认所有组分为0
    }

    // 设置TimeDis默认值
    timeDis.nLogNum = 20; // 默认日志数量
    timeDis.dTB = 1.0;    // 默认时间基准值

    while (getline(fin, line)) {
        line = trim(line);

        if (line.empty() || line[0] == '/') { // 跳过空行和注释行
            continue;
        }

        if (line[0] == '[') {
            section = line.substr(1, line.find(']') - 1);
            transform(section.begin(), section.end(), section.begin(), ::toupper);
            currentWell = -1;
            currentLimit = nullptr;
            continue;
        }

        vector<string> tokens = split(line, ' ');

        if (tokens.empty()) {
            continue;
        }

        try {
            if (section == "RESERVOIRBASEDATA") {
                if (tokens.size() == 9) { // 基础参数 + 油藏类型
                    for (int i = 0; i < 8; ++i) {
                        pBaseData[i] = atof(tokens[i].c_str());
                    }

                    // 读取油藏类型
                    reservoirType = (RESERVOIR_TYPE)atoi(tokens[8].c_str());
                } else if (tokens.size() == 8) { // 仅基础参数，默认为油藏
                    for (int i = 0; i < 8; ++i) {
                        pBaseData[i] = atof(tokens[i].c_str());
                    }

                    reservoirType = OIL_RESERVOIR; // 默认为油藏
                } else {
                    throw runtime_error("油藏基础参数个数不正确");
                }
            } else if (section == "GASRESERVOIRDATA" && reservoirType == GAS_RESERVOIR) {
                if (tokens[0] == "Temperature" && tokens.size() >= 2) {
                    reservoirTemp = atof(tokens[1].c_str());
                } else if (tokens[0] == "Components") {
                    // 等待下一行读取组分数据
                } else if (tokens.size() >= 32) { // 气体组分数据，应有32个值
                    // 确保有足够的数字用于32个组分
                    int compCount = min(32, (int)tokens.size());

                    for (int i = 0; i < compCount; ++i) {
                        gasComponents[i] = atof(tokens[i].c_str());
                    }
                }
            } else if (section == "BOUNDSHAPE") {
                if (tokens[0] == "type") {
                    if (tokens.size() < 2) {
                        throw runtime_error("缺少边界类型");
                    }

                    if (tokens[1] == "CIRCLE") {
                        shape.boundShape = CIRCLE;
                    } else if (tokens[1] == "POLYGON") {
                        shape.boundShape = POLYGON;
                    } else {
                        throw runtime_error("未知边界类型: " + tokens[1]);
                    }
                } else if (shape.boundShape == CIRCLE) {
                    // 处理圆形边界参数
                    if (tokens[0] == "center" && tokens.size() >= 3) {
                        shape.cCenter.x = atof(tokens[1].c_str());
                        shape.cCenter.y = atof(tokens[2].c_str());
                    } else if (tokens[0] == "radius" && tokens.size() >= 2) {
                        shape.dRadius = atof(tokens[1].c_str());
                    } else if (tokens[0] == "boundType" && tokens.size() >= 2) {
                        shape.boundType = StringToBoundType(tokens[1]);
                    } else if (tokens[0] == "dConstPressure" && tokens.size() >= 2) {
                        shape.dConstPressure = atof(tokens[1].c_str());
                    }
                } else if (shape.boundShape == POLYGON) {
                    // 处理多边形边界参数
                    if (tokens[0] == "segments" && tokens.size() >= 2) {
                        shape.nNumSegs = atoi(tokens[1].c_str());
                    } else if (tokens.size() >= 5) { // 至少需要4个坐标值和1个边界类型
                        Segment seg;
                        seg.p1 = Point(atof(tokens[0].c_str()), atof(tokens[1].c_str()));
                        seg.p2 = Point(atof(tokens[2].c_str()), atof(tokens[3].c_str()));
                        seg.boundType = StringToBoundType(tokens[4]);

                        // 如果是定压力边界，还需要读取压力值
                        if (tokens.size() >= 6 && seg.boundType == CONST_PRESSURE) {
                            seg.dConstPressure = atof(tokens[5].c_str());
                        }

                        shape.vecSegments.push_back(seg);
                    }
                }
            } else if (section == "WELLS") {
                // 开始读取一个新井
                if (tokens.size() >= 4 && currentWell == -1) {
                    CBoundWell well;
                    well.wellType = StringToWellType(tokens[0]);

                    // 根据井类型读取不同参数
                    if (well.wellType == VERTICAL) {
                        if (tokens.size() < 4) {
                            throw runtime_error("垂直井参数不足");
                        }

                        well.cCenter.x = atof(tokens[1].c_str());
                        well.cCenter.y = atof(tokens[2].c_str());
                        well.dRadius = atof(tokens[3].c_str());
                        tempWells.push_back(well);
                        currentWell = tempWells.size() - 1;
                    } else if (well.wellType == VERTICAL_FRACTURED) {
                        // 垂直压裂井参数
                        if (tokens.size() < 6) {
                            throw runtime_error("垂直压裂井参数不足");
                        }

                        well.cCenter.x = atof(tokens[1].c_str());
                        well.cCenter.y = atof(tokens[2].c_str());
                        well.dRadius = atof(tokens[3].c_str());
                        well.dBeta = atof(tokens[4].c_str());
                        well.dHf = atof(tokens[5].c_str());
                        tempWells.push_back(well);
                        currentWell = tempWells.size() - 1;
                    } else if (well.wellType == MULTI_FRACED_HORIZONTAL) {
                        // 多段压裂水平井参数，只需要读取cCenter, dRadius, dBeta, dLength, nFracNum, dHf
                        if (tokens.size() < 7) {
                            throw runtime_error("多段压裂水平井参数不足");
                        }

                        well.cCenter.x = atof(tokens[1].c_str());
                        well.cCenter.y = atof(tokens[2].c_str());
                        well.dRadius = atof(tokens[3].c_str());
                        well.dBeta = atof(tokens[4].c_str());
                        well.dLength = atof(tokens[5].c_str());
                        well.nFracNum = atoi(tokens[6].c_str());

                        if (tokens.size() >= 8) {
                            well.dHf = atof(tokens[7].c_str());
                        } else {
                            well.dHf = 100.0;    // 默认裂缝半长
                        }

                        // 初始化裂缝数组
                        if (well.nFracNum > 0) {
                            well.vecFracs.resize(2 * well.nFracNum);
                            // 裂缝的具体位置会在模型中自动计算，这里不需要读取
                        }

                        tempWells.push_back(well);
                        currentWell = tempWells.size() - 1;
                    }
                } else if (currentWell >= 0) {
                    CBoundWell &w = tempWells[currentWell];

                    if (tokens[0] == "bisMultFlow" && tokens.size() >= 2) {
                        w.bisMultFlow = (atoi(tokens[1].c_str()) != 0);
                    } else if (tokens[0] == "dSkin" && tokens.size() >= 2) {
                        w.dSkin = atof(tokens[1].c_str());
                    } else if (tokens[0] == "dC" && tokens.size() >= 2) {
                        w.dC = atof(tokens[1].c_str());
                    } else if (tokens[0] == "TimeQ" && tokens.size() >= 2) {
                        int numTimes = atoi(tokens[1].c_str());
                        w.nTimeNumQ = numTimes;
                        // w.pdTimeQ = new double[numTimes];
                        // w.pdQ = new double[numTimes];
                    } else if (tokens.size() == 2 && w.pdTimeQ != nullptr && w.pdQ != nullptr) {
                        // 这里假设TimeQ后的每行数据都是时间和流量对
                        static int timeIdx = 0;

                        if (timeIdx < w.nTimeNumQ) {
                            w.pdTimeQ[timeIdx] = atof(tokens[0].c_str());
                            w.pdQ[timeIdx] = atof(tokens[1].c_str());
                            timeIdx++;

                            // 如果读完了所有时间流量对，重置计数器
                            if (timeIdx == w.nTimeNumQ) {
                                timeIdx = 0;
                            }
                        }
                    } else if (tokens[0] == "FracFc" && w.wellType == VERTICAL_FRACTURED) {
                        if (tokens.size() - 1 != 2) {
                            throw runtime_error("导流能力数量不匹配");
                        }

                        for (int i = 0; i < 2; ++i) {
                            w.wFrac[i].dFc = atof(tokens[i + 1].c_str());
                        }
                    } else if (tokens[0] == "FracFc" && w.wellType == MULTI_FRACED_HORIZONTAL) {
                        if (tokens.size() - 1 != w.nFracNum) {
                            throw runtime_error("导流能力数量不匹配");
                        }

                        for (int i = 0; i < w.nFracNum; ++i) {
                            w.vecFracs[2 * i].dFc = atof(tokens[i + 1].c_str());
                            w.vecFracs[2 * i + 1].dFc = atof(tokens[i + 1].c_str());
                        }
                    }
                }
            } else if (section == "LIMITS") {
                // 修改LIMITS部分，支持更多参数
                if (tokens[0] == "segments" && tokens.size() >= 2) {
                    tempLimits.push_back(CBoundPolygon());
                    currentLimit = &tempLimits.back();
                    currentLimit->nNumSegs = atoi(tokens[1].c_str());
                } else if (currentLimit && tokens[0] == "bAnchor" && tokens.size() >= 2) {
                    currentLimit->bAnchor = (atoi(tokens[1].c_str()) != 0);
                } else if (currentLimit && tokens[0] == "dComKr" && tokens.size() >= 2) {
                    currentLimit->dComKr = atof(tokens[1].c_str());
                } else if (currentLimit && tokens[0] == "dComW" && tokens.size() >= 2) {
                    currentLimit->dComW = atof(tokens[1].c_str());
                } else if (currentLimit && tokens.size() == 4) { // 复合区的线段只需要4个坐标值
                    Segment seg;
                    seg.p1 = Point(atof(tokens[0].c_str()), atof(tokens[1].c_str()));
                    seg.p2 = Point(atof(tokens[2].c_str()), atof(tokens[3].c_str()));
                    // 复合区不需要指定边界类型
                    currentLimit->vecSegments.push_back(seg);
                }
                // 向后兼容旧格式（前两个数字分别表示nNumSegs和bAnchor）
                else if (tokens.size() == 2 && isdigit(tokens[0][0]) && isdigit(tokens[1][0])) {
                    tempLimits.push_back(CBoundPolygon());
                    currentLimit = &tempLimits.back();
                    currentLimit->nNumSegs = atoi(tokens[0].c_str());
                    currentLimit->bAnchor = (atoi(tokens[1].c_str()) != 0);
                }
            } else if (section == "FAULTS") {
                if (tokens.size() == 4) {
                    Frac f;
                    f.p1 = Point(atof(tokens[0].c_str()), atof(tokens[1].c_str()));
                    f.p2 = Point(atof(tokens[2].c_str()), atof(tokens[3].c_str()));
                    tempFaults.push_back(f);
                }
            } else if (section == "FRACS") {
                if (tokens.size() >= 4) {
                    Frac f;
                    f.p1 = Point(atof(tokens[0].c_str()), atof(tokens[1].c_str()));
                    f.p2 = Point(atof(tokens[2].c_str()), atof(tokens[3].c_str()));

                    // 读取储能比和导流能力
                    if (tokens.size() >= 5) {
                        f.dW = atof(tokens[4].c_str());
                    } else {
                        f.dW = 1.0;    // 默认储能比为1.0
                    }

                    if (tokens.size() >= 6) {
                        f.dFc = atof(tokens[5].c_str());
                    } else {
                        f.dFc = 100.0;    // 默认导流能力为100.0
                    }

                    tempFracs.push_back(f);
                }
            } else if (section == "TIMEDIS") {
                if (tokens[0] == "nLogNum" && tokens.size() >= 2) {
                    timeDis.nLogNum = atoi(tokens[1].c_str());
                } else if (tokens[0] == "dTB" && tokens.size() >= 2) {
                    timeDis.dTB = atof(tokens[1].c_str());
                }
            }
        } catch (const exception &e) {
            throw runtime_error("解析错误[" + section + "]: " + e.what() + ", 行: " + line);
        }
    }

    // 最终数据校验
    if (shape.boundShape == CIRCLE) {
        if (shape.dRadius <= 0) {
            throw runtime_error("圆形边界半径必须大于0");
        }
    } else if (shape.boundShape == POLYGON) {
        if (shape.vecSegments.size() != shape.nNumSegs) {
            throw runtime_error("多边形边界线段数量不匹配");
        }
    }

    // 给每个复合区设置默认值(如果未指定)，并验证线段数量
    for (int i = 0; i < tempLimits.size(); ++i) {
        if (tempLimits[i].dComKr <= 0) {
            tempLimits[i].dComKr = 1.0;    // 默认渗透率比为1.0
        }

        if (tempLimits[i].dComW <= 0) {
            tempLimits[i].dComW = 1.0;    // 默认孔隙度比为1.0
        }

        // 验证线段数量
        if (tempLimits[i].vecSegments.size() != tempLimits[i].nNumSegs) {
            throw runtime_error("复合区线段数量与指定数量不匹配");
        }
    }

    // 对气藏参数进行验证
    if (reservoirType == GAS_RESERVOIR) {

        // 检查气体组分总和是否接近1.0
        double sum = 0.0;

        for (int i = 0; i < 32; i++) {
            sum += gasComponents[i];
        }

        if (sum < (100.0 - 1e-5) || sum > (100.0 + 1e-5)) {
            // 使用stringstream替代to_string，以兼容VS2010
            std::stringstream ss;
            ss << sum;
            throw runtime_error("气体组分总和应接近100.0%，当前值为: " + ss.str() + "%");
        }
    }

    wells = tempWells;
    limits = tempLimits;
    faults = tempFaults;
    fracs = tempFracs;
}

void PrintWell(const CBoundWell &w) {
    cout << "井类型: ";

    switch (w.wellType) {
        case VERTICAL:
            cout << "垂直井";
            break;

        case VERTICAL_FRACTURED:
            cout << "垂直压裂井";
            break;

        case MULTI_FRACED_HORIZONTAL:
            cout << "多段压裂水平井";
            break;
    }

    cout << "\n位置: (" << w.cCenter.x << ", " << w.cCenter.y << ")"
         << "\n半径: " << w.dRadius
         << "\n流动模式: " << (w.bisMultFlow ? "多级" : "单级")
         << "\n时间流量数据(" << w.nTimeNumQ << "个):\n";

    for (int i = 0; i < w.nTimeNumQ; ++i) {
        cout << "  " << w.pdTimeQ[i] << "hr: " << w.pdQ[i] << " m^3/d\n";
    }
}

// 气体组分名称数组(包含英文名称和化学式)
const char* gasComponentNames[32] = {
    "甲烷(Methane, CH4)", "乙烷(Ethane, C2H6)", "丙烷(Propane, C3H8)",
    "异丁烷(Isobutane, i-C4H10)", "正丁烷(n-Butane, n-C4H10)",
    "异戊烷(Isopentane, i-C5H12)", "正戊烷(n-Pentane, n-C5H12)",
    "己烷(Hexane, C6H14)", "庚烷(Heptane, C7H16)", "辛烷(Octane, C8H18)",
    "壬烷(Nonane, C9H20)", "葵烷(Decane, C10H22)",
    "十一烷(Undecane, C11H24)", "十二烷(Dodecane, C12H26)",
    "十三烷(Tridecane, C13H28)", "十四烷(Tetradecane, C14H30)",
    "十五烷(Pentadecane, C15H32)", "十六烷(Hexadecane, C16H34)",
    "十七烷(Heptadecane, C17H36)", "十八烷(Octadecane, C18H38)",
    "十九烷(Nonadecane, C19H40)", "廿烷(Icosane, C20H42)",
    "廿一烷(Henicosane, C21H44)", "廿一以上烷(C21+)",
    "硫化氢(Hydrogen Sulfide, H2S)", "二氧化碳(Carbon Dioxide, CO2)",
    "氮气(Nitrogen, N2)", "水(Water, H2O)", "氢气(Hydrogen, H2)",
    "氦气(Helium, He)", "氧气(Oxygen, O2)", "一氧化碳(Carbon Monoxide, CO)"
};

// 输出油藏基本信息及类型信息
void OutputReservoirInfo(const double pBaseData[8],
                         RESERVOIR_TYPE reservoirType,
                         double reservoirTemp,
                         const double gasComponents[32],
                         const std::string reservoirPara[8],
                         std::ostream &out = std::cout) {
    out << "====================== 油藏基本信息 ======================" << std::endl;

    // 使用预定义的参数名称
    for (int i = 0; i < 8; i++) {
        out << reservoirPara[i] << ": " << pBaseData[i];

        // 根据参数类型添加适当的单位
        if (i == 0) {
            out << " MPa";    // 初始压力
        } else if (i == 1) {
            out << " mD";    // 渗透率
        } else if (i == 2) {
            out << " m";    // 厚度
        } else if (i == 3) {
            out << " mPa·s";    // 油粘度
        } else if (i == 5) {
            out << " (分数)";    // 孔隙度
        } else if (i == 6) {
            out << " MPa^(-1)";    // 岩石压缩系数
        }

        out << std::endl;
    }

    // 输出油藏类型
    out << "油藏类型: ";

    if (reservoirType == OIL_RESERVOIR) {
        out << "油藏" << std::endl;
    } else if (reservoirType == GAS_RESERVOIR) {
        out << "气藏" << std::endl;
        out << "气藏温度: " << reservoirTemp << " K" << std::endl;

        // 输出气体组分
        out << "气体组分:" << std::endl;
        double sum = 0.0;

        for (int i = 0; i < 32; i++) {
            if (gasComponents[i] > 1.0e-5) {
                // 只输出非零组分
                // 使用名称数组，显示组分名称
                out << "  " << gasComponentNames[i] << ": " << gasComponents[i] << "%" << std::endl;
                sum += gasComponents[i];
            }
        }

        out << "组分总和: " << sum << "%" << std::endl;
    }

    out << "==========================================================" << std::endl;
}

DWORD WINAPI MessagePumpThread(LPVOID lpParam) {
    MSG msg;

    while (g_calculationRunning) {
        // 处理所有消息
        while (PeekMessage(&msg, NULL, 0, 0, PM_REMOVE)) {
            TranslateMessage(&msg);
            DispatchMessage(&msg);
        }

        // 显示进度条
        ShowProgressBar(g_timeStep, g_totalSteps);

        // 短暂休眠，减少CPU使用率
        Sleep(100);
    }

    return 0;
}

int main() {
    std::string reservoirPara[8] = {"初始压力", "渗透率", "厚度", "黏度", "体积系数", "孔隙度", "综合压缩系数", "渗透率各向异性系数"};

    try {
        CBoundShape shape;
        vector<CBoundWell> wells;
        vector<CBoundPolygon> limits;
        vector<Frac> faults, fracs;
        double pBaseData[8] = {0};
        TimeDis timeDis;
        RESERVOIR_TYPE reservoirType;
        double reservoirTemp;
        double gasComponents[32];

        ReadParameters("input/case_oil_faults.txt", shape, wells, limits, faults, fracs, pBaseData, timeDis, reservoirType,
                       reservoirTemp, gasComponents);

        // 详细油藏信息输出
        OutputReservoirInfo(pBaseData, reservoirType, reservoirTemp, gasComponents, reservoirPara);

        cout << "\n===== 井数据 =====" << endl;

        for (size_t i = 0; i < wells.size(); ++i) {
            cout << "\n井 " << i + 1 << ":\n";
            PrintWell(wells[i]);
        }

        /************************************************************************/
        /* load netgen mesh generation library and generate mesh                */
        /************************************************************************/
        HMODULE hMod = LoadLibrary("meshgene.dll");

        if (NULL == hMod) {
            std::cout << "meshgene.dll加载失败！\n";
            return 1;
        }

        typedef Ng_Mesh *(*NgMeshGen2DByIn2d)(std::vector<CBoundWell>& wells, const std::vector<CBoundPolygon> limits,
                                              const std::vector<Frac>& faults, const std::vector<Frac>& fracs, const CBoundShape & shape);
        typedef Ng_Mesh *(*NgMeshGen2D)(std::string strFileIn2D);
        typedef void *(*Ng2VTK)(Ng_Mesh * mesh);
        // Get the address of the Ng_Mesh function
        NgMeshGen2DByIn2d meshgene1 = (NgMeshGen2DByIn2d)GetProcAddress(hMod, "NgMeshGen2DByIn2d");
        NgMeshGen2D meshgene2 = (NgMeshGen2D)GetProcAddress(hMod, "NgMeshGen2D");

        if (NULL == meshgene1) {
            FreeLibrary(hMod);
            std::cout << "meshgene文件加载函数地址获取失败！\n";
            return 1;
        }

        Ng_Mesh *mesh = meshgene1(wells, limits, faults, fracs, shape);

        // Ng_Mesh *mesh=meshgene2("netgen-1.in2d");
        Ng2VTK ng = (Ng2VTK)GetProcAddress(hMod, "Ng2VTK");
        ng(mesh);
        // return 0;

        HMODULE hMod_solver = LoadLibrary("singlePhaseSolverDll.dll");

        if (NULL == hMod_solver) {
            std::cout << "singlePhaseSolverDll.dll加载失败！\n";
            return 1;
        }

        //      TimeDis timeDis;
        //      timeDis.dTB = 1.0e-5; //
        // 初始化全局进度变量
        g_timeStep = 0;
        g_totalSteps = 100;
        g_calculationRunning = true;

        // 创建消息处理线程
        std::cout << "开始计算，请稍候...\n";
        HANDLE hThread = CreateThread(NULL, 0, MessagePumpThread, NULL, 0, NULL);

        std::vector<std::pair<double, std::vector<double>>> vecFieldPre;
        int timeStep, totalSteps;

        if (reservoirType == OIL_RESERVOIR) {
            typedef int (*Solver)(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);
            Solver solverFun = (Solver)GetProcAddress(hMod_solver, "singlePhaseSolverNgmesh");

            if (NULL == solverFun) {
                FreeLibrary(hMod_solver);
                std::cout << "singlePhaseSolverNgmesh failed!\n";
                return 1;
            }

            //    std::vector<std::vector<Point>> vecBP;
            // std::vector<std::pair<double, std::vector<double>>> vecFieldPre;
            solverFun(timeDis, pBaseData, mesh, wells, limits, faults, fracs, shape, vecFieldPre, g_timeStep, g_totalSteps);
            // 更新全局变量，确保显示100%进度
        } else if (reservoirType == GAS_RESERVOIR) {
            typedef int (*Solver)(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);
            Solver solverFun = (Solver)GetProcAddress(hMod_solver, "singlePhaseGasSolverNgmesh");

            if (NULL == solverFun) {
                FreeLibrary(hMod_solver);
                std::cout << "singlePhaseGasSolverNgmesh failed!\n";
                return 1;
            }

            //    std::vector<std::vector<Point>> vecBP;
            solverFun(gasComponents, reservoirTemp, timeDis, pBaseData, mesh, wells, limits, faults, fracs, shape, vecFieldPre,
                      g_timeStep, g_totalSteps);
        }

        // 停止消息处理线程
        g_calculationRunning = false;

        if (hThread != NULL) {
            WaitForSingleObject(hThread, INFINITE);
            CloseHandle(hThread);
        }

        // 显示100%完成的进度条
        ShowProgressBar(g_totalSteps, g_totalSteps);
        std::cout << std::endl
                  << "计算完成！" << std::endl;

        // export bottom hole pressure history
        std::cout << "Export log-log plot and history plot data to files\n";
        std::ofstream historyFile;
        historyFile.precision(12);
        historyFile.open("output/history.txt");

        if (!historyFile.is_open()) {
            std::cerr << "Error opening file: " << "history.txt" << std::endl;
            return 1;
        }

        int nWell = 0;

        for (int i = 0; i < wells[nWell].m_PressureData.size(); i++) {
            historyFile << wells[nWell].m_PressureData[i].x << "\t" << wells[nWell].m_PressureData[i].y << "\t" << std::endl;
        }

        historyFile.close();

        // calculate the pressure derivative
        typedef bool (*PreLog)(const std::vector<Point>& pressure, const int& nSection, double* TimeQ, double* dQ, int nTimeQ,
                               std::vector<Point>& logPre);
        PreLog preLogFun = (PreLog)GetProcAddress(hMod_solver, "logLogPre");

        if (NULL == preLogFun) {
            FreeLibrary(hMod_solver);
            std::cout << "preLogFun failed!\n";
            return 1;
        }

        int nSection = 3;
        std::vector<Point> logPre;
        preLogFun(wells[nWell].m_PressureData, nSection, wells[nWell].pdTimeQ, wells[nWell].pdQ, wells[nWell].nTimeNumQ,
                  logPre);

        // export the pressure derivative
        std::ofstream logFile;
        logFile.precision(8);
        logFile.open("output/logPre.txt");

        if (!logFile.is_open()) {
            std::cerr << "Error opening file: " << "logPre.txt" << std::endl;
            return 1;
        }

        for (int i = 0; i < logPre.size(); i++) {
            logFile << logPre[i].x << "\t" << logPre[i].y << "\t" << logPre[i].z << "\t" << logPre[i].pointData[1] << std::endl;
        }

        logFile.close();

        // export the semi-log data
        std::ofstream semiLogFile;
        semiLogFile.precision(8);
        semiLogFile.open("output/semiLogPre.txt");

        if (!semiLogFile.is_open()) {
            std::cerr << "Error opening file: " << "semiLogPre.txt" << std::endl;
            return 1;
        }

        for (int i = 0; i < logPre.size(); i++) {
            semiLogFile << logPre[i].x << "\t" << logPre[i].pointData[0] << "\t" << std::endl;
        }

        semiLogFile.close();

        // export the rate of selected well
        std::ofstream rateFile;
        rateFile.precision(8);
        rateFile.open("output/rate.txt");

        if (!rateFile.is_open()) {
            std::cerr << "Error opening file: " << "rate.txt" << std::endl;
            return 1;
        }

        for (int i = 0; i < wells[nWell].nTimeNumQ; i++) {
            rateFile << wells[nWell].pdTimeQ[i] << "\t" << wells[nWell].pdQ[i] << std::endl;
        }

        rateFile.close();

        // export with pressure field with vtk format (version 4.2)
        for (int i = 0; i < vecFieldPre.size(); ++i) {
            if (i % 10 == 0) {
                Ng2VTKWithData(mesh, vecFieldPre[i].second, "output/" + std::to_string(long double(i)) + ".vtk");
            }
        }

        FreeLibrary(hMod);

        system("run_python.bat");
    } catch (const exception &e) {
        cerr << "错误: " << e.what() << endl;
        return 1;
    }

    return 0;
}