diff --git a/.gitignore b/.gitignore index 0a7e6d2..a93ac31 100644 --- a/.gitignore +++ b/.gitignore @@ -119,4 +119,5 @@ ML/nmWTAI-ML/data ML/nmWTAI-ML/models ML/nmWTAI-ML/results __pycache__ -ML/Training \ No newline at end of file +ML/Training/Debug +ML/Training/Release \ No newline at end of file diff --git a/ML/Training/Runner/RunnerIO.cpp b/ML/Training/Runner/RunnerIO.cpp new file mode 100644 index 0000000..268b24d --- /dev/null +++ b/ML/Training/Runner/RunnerIO.cpp @@ -0,0 +1,196 @@ +#include "RunnerIO.h" +#include +#include +#include +#include + +static uint32_t makeMagic(const char a, const char b, const char c, const char d) +{ + return (uint32_t)(uint8_t)a + | ((uint32_t)(uint8_t)b << 8) + | ((uint32_t)(uint8_t)c << 16) + | ((uint32_t)(uint8_t)d << 24); +} + +static bool readU32(std::ifstream& fs, uint32_t& v) +{ + fs.read((char*)&v, sizeof(v)); + return fs.good(); +} + +static bool writeU32(std::ofstream& fs, uint32_t v) +{ + fs.write((const char*)&v, sizeof(v)); + return fs.good(); +} + +static bool readDouble(std::ifstream& fs, double& v) +{ + fs.read((char*)&v, sizeof(v)); + return fs.good(); +} + +static inline bool writeDoubles(std::ofstream& fs, const double* p, size_t n) +{ + if (n == 0) return true; + fs.write((const char*)p, (std::streamsize)(n * sizeof(double))); + return fs.good(); +} + +static bool writeDouble(std::ofstream& fs, double v) +{ + fs.write((const char*)&v, sizeof(v)); + return fs.good(); +} + +namespace RunnerIO +{ + bool fileExists(const std::string& path) + { + std::ifstream fs(path.c_str(), std::ios::binary); + return fs.good(); + } + + void printParams(const RunnerParams& p) + { + std::cout << "Params:" + << " k=" << p.k + << ", skin=" << p.skin + << ", C=" << p.wellboreC + << ", phi=" << p.phi + << ", h=" << p.h + << ", Cf=" << p.Cf; + + if (!p.timeQ.empty() && !p.q.empty() && p.timeQ.size() == p.q.size()) { + std::cout << ", sectionIndex=" << p.sectionIndex + << ", nQ=" << p.timeQ.size(); + } else { + std::cout << ", schedule=DEFAULT(from scene)"; + } + std::cout << std::endl; + } + + bool readParamsBin(const std::string& path, RunnerParams& out) + { + std::ifstream fs(path.c_str(), std::ios::binary); + if (!fs) + { + std::cerr << "RunnerIO::readParamsBin: cannot open " << path << std::endl; + return false; + } + + const uint32_t MAGIC = makeMagic('P', 'R', 'M', '1'); + uint32_t magic = 0, ver = 0; + if (!readU32(fs, magic) || !readU32(fs, ver)) + { + std::cerr << "RunnerIO::readParamsBin: header read failed\n"; + return false; + } + if (magic != MAGIC || ver != 1) + { + std::cerr << "RunnerIO::readParamsBin: magic/version mismatch\n"; + return false; + } + + if (!readDouble(fs, out.k)) return false; + if (!readDouble(fs, out.skin)) return false; + if (!readDouble(fs, out.wellboreC)) return false; + if (!readDouble(fs, out.phi)) return false; + if (!readDouble(fs, out.h)) return false; + if (!readDouble(fs, out.Cf)) return false; + + out.sectionIndex = 0; + out.timeQ.clear(); + out.q.clear(); + + int nextByte = fs.peek(); + if (nextByte == EOF) { + return true; + } + + uint32_t sec = 0, nQ = 0; + if (!readU32(fs, sec)) { + std::cerr << "RunnerIO::readParamsBin: schedule extension present but missing sectionIndex\n"; + return false; + } + if (!readU32(fs, nQ)) { + std::cerr << "RunnerIO::readParamsBin: schedule extension present but missing nQ\n"; + return false; + } + if (nQ < 2 || nQ > 100000) { + std::cerr << "RunnerIO::readParamsBin: invalid nQ=" << nQ << "\n"; + return false; + } + + out.sectionIndex = sec; + out.timeQ.resize(nQ); + out.q.resize(nQ); + + for (uint32_t i = 0; i < nQ; ++i) { + if (!readDouble(fs, out.timeQ[i])) return false; + } + for (uint32_t i = 0; i < nQ; ++i) { + if (!readDouble(fs, out.q[i])) return false; + } + + return fs.good(); + } + + bool writeResultBin(const std::string& path, const RunnerResult& r) + { + std::ofstream fs(path.c_str(), std::ios::binary); + if (!fs) + { + std::cerr << "RunnerIO::writeResultBin: cannot open " << path << std::endl; + return false; + } + + //使用普通静态变量 + static std::vector buf; + buf.resize(4 * 1024 * 1024); + fs.rdbuf()->pubsetbuf(&buf[0], (std::streamsize)buf.size()); + + const uint32_t MAGIC = makeMagic('R', 'S', 'B', '1'); + if (!writeU32(fs, MAGIC)) return false; + if (!writeU32(fs, 1)) return false; + + if (!writeU32(fs, (uint32_t)r.nWells)) return false; + if (!writeU32(fs, (uint32_t)r.nSteps)) return false; + + // t + if (r.nSteps > 0) { + if (r.t.size() < r.nSteps) return false; + if (!writeDoubles(fs, &r.t[0], r.nSteps)) return false; // 使用 &r.t[0] 替代 .data() + } + + // pw + for (unsigned int w = 0; w < r.nWells; ++w) + { + if (r.nSteps > 0) { + if (w >= r.pw.size()) return false; + if (r.pw[w].size() < r.nSteps) return false; + if (!writeDoubles(fs, &r.pw[w][0], r.nSteps)) return false; // 使用 &r.pw[w][0] + } + } + + // loglog + for (unsigned int w = 0; w < r.nWells; ++w) + { + uint32_t nLogLog = 0; + if (w < r.loglog_t.size()) nLogLog = (uint32_t)r.loglog_t[w].size(); + if (!writeU32(fs, nLogLog)) return false; + + if (nLogLog == 0) continue; + + if (w >= r.loglog_p.size() || w >= r.loglog_deriv.size()) return false; + if (r.loglog_p[w].size() != nLogLog) return false; + if (r.loglog_deriv[w].size() != nLogLog) return false; + + if (!writeDoubles(fs, &r.loglog_t[w][0], nLogLog)) return false; + if (!writeDoubles(fs, &r.loglog_p[w][0], nLogLog)) return false; + if (!writeDoubles(fs, &r.loglog_deriv[w][0], nLogLog)) return false; + } + + return fs.good(); + } +} \ No newline at end of file diff --git a/ML/Training/Runner/RunnerIO.h b/ML/Training/Runner/RunnerIO.h new file mode 100644 index 0000000..8c84bdb --- /dev/null +++ b/ML/Training/Runner/RunnerIO.h @@ -0,0 +1,48 @@ +#pragma once +#include +#include +#include + +struct RunnerParams +{ + // 6 params + double k; + double skin; + double wellboreC; + double phi; + double h; + double Cf; + + // optional schedule extension + uint32_t sectionIndex; // wellFlowSectionIndex + std::vector timeQ; // hours + std::vector q; // m^3/d + + RunnerParams() + : k(0), skin(0), wellboreC(0), phi(0), h(0), Cf(0), sectionIndex(0) + {} +}; + +struct RunnerResult +{ + unsigned int nWells; + unsigned int nSteps; + std::vector t; + std::vector< std::vector > pw; + + std::vector< std::vector > loglog_t; + std::vector< std::vector > loglog_p; + std::vector< std::vector > loglog_deriv; + + RunnerResult() : nWells(0), nSteps(0) {} +}; + +namespace RunnerIO +{ + bool fileExists(const std::string& path); + + bool readParamsBin(const std::string& path, RunnerParams& out); + bool writeResultBin(const std::string& path, const RunnerResult& r); + + void printParams(const RunnerParams& p); +} diff --git a/ML/Training/Runner/runner_main.cpp b/ML/Training/Runner/runner_main.cpp new file mode 100644 index 0000000..c29f7e2 --- /dev/null +++ b/ML/Training/Runner/runner_main.cpp @@ -0,0 +1,444 @@ +#include "singlePhaseSolver.h" +#include +#include +#include +#include +#include +#include +#include + +#include "pch.h" +#include "DatasetIO.h" +#include "RunnerIO.h" + +// HX_NWTM +typedef void (*HX_NWTM_MODEL_Func)(HX_NWTM_MODEL_OUTPUT&, const HX_NWTM_MODEL_INPUT&, std::string); + +// log-log pre +typedef bool (*PreLogFun)(const std::vector&, const int&, double*, double*, int, std::vector&); + +static bool fileExistsA(const std::string& p) +{ + DWORD attr = GetFileAttributesA(p.c_str()); + return (attr != INVALID_FILE_ATTRIBUTES) && !(attr & FILE_ATTRIBUTE_DIRECTORY); +} + +static std::string getExeDir() +{ + char buf[MAX_PATH] = { 0 }; + DWORD len = GetModuleFileNameA(NULL, buf, MAX_PATH); + if (len == 0) return "."; + std::string path(buf, len); + size_t pos = path.find_last_of("\\/"); + if (pos != std::string::npos) path.resize(pos); + return path; +} + +// =============== 组装模型输入 =============== +static HX_NWTM_MODEL_INPUT buildModelInputFromDataset(const PebiScene& scene, const HX_NWTM_GRID_OUTPUT2& gridOutput2) +{ + HX_NWTM_MODEL_INPUT input(gridOutput2); + input.T = scene.solverType; + + input.Rate.t = scene.Rate.t; + input.Rate.qo = scene.Rate.qo; + input.Rate.qg = scene.Rate.qg; + input.Rate.qw = scene.Rate.qw; + + input.CS.C = scene.CS.C; + input.CS.S = scene.CS.S; + + input.PVT.p = scene.PVT.p; + input.PVT.pb = scene.PVT.pb; + input.PVT.Rso = scene.PVT.Rso; + input.PVT.Bo = scene.PVT.Bo; + input.PVT.Co = scene.PVT.Co; + input.PVT.miuo = scene.PVT.miuo; + input.PVT.rouo = scene.PVT.rouo; + input.PVT.Rv = scene.PVT.Rv; + input.PVT.Bg = scene.PVT.Bg; + input.PVT.Cg = scene.PVT.Cg; + input.PVT.miug = scene.PVT.miug; + input.PVT.roug = scene.PVT.roug; + input.PVT.Z = scene.PVT.Z; + input.PVT.Rsw = scene.PVT.Rsw; + input.PVT.Bw = scene.PVT.Bw; + input.PVT.Cw = scene.PVT.Cw; + input.PVT.miuw = scene.PVT.miuw; + input.PVT.rouw = scene.PVT.rouw; + input.PVT.V = scene.PVT.V; + input.PVT.k_kinitial = scene.PVT.k_kinitial; + input.PVT.Cf_Cfinitial = scene.PVT.Cf_Cfinitial; + input.PVT.So = scene.PVT.So; + input.PVT.Kro = scene.PVT.Kro; + input.PVT.Sg = scene.PVT.Sg; + input.PVT.Krg = scene.PVT.Krg; + input.PVT.Sw = scene.PVT.Sw; + input.PVT.Krw = scene.PVT.Krw; + + input.Base.Pi = scene.Base.Pi; + input.Base.Cti = scene.Base.Cti; + input.Base.Cf = scene.Base.Cf; + input.Base.Soi = scene.Base.Soi; + input.Base.Sgi = scene.Base.Sgi; + input.Base.Swi = scene.Base.Swi; + input.Base.d = scene.Base.d; + input.Base.dt_Min = scene.Base.dt_Min; + input.Base.dt_Max = scene.Base.dt_Max; + + size_t nCells = gridOutput2.Trinodexy.size(); + input.Base.k = dVec1(nCells, scene.Base.k_ref); + input.Base.phi = dVec1(nCells, scene.Base.phi_ref); + input.Base.h = dVec1(nCells, scene.Base.h_ref); + + return input; +} + +static inline size_t inferWellCount(const HX_NWTM_MODEL_INPUT& in) +{ + size_t n = 0; + n = (std::max)(n, in.Rate.t.size()); + n = (std::max)(n, in.CS.C.size()); + n = (std::max)(n, in.CS.S.size()); + if (n == 0) n = 1; + return n; +} + +// =============== 应用采样参数 =============== +static void applySampledParamsAndMaybeOverrideRate(HX_NWTM_MODEL_INPUT& in, const RunnerParams& p) +{ + const size_t nWells = inferWellCount(in); + + // Ensure well vector sizes (不会清空 capacity) + if (in.Rate.t.size() < nWells) in.Rate.t.resize(nWells); + if (in.Rate.qo.size() < nWells) in.Rate.qo.resize(nWells); + if (in.Rate.qg.size() < nWells) in.Rate.qg.resize(nWells); + if (in.Rate.qw.size() < nWells) in.Rate.qw.resize(nWells); + + if (in.CS.S.size() < nWells) in.CS.S.resize(nWells, 0.0); + if (in.CS.C.size() < nWells) in.CS.C.resize(nWells, 0.0); + + // 覆盖 skin / wellboreC（fill 更快） + std::fill(in.CS.S.begin(), in.CS.S.end(), p.skin); + std::fill(in.CS.C.begin(), in.CS.C.end(), p.wellboreC); + + // 覆盖 k/phi/h/Cf（resize + fill，复用 capacity） + const size_t nCells = in.GRID.Trinodexy.size(); + in.Base.k.resize(nCells); + in.Base.phi.resize(nCells); + in.Base.h.resize(nCells); + std::fill(in.Base.k.begin(), in.Base.k.end(), p.k); + std::fill(in.Base.phi.begin(), in.Base.phi.end(), p.phi); + std::fill(in.Base.h.begin(), in.Base.h.end(), p.h); + in.Base.Cf = p.Cf; + + // 制度覆盖：params.bin 带 schedule 才覆盖 + if (!p.timeQ.empty() && !p.q.empty() && p.timeQ.size() == p.q.size()) + { + const size_t nQ = p.timeQ.size(); + for (size_t w = 0; w < nWells; ++w) + { + dVec1& tt = in.Rate.t[w]; + dVec1& qo = in.Rate.qo[w]; + dVec1& qg = in.Rate.qg[w]; + dVec1& qw = in.Rate.qw[w]; + + tt.resize(nQ); + qo.resize(nQ); + qg.resize(nQ); + qw.resize(nQ); + + std::copy(p.timeQ.begin(), p.timeQ.end(), tt.begin()); + std::copy(p.q.begin(), p.q.end(), qo.begin()); + std::fill(qg.begin(), qg.end(), 0.0); + std::fill(qw.begin(), qw.end(), 0.0); + } + } + +} + +// =============== 转换求解结果=============== +static void toRunnerResult(const HX_NWTM_MODEL_OUTPUT& out, RunnerResult& rr) +{ + rr.nSteps = (unsigned int)out.t.size(); + rr.t = out.t; + + rr.nWells = (unsigned int)out.pw.size(); + rr.pw.resize(rr.nWells); + + for (unsigned int w = 0; w < rr.nWells; ++w) { + rr.pw[w] = out.pw[w]; + if (rr.pw[w].size() < rr.nSteps) rr.pw[w].resize(rr.nSteps, 0.0); + } +} + +// =============== 计算双对数曲线（传入已加载的 preLogFun） =============== +static bool computeLogLogCurves(RunnerResult& rr, const PebiScene& scene, const RunnerParams& params, PreLogFun preLogFun) +{ + + // 无论有没有 DLL，都先保证结构完整，避免 writeResultBin 越界/串数据 + rr.loglog_t.assign(rr.nWells, std::vector()); + rr.loglog_p.assign(rr.nWells, std::vector()); + rr.loglog_deriv.assign(rr.nWells, std::vector()); + + if (!preLogFun) { + return true; // 没 dll：loglog 为空，但结构正确 + } + + for (unsigned int w = 0; w < rr.nWells; ++w) + { + std::vector wellPressureData; + wellPressureData.resize(rr.nSteps); + for (unsigned int i = 0; i < rr.nSteps; ++i) { + wellPressureData[i].x = rr.t[i]; + wellPressureData[i].y = rr.pw[w][i]; + wellPressureData[i].z = 0.0; + } + + int sectionIndex = 0; + std::vector timeQ; + std::vector q; + + // 1) params schedule 优先 + if (!params.timeQ.empty() && !params.q.empty() && params.timeQ.size() == params.q.size()) + { + sectionIndex = (int)params.sectionIndex; + timeQ = params.timeQ; + q = params.q; + } + else + { + // 2) scene 默认 + if (!scene.wellFlowSectionIndex.empty()) { + size_t idx = (w < (unsigned int)scene.wellFlowSectionIndex.size()) ? w : 0; + sectionIndex = scene.wellFlowSectionIndex[idx]; + } + + if (w < scene.Rate.t.size()) timeQ = scene.Rate.t[w]; + if (w < scene.Rate.qo.size()) q = scene.Rate.qo[w]; + } + + if (timeQ.size() < 2 || q.size() != timeQ.size()) { + return false; + } + + std::vector logPreResult; + bool success = preLogFun( + wellPressureData, + sectionIndex, + timeQ.empty() ? NULL : &timeQ[0], // 安全处理空vector + q.empty() ? NULL : &q[0], + (int)timeQ.size(), + logPreResult + ); + + if (!success || logPreResult.size() <= 2) continue; + + std::vector& lt = rr.loglog_t[w]; + std::vector& lp = rr.loglog_p[w]; + std::vector& ld = rr.loglog_deriv[w]; + + lt.reserve(logPreResult.size() - 2); + lp.reserve(logPreResult.size() - 2); + ld.reserve(logPreResult.size() - 2); + + for (size_t i = 1; i + 1 < logPreResult.size(); ++i) { + const Point& pt = logPreResult[i]; + lt.push_back(pt.x); + lp.push_back(pt.y); + ld.push_back(pt.z); + } + } + + return true; +} + +// =============== server 模式：循环处理 paramsPath/resultPath =============== +static int runServer(const std::string& datasetPath, + const std::string& dllPath, + const std::string& licPath) +{ + if (!fileExistsA(datasetPath)) { std::cerr << "ERROR: dataset not found\n"; return 10; } + if (!fileExistsA(dllPath)) { std::cerr << "ERROR: dll not found\n"; return 12; } + if (!fileExistsA(licPath)) { std::cerr << "ERROR: license not found\n"; return 13; } + + // 1) load dataset once + PebiScene scene; + HX_NWTM_GRID_OUTPUT2 gridOutput2; + if (!DatasetIO::loadDataset(datasetPath, scene, gridOutput2, true)) { + std::cerr << "ERROR: loadDataset failed\n"; + return 30; + } + + // 2) prepare base input once + HX_NWTM_MODEL_INPUT modelInput = buildModelInputFromDataset(scene, gridOutput2); + + // 3) load HX solver once + HMODULE hx = LoadLibraryA(dllPath.c_str()); + if (!hx) { std::cerr << "ERROR: LoadLibrary(HX) failed\n"; return 40; } + + HX_NWTM_MODEL_Func HX_NWTM_MODEL = (HX_NWTM_MODEL_Func)GetProcAddress(hx, "HX_NWTM_MODEL"); + if (!HX_NWTM_MODEL) { std::cerr << "ERROR: GetProcAddress(HX_NWTM_MODEL) failed\n"; FreeLibrary(hx); return 41; } + + // 4) load loglog dll once + PreLogFun preLogFun = NULL; + HMODULE hMod_solver = LoadLibrary(L"singlePhaseSolverDll.dll"); + if (hMod_solver) { + preLogFun = (PreLogFun)GetProcAddress(hMod_solver, "logLogPre"); + if (!preLogFun) { + FreeLibrary(hMod_solver); + hMod_solver = NULL; + } + } + + std::string line; + RunnerResult rr; + + // 输入格式：每行 "paramsPath resultPath" + while (std::getline(std::cin, line)) + { + if (line.empty()) continue; + std::istringstream iss(line); + std::string paramsPath, resultPath; + if (!(iss >> paramsPath >> resultPath)) { + std::cerr << "ERROR: invalid input line, expected: paramsPath resultPath\n"; + continue; + } + + RunnerParams params; + if (!RunnerIO::readParamsBin(paramsPath, params)) { + std::cerr << "ERROR: read params failed: " << paramsPath << "\n"; + continue; + } + + // 覆盖参数（复用内存） + applySampledParamsAndMaybeOverrideRate(modelInput, params); + + HX_NWTM_MODEL_OUTPUT modelOutput; + try { + HX_NWTM_MODEL(modelOutput, modelInput, licPath); + } + catch (...) { + std::cerr << "ERROR: solver exception\n"; + continue; + } + + toRunnerResult(modelOutput, rr); + + if (!computeLogLogCurves(rr, scene, params, preLogFun)) { + std::cerr << "ERROR: computeLogLogCurves failed\n"; + continue; + } + + if (!RunnerIO::writeResultBin(resultPath, rr)) { + std::cerr << "ERROR: write result failed: " << resultPath << "\n"; + continue; + } + + // 方便 Python 端读取：成功就输出一行 OK + std::cout << "OK " << resultPath << "\n"; + std::cout.flush(); + } + + if (hMod_solver) FreeLibrary(hMod_solver); + FreeLibrary(hx); + return 0; +} + +// =============== 单次模式 =============== +int main(int argc, char** argv) +{ + std::string exeDir = getExeDir(); + + // server: runner.exe --server dataset.bin HX_NWTM.dll lic.dat + if (argc >= 2 && std::string(argv[1]) == "--server") + { + if (argc < 5) { + std::cerr << "Usage: runner.exe --server \n"; + return 2; + } + return runServer(argv[2], argv[3], argv[4]); + } + + // 默认单次模式 + std::string dataDir = exeDir + "\\..\\..\\nmWTAI-ML\\data\\temp"; + + std::string datasetPath = dataDir + "\\dataset.bin"; + std::string paramsPath = dataDir + "\\params.bin"; + std::string resultPath = dataDir + "\\result.bin"; + std::string dllPath = exeDir + "\\HX_NWTM.dll"; + std::string licPath = exeDir + "\\..\\..\\..\\Bin\\Res\\license\\HXNWTM_license.dat"; + + if (argc >= 2) datasetPath = argv[1]; + if (argc >= 3) paramsPath = argv[2]; + if (argc >= 4) resultPath = argv[3]; + if (argc >= 5) dllPath = argv[4]; + if (argc >= 6) licPath = argv[5]; + + if (!fileExistsA(datasetPath)) { std::cerr << "ERROR: dataset not found\n"; return 10; } + if (!fileExistsA(paramsPath)) { std::cerr << "ERROR: params not found\n"; return 11; } + if (!fileExistsA(dllPath)) { std::cerr << "ERROR: dll not found\n"; return 12; } + if (!fileExistsA(licPath)) { std::cerr << "ERROR: license not found\n"; return 13; } + + RunnerParams params; + if (!RunnerIO::readParamsBin(paramsPath, params)) { + std::cerr << "ERROR: read params failed\n"; + return 20; + } + + PebiScene scene; + HX_NWTM_GRID_OUTPUT2 gridOutput2; + if (!DatasetIO::loadDataset(datasetPath, scene, gridOutput2, true)) { + std::cerr << "ERROR: loadDataset failed\n"; + return 30; + } + + HX_NWTM_MODEL_INPUT modelInput = buildModelInputFromDataset(scene, gridOutput2); + applySampledParamsAndMaybeOverrideRate(modelInput, params); + + HMODULE hx = LoadLibraryA(dllPath.c_str()); + if (!hx) { std::cerr << "ERROR: LoadLibrary failed\n"; return 40; } + + HX_NWTM_MODEL_Func HX_NWTM_MODEL = (HX_NWTM_MODEL_Func)GetProcAddress(hx, "HX_NWTM_MODEL"); + if (!HX_NWTM_MODEL) { std::cerr << "ERROR: GetProcAddress failed\n"; FreeLibrary(hx); return 41; } + + // loglog dll + PreLogFun preLogFun = NULL; + HMODULE hMod_solver = LoadLibrary(L"singlePhaseSolverDll.dll"); + if (hMod_solver) { + preLogFun = (PreLogFun)GetProcAddress(hMod_solver, "logLogPre"); + if (!preLogFun) { FreeLibrary(hMod_solver); hMod_solver = NULL; } + } + + HX_NWTM_MODEL_OUTPUT modelOutput; + try { + HX_NWTM_MODEL(modelOutput, modelInput, licPath); + } + catch (...) { + std::cerr << "ERROR: solver exception\n"; + if (hMod_solver) FreeLibrary(hMod_solver); + FreeLibrary(hx); + return 50; + } + + RunnerResult rr; + toRunnerResult(modelOutput, rr); + + if (!computeLogLogCurves(rr, scene, params, preLogFun)) { + std::cerr << "ERROR: computeLogLogCurves failed\n"; + if (hMod_solver) FreeLibrary(hMod_solver); + FreeLibrary(hx); + return 55; + } + + if (!RunnerIO::writeResultBin(resultPath, rr)) { + std::cerr << "ERROR: write result failed\n"; + if (hMod_solver) FreeLibrary(hMod_solver); + FreeLibrary(hx); + return 60; + } + + if (hMod_solver) FreeLibrary(hMod_solver); + FreeLibrary(hx); + return 0; +} \ No newline at end of file diff --git a/ML/Training/Training/DatasetIO.cpp b/ML/Training/Training/DatasetIO.cpp new file mode 100644 index 0000000..7e18aeb --- /dev/null +++ b/ML/Training/Training/DatasetIO.cpp @@ -0,0 +1,780 @@ +#include "DatasetIO.h" + +// ---------------- public ---------------- + +bool DatasetIO::fileExists(const std::string& filename) +{ + std::ifstream fs(filename.c_str(), std::ios::binary); + return fs.good(); +} + +bool DatasetIO::saveDataset(const std::string& filename, + const PebiScene& scene, + const HX_NWTM_GRID_INPUT& gridInput, + const HX_NWTM_GRID_OUTPUT2& grid) +{ + std::ofstream fs(filename.c_str(), std::ios::binary); + if (!fs) { + std::cerr << "DatasetIO::saveDataset: 无法打开文件: " << filename << std::endl; + return false; + } + + DatasetHeader hdr; + hdr.magic = (unsigned int)MAGIC; + hdr.version = (unsigned int)VERSION; + hdr.reserved = 0; + hdr.sceneKey = computeSceneKey64(gridInput); + hdr.nCells = (unsigned int)grid.Trinodexy.size(); + + // 井数：优先用 Rate.t.size()，没有就用 wellName.size() + unsigned int wells = 0; + if (!scene.Rate.t.empty()) wells = (unsigned int)scene.Rate.t.size(); + else wells = (unsigned int)scene.wellName.size(); + + hdr.nWells = wells; + hdr.solverType = scene.solverType; + + if (!writeHeader(fs, hdr)) { + std::cerr << "DatasetIO::saveDataset: 写 header 失败\n"; + return false; + } + + // 写 scene（为了 dataset 自包含） + if (!writeScene(fs, scene)) { + std::cerr << "DatasetIO::saveDataset: 写 scene 失败\n"; + return false; + } + + // 写 gridOutput2 + if (!writeGridOutput2(fs, grid)) { + std::cerr << "DatasetIO::saveDataset: 写 gridOutput2 失败\n"; + return false; + } + + return fs.good(); +} + +bool DatasetIO::readDatasetSceneKey(const std::string& filename, + unsigned long long& outSceneKey, + unsigned int& outCells, + unsigned int& outWells, + int& outSolverType) +{ + std::ifstream fs(filename.c_str(), std::ios::binary); + if (!fs) return false; + + DatasetHeader hdr; + if (!readHeader(fs, hdr)) return false; + if (hdr.magic != (unsigned int)MAGIC) return false; + if (hdr.version != (unsigned int)VERSION) return false; + + outSceneKey = hdr.sceneKey; + outCells = hdr.nCells; + outWells = hdr.nWells; + outSolverType = hdr.solverType; + return true; +} + +bool DatasetIO::loadDataset(const std::string& filename, + PebiScene& outScene, + HX_NWTM_GRID_OUTPUT2& outGrid, + bool strictVerify) +{ + std::ifstream fs(filename.c_str(), std::ios::binary); + if (!fs) { + std::cerr << "DatasetIO::loadDataset: 无法打开文件: " << filename << std::endl; + return false; + } + + DatasetHeader hdr; + if (!readHeader(fs, hdr)) { + std::cerr << "DatasetIO::loadDataset: 读 header 失败\n"; + return false; + } + if (hdr.magic != (unsigned int)MAGIC) { + std::cerr << "DatasetIO::loadDataset: magic 不匹配\n"; + return false; + } + if (hdr.version != (unsigned int)VERSION) { + std::cerr << "DatasetIO::loadDataset: version 不匹配\n"; + return false; + } + + if (!readScene(fs, outScene)) { + std::cerr << "DatasetIO::loadDataset: 读 scene 失败\n"; + return false; + } + + if (!readGridOutput2(fs, outGrid)) { + std::cerr << "DatasetIO::loadDataset: 读 gridOutput2 失败\n"; + return false; + } + + // 快速检查 cell 数 + if (hdr.nCells != (unsigned int)outGrid.Trinodexy.size()) { + std::cerr << "DatasetIO::loadDataset: nCells 不匹配，dataset 可能损坏\n"; + return false; + } + + // 严格校验：用 scene 重建 gridInput 再算一次 sceneKey，对比 header，防止写错/损坏 + if (strictVerify) { + HX_NWTM_GRID_INPUT in = buildGridInputFromScene(outScene); + unsigned long long key = computeSceneKey64(in); + if (key != hdr.sceneKey) { + std::cerr << "DatasetIO::loadDataset: sceneKey 校验失败（dataset 内容与 header 不一致）\n"; + return false; + } + } + + return fs.good(); +} + +// ---------------- header io ---------------- + +bool DatasetIO::writeHeader(std::ofstream& fs, const DatasetHeader& hdr) +{ + if (!writeUInt(fs, hdr.magic)) return false; + if (!writeUInt(fs, hdr.version)) return false; + if (!writeUInt(fs, hdr.reserved)) return false; + if (!writeULL(fs, hdr.sceneKey)) return false; + if (!writeUInt(fs, hdr.nCells)) return false; + if (!writeUInt(fs, hdr.nWells)) return false; + if (!writeInt(fs, hdr.solverType)) return false; + return true; +} + +bool DatasetIO::readHeader(std::ifstream& fs, DatasetHeader& hdr) +{ + if (!readUInt(fs, hdr.magic)) return false; + if (!readUInt(fs, hdr.version)) return false; + if (!readUInt(fs, hdr.reserved)) return false; + if (!readULL(fs, hdr.sceneKey)) return false; + if (!readUInt(fs, hdr.nCells)) return false; + if (!readUInt(fs, hdr.nWells)) return false; + if (!readInt(fs, hdr.solverType)) return false; + return true; +} + +// ---------------- scene -> gridInput ---------------- + +HX_NWTM_GRID_INPUT DatasetIO::buildGridInputFromScene(const PebiScene& scene) +{ + HX_NWTM_GRID_INPUT input; + input.D = scene.D; + input.GridControl = scene.GridControl; + input.Boundary = scene.Boundary; + input.VerticalWell = scene.VerticalWell; + input.HorizontalWell = scene.HorizontalWell; + input.FractureVerticalWell = scene.FractureVerticalWell; + input.MultistageFracturedHorizontalWell = scene.MultistageFracturedHorizontalWell; + input.InclinedWell = scene.InclinedWell; + input.Fault = scene.Fault; + return input; +} + +// ---------------- primitive io ---------------- + +bool DatasetIO::writeUInt(std::ofstream& fs, unsigned int v) { fs.write((const char*)&v, sizeof(v)); return fs.good(); } +bool DatasetIO::writeULL(std::ofstream& fs, unsigned long long v) { fs.write((const char*)&v, sizeof(v)); return fs.good(); } +bool DatasetIO::writeInt(std::ofstream& fs, int v) { fs.write((const char*)&v, sizeof(v)); return fs.good(); } +bool DatasetIO::writeDouble(std::ofstream& fs, double v) { fs.write((const char*)&v, sizeof(v)); return fs.good(); } + +bool DatasetIO::readUInt(std::ifstream& fs, unsigned int& v) { fs.read((char*)&v, sizeof(v)); return fs.good(); } +bool DatasetIO::readULL(std::ifstream& fs, unsigned long long& v) { fs.read((char*)&v, sizeof(v)); return fs.good(); } +bool DatasetIO::readInt(std::ifstream& fs, int& v) { fs.read((char*)&v, sizeof(v)); return fs.good(); } +bool DatasetIO::readDouble(std::ifstream& fs, double& v) { fs.read((char*)&v, sizeof(v)); return fs.good(); } + +bool DatasetIO::writeSafeSize(std::ofstream& fs, size_t sz) +{ + unsigned int u = (unsigned int)sz; + return writeUInt(fs, u); +} + +bool DatasetIO::readSafeSize(std::ifstream& fs, size_t& sz) +{ + unsigned int u = 0; + if (!readUInt(fs, u)) return false; + + // 防御：避免异常超大导致内存爆炸 + const unsigned int MAX_SIZE = 100000000U; + if (u > MAX_SIZE) return false; + + sz = (size_t)u; + return true; +} + +// ---------------- string / std::vector IO ---------------- +// 说明：dataset.bin 自己定义字符串格式：uint32 字节长度 + raw bytes（UTF-8/ASCII） +// 0xFFFFFFFF 表示 null 字符串 + +bool DatasetIO::writeString(std::ofstream& fs, const std::string& s) +{ + unsigned int len = (unsigned int)s.size(); + if (!writeUInt(fs, len)) return false; + if (len > 0) { + fs.write(s.data(), len); + if (!fs.good()) return false; + } + return true; +} + +bool DatasetIO::readString(std::ifstream& fs, std::string& s) +{ + unsigned int len = 0; + if (!readUInt(fs, len)) return false; + + if (len == 0xFFFFFFFFU) { s.clear(); return true; } + + // 防御 + const unsigned int MAX_STR = 16 * 1024 * 1024U; + if (len > MAX_STR) return false; + + s.assign(len, '\0'); + if (len > 0) { + fs.read(&s[0], len); + if (!fs.good()) return false; + } + return true; +} + +bool DatasetIO::writeStdVecD(std::ofstream& fs, const std::vector& v) +{ + if (!writeSafeSize(fs, v.size())) return false; + if (!v.empty()) { + fs.write((const char*)&v[0], (std::streamsize)(v.size() * sizeof(double))); + if (!fs.good()) return false; + } + return true; +} + +bool DatasetIO::readStdVecD(std::ifstream& fs, std::vector& v) +{ + size_t sz = 0; + if (!readSafeSize(fs, sz)) return false; + v.resize(sz); + if (sz > 0) { + fs.read((char*)&v[0], (std::streamsize)(sz * sizeof(double))); + if (!fs.good()) return false; + } + return true; +} + +bool DatasetIO::writeStdVecI(std::ofstream& fs, const std::vector& v) +{ + if (!writeSafeSize(fs, v.size())) return false; + if (!v.empty()) { + fs.write((const char*)&v[0], (std::streamsize)(v.size() * sizeof(int))); + if (!fs.good()) return false; + } + return true; +} + +bool DatasetIO::readStdVecI(std::ifstream& fs, std::vector& v) +{ + size_t sz = 0; + if (!readSafeSize(fs, sz)) return false; + v.resize(sz); + if (sz > 0) { + fs.read((char*)&v[0], (std::streamsize)(sz * sizeof(int))); + if (!fs.good()) return false; + } + return true; +} + +bool DatasetIO::writeStdVec2D(std::ofstream& fs, const std::vector >& v) +{ + if (!writeSafeSize(fs, v.size())) return false; + for (size_t i = 0; i < v.size(); ++i) { + if (!writeStdVecD(fs, v[i])) return false; + } + return true; +} + +bool DatasetIO::readStdVec2D(std::ifstream& fs, std::vector >& v) +{ + size_t rows = 0; + if (!readSafeSize(fs, rows)) return false; + v.resize(rows); + for (size_t i = 0; i < rows; ++i) { + if (!readStdVecD(fs, v[i])) return false; + } + return true; +} + +bool DatasetIO::writeStdVec3D(std::ofstream& fs, const std::vector > >& v) +{ + if (!writeSafeSize(fs, v.size())) return false; + for (size_t i = 0; i < v.size(); ++i) { + if (!writeStdVec2D(fs, v[i])) return false; + } + return true; +} + +bool DatasetIO::readStdVec3D(std::ifstream& fs, std::vector > >& v) +{ + size_t depth = 0; + if (!readSafeSize(fs, depth)) return false; + v.resize(depth); + for (size_t i = 0; i < depth; ++i) { + if (!readStdVec2D(fs, v[i])) return false; + } + return true; +} + +// ---------------- scene IO ---------------- + +bool DatasetIO::writeScene(std::ofstream& fs, const PebiScene& s) +{ + // 基本版本 + if (!writeInt(fs, s.version)) return false; + + // 网格几何 + if (!writeInt(fs, s.D)) return false; + if (!writeDouble(fs, s.GridControl)) return false; + if (!writeStdVec2D(fs, s.Boundary)) return false; + if (!writeStdVec2D(fs, s.VerticalWell)) return false; + if (!writeStdVec2D(fs, s.HorizontalWell)) return false; + if (!writeStdVec2D(fs, s.FractureVerticalWell)) return false; + if (!writeStdVec3D(fs, s.MultistageFracturedHorizontalWell)) return false; + if (!writeStdVec2D(fs, s.InclinedWell)) return false; + if (!writeStdVec2D(fs, s.Fault)) return false; + + // 井顺序 + if (!writeStdVecI(fs, s.wellType)) return false; + + // 井名 + if (!writeSafeSize(fs, s.wellName.size())) return false; + for (size_t i = 0; i < s.wellName.size(); ++i) { + if (!writeString(fs, s.wellName[i])) return false; + } + + // 求解器 + if (!writeInt(fs, s.solverType)) return false; + + // Rate + // 注意：scene.Rate.t 是 vector> + // 为了兼容，这里存原始数据，不做转换 + if (!writeStdVec2D(fs, s.Rate.t)) return false; + if (!writeStdVec2D(fs, s.Rate.qo)) return false; + if (!writeStdVec2D(fs, s.Rate.qg)) return false; + if (!writeStdVec2D(fs, s.Rate.qw)) return false; + + // CS + if (!writeStdVecD(fs, s.CS.C)) return false; + if (!writeStdVecD(fs, s.CS.S)) return false; + + // 流动段 + if (!writeStdVecI(fs, s.wellFlowSectionIndex)) return false; + + // PVT + if (!writeStdVecD(fs, s.PVT.p)) return false; + if (!writeDouble(fs, s.PVT.pb)) return false; + + if (!writeStdVecD(fs, s.PVT.Rso)) return false; + if (!writeStdVecD(fs, s.PVT.Bo)) return false; + if (!writeStdVecD(fs, s.PVT.Co)) return false; + if (!writeStdVecD(fs, s.PVT.miuo)) return false; + if (!writeStdVecD(fs, s.PVT.rouo)) return false; + + if (!writeStdVecD(fs, s.PVT.Rv)) return false; + if (!writeStdVecD(fs, s.PVT.Bg)) return false; + if (!writeStdVecD(fs, s.PVT.Cg)) return false; + if (!writeStdVecD(fs, s.PVT.miug)) return false; + if (!writeStdVecD(fs, s.PVT.roug)) return false; + if (!writeStdVecD(fs, s.PVT.Z)) return false; + + if (!writeStdVecD(fs, s.PVT.Rsw)) return false; + if (!writeStdVecD(fs, s.PVT.Bw)) return false; + if (!writeStdVecD(fs, s.PVT.Cw)) return false; + if (!writeStdVecD(fs, s.PVT.miuw)) return false; + if (!writeStdVecD(fs, s.PVT.rouw)) return false; + + if (!writeStdVecD(fs, s.PVT.V)) return false; + if (!writeStdVecD(fs, s.PVT.k_kinitial)) return false; + if (!writeStdVecD(fs, s.PVT.Cf_Cfinitial)) return false; + + if (!writeStdVecD(fs, s.PVT.So)) return false; + if (!writeStdVecD(fs, s.PVT.Kro)) return false; + if (!writeStdVecD(fs, s.PVT.Sg)) return false; + if (!writeStdVecD(fs, s.PVT.Krg)) return false; + if (!writeStdVecD(fs, s.PVT.Sw)) return false; + if (!writeStdVecD(fs, s.PVT.Krw)) return false; + + // Base + if (!writeDouble(fs, s.Base.Pi)) return false; + if (!writeDouble(fs, s.Base.Cti)) return false; + if (!writeDouble(fs, s.Base.Cf)) return false; + if (!writeDouble(fs, s.Base.Soi)) return false; + if (!writeDouble(fs, s.Base.Sgi)) return false; + if (!writeDouble(fs, s.Base.Swi)) return false; + + if (!writeDouble(fs, s.Base.d)) return false; + if (!writeDouble(fs, s.Base.dt_Min)) return false; + if (!writeDouble(fs, s.Base.dt_Max)) return false; + + if (!writeDouble(fs, s.Base.k_ref)) return false; + if (!writeDouble(fs, s.Base.phi_ref)) return false; + if (!writeDouble(fs, s.Base.h_ref)) return false; + + return true; +} + +bool DatasetIO::readScene(std::ifstream& fs, PebiScene& s) +{ + if (!readInt(fs, s.version)) return false; + + // 网格几何 + if (!readInt(fs, s.D)) return false; + if (!readDouble(fs, s.GridControl)) return false; + if (!readStdVec2D(fs, s.Boundary)) return false; + if (!readStdVec2D(fs, s.VerticalWell)) return false; + if (!readStdVec2D(fs, s.HorizontalWell)) return false; + if (!readStdVec2D(fs, s.FractureVerticalWell)) return false; + if (!readStdVec3D(fs, s.MultistageFracturedHorizontalWell)) return false; + if (!readStdVec2D(fs, s.InclinedWell)) return false; + if (!readStdVec2D(fs, s.Fault)) return false; + + // 井顺序 + if (!readStdVecI(fs, s.wellType)) return false; + + // 井名 + size_t nName = 0; + if (!readSafeSize(fs, nName)) return false; + s.wellName.resize(nName); + for (size_t i = 0; i < nName; ++i) { + if (!readString(fs, s.wellName[i])) return false; + } + + // 求解器 + if (!readInt(fs, s.solverType)) return false; + + // Rate + if (!readStdVec2D(fs, s.Rate.t)) return false; + if (!readStdVec2D(fs, s.Rate.qo)) return false; + if (!readStdVec2D(fs, s.Rate.qg)) return false; + if (!readStdVec2D(fs, s.Rate.qw)) return false; + + // CS + if (!readStdVecD(fs, s.CS.C)) return false; + if (!readStdVecD(fs, s.CS.S)) return false; + + // 井流量段索引 + if (!readStdVecI(fs, s.wellFlowSectionIndex)) return false; + + // PVT + if (!readStdVecD(fs, s.PVT.p)) return false; + if (!readDouble(fs, s.PVT.pb)) return false; + + if (!readStdVecD(fs, s.PVT.Rso)) return false; + if (!readStdVecD(fs, s.PVT.Bo)) return false; + if (!readStdVecD(fs, s.PVT.Co)) return false; + if (!readStdVecD(fs, s.PVT.miuo)) return false; + if (!readStdVecD(fs, s.PVT.rouo)) return false; + + if (!readStdVecD(fs, s.PVT.Rv)) return false; + if (!readStdVecD(fs, s.PVT.Bg)) return false; + if (!readStdVecD(fs, s.PVT.Cg)) return false; + if (!readStdVecD(fs, s.PVT.miug)) return false; + if (!readStdVecD(fs, s.PVT.roug)) return false; + if (!readStdVecD(fs, s.PVT.Z)) return false; + + if (!readStdVecD(fs, s.PVT.Rsw)) return false; + if (!readStdVecD(fs, s.PVT.Bw)) return false; + if (!readStdVecD(fs, s.PVT.Cw)) return false; + if (!readStdVecD(fs, s.PVT.miuw)) return false; + if (!readStdVecD(fs, s.PVT.rouw)) return false; + + if (!readStdVecD(fs, s.PVT.V)) return false; + if (!readStdVecD(fs, s.PVT.k_kinitial)) return false; + if (!readStdVecD(fs, s.PVT.Cf_Cfinitial)) return false; + + if (!readStdVecD(fs, s.PVT.So)) return false; + if (!readStdVecD(fs, s.PVT.Kro)) return false; + if (!readStdVecD(fs, s.PVT.Sg)) return false; + if (!readStdVecD(fs, s.PVT.Krg)) return false; + if (!readStdVecD(fs, s.PVT.Sw)) return false; + if (!readStdVecD(fs, s.PVT.Krw)) return false; + + // Base + if (!readDouble(fs, s.Base.Pi)) return false; + if (!readDouble(fs, s.Base.Cti)) return false; + if (!readDouble(fs, s.Base.Cf)) return false; + if (!readDouble(fs, s.Base.Soi)) return false; + if (!readDouble(fs, s.Base.Sgi)) return false; + if (!readDouble(fs, s.Base.Swi)) return false; + + if (!readDouble(fs, s.Base.d)) return false; + if (!readDouble(fs, s.Base.dt_Min)) return false; + if (!readDouble(fs, s.Base.dt_Max)) return false; + + if (!readDouble(fs, s.Base.k_ref)) return false; + if (!readDouble(fs, s.Base.phi_ref)) return false; + if (!readDouble(fs, s.Base.h_ref)) return false; + + return true; +} + +// ---------------- gridOutput2 IO ---------------- + +bool DatasetIO::writeVec1D(std::ofstream& fs, const dVec1& v) +{ + if (!writeSafeSize(fs, v.size())) return false; + for (size_t i = 0; i < v.size(); ++i) if (!writeDouble(fs, v[i])) return false; + return true; +} + +bool DatasetIO::writeIVec1D(std::ofstream& fs, const iVec1& v) +{ + if (!writeSafeSize(fs, v.size())) return false; + for (size_t i = 0; i < v.size(); ++i) if (!writeInt(fs, v[i])) return false; + return true; +} + +bool DatasetIO::writeVec2D(std::ofstream& fs, const dVec2& v) +{ + if (!writeSafeSize(fs, v.size())) return false; + for (size_t i = 0; i < v.size(); ++i) if (!writeVec1D(fs, v[i])) return false; + return true; +} + +bool DatasetIO::writeIVec2D(std::ofstream& fs, const iVec2& v) +{ + if (!writeSafeSize(fs, v.size())) return false; + for (size_t i = 0; i < v.size(); ++i) if (!writeIVec1D(fs, v[i])) return false; + return true; +} + +bool DatasetIO::writeVec3D(std::ofstream& fs, const dVec3& v) +{ + if (!writeSafeSize(fs, v.size())) return false; + for (size_t i = 0; i < v.size(); ++i) if (!writeVec2D(fs, v[i])) return false; + return true; +} + +bool DatasetIO::readVec1D(std::ifstream& fs, dVec1& v) +{ + size_t sz = 0; + if (!readSafeSize(fs, sz)) return false; + v.resize(sz); + for (size_t i = 0; i < sz; ++i) if (!readDouble(fs, v[i])) return false; + return true; +} + +bool DatasetIO::readIVec1D(std::ifstream& fs, iVec1& v) +{ + size_t sz = 0; + if (!readSafeSize(fs, sz)) return false; + v.resize(sz); + for (size_t i = 0; i < sz; ++i) if (!readInt(fs, v[i])) return false; + return true; +} + +bool DatasetIO::readVec2D(std::ifstream& fs, dVec2& v) +{ + size_t rows = 0; + if (!readSafeSize(fs, rows)) return false; + v.resize(rows); + for (size_t i = 0; i < rows; ++i) if (!readVec1D(fs, v[i])) return false; + return true; +} + +bool DatasetIO::readIVec2D(std::ifstream& fs, iVec2& v) +{ + size_t rows = 0; + if (!readSafeSize(fs, rows)) return false; + v.resize(rows); + for (size_t i = 0; i < rows; ++i) if (!readIVec1D(fs, v[i])) return false; + return true; +} + +bool DatasetIO::readVec3D(std::ifstream& fs, dVec3& v) +{ + size_t depth = 0; + if (!readSafeSize(fs, depth)) return false; + v.resize(depth); + for (size_t i = 0; i < depth; ++i) if (!readVec2D(fs, v[i])) return false; + return true; +} + +bool DatasetIO::writeGridOutput2(std::ofstream& fs, const HX_NWTM_GRID_OUTPUT2& g) +{ + // top + if (!writeVec2D(fs, g.Trinodexy)) return false; + if (!writeVec1D(fs, g.Area)) return false; + if (!writeVec2D(fs, g.D)) return false; + + // ZhiJingNeiBianJie + if (!writeInt(fs, g.ZhiJingNeiBianJie.n)) return false; + if (!writeVec2D(fs, g.ZhiJingNeiBianJie.XiLinw)) return false; + if (!writeVec2D(fs, g.ZhiJingNeiBianJie.lw)) return false; + if (!writeVec2D(fs, g.ZhiJingNeiBianJie.dw)) return false; + if (!writeVec1D(fs, g.ZhiJingNeiBianJie.rw)) return false; + if (!writeIVec2D(fs, g.ZhiJingNeiBianJie.inwell)) return false; + + // LieFengJingNeiBianJie + if (!writeInt(fs, g.LieFengJingNeiBianJie.n)) return false; + if (!writeVec2D(fs, g.LieFengJingNeiBianJie.XiLinf)) return false; + if (!writeVec2D(fs, g.LieFengJingNeiBianJie.lf)) return false; + if (!writeVec2D(fs, g.LieFengJingNeiBianJie.df)) return false; + if (!writeVec1D(fs, g.LieFengJingNeiBianJie.xf)) return false; + if (!writeIVec2D(fs, g.LieFengJingNeiBianJie.infra)) return false; + + // DuoJiYaLieShuiPingJingNeiBianJie + if (!writeInt(fs, g.DuoJiYaLieShuiPingJingNeiBianJie.n)) return false; + if (!writeVec2D(fs, g.DuoJiYaLieShuiPingJingNeiBianJie.XiLinh)) return false; + if (!writeVec2D(fs, g.DuoJiYaLieShuiPingJingNeiBianJie.lh)) return false; + if (!writeVec2D(fs, g.DuoJiYaLieShuiPingJingNeiBianJie.dh)) return false; + if (!writeVec2D(fs, g.DuoJiYaLieShuiPingJingNeiBianJie.dsxf)) return false; + if (!writeIVec2D(fs, g.DuoJiYaLieShuiPingJingNeiBianJie.inhor)) return false; + if (!writeIVec1D(fs, g.DuoJiYaLieShuiPingJingNeiBianJie.nhor)) return false; + + // WaiBianJie + if (!writeInt(fs, g.WaiBianJie.n)) return false; + if (!writeVec2D(fs, g.WaiBianJie.WaiBianh)) return false; + if (!writeVec2D(fs, g.WaiBianJie.WaiBianl)) return false; + if (!writeVec2D(fs, g.WaiBianJie.WaiBiand)) return false; + + // NeiBuDuanCeng + if (!writeInt(fs, g.NeiBuDuanCeng.n)) return false; + if (!writeVec2D(fs, g.NeiBuDuanCeng.faultb1)) return false; + if (!writeVec2D(fs, g.NeiBuDuanCeng.faultb2)) return false; + if (!writeVec2D(fs, g.NeiBuDuanCeng.faultl1)) return false; + if (!writeVec2D(fs, g.NeiBuDuanCeng.faultd1)) return false; + + // YuChuLiJuZhen + if (!writeIVec1D(fs, g.YuChuLiJuZhen.ia)) return false; + if (!writeIVec1D(fs, g.YuChuLiJuZhen.ja)) return false; + if (!writeIVec2D(fs, g.YuChuLiJuZhen.nzeros)) return false; + if (!writeInt(fs, g.YuChuLiJuZhen.numk)) return false; + + return true; +} + +bool DatasetIO::readGridOutput2(std::ifstream& fs, HX_NWTM_GRID_OUTPUT2& g) +{ + if (!readVec2D(fs, g.Trinodexy)) return false; + if (!readVec1D(fs, g.Area)) return false; + if (!readVec2D(fs, g.D)) return false; + + if (!readInt(fs, g.ZhiJingNeiBianJie.n)) return false; + if (!readVec2D(fs, g.ZhiJingNeiBianJie.XiLinw)) return false; + if (!readVec2D(fs, g.ZhiJingNeiBianJie.lw)) return false; + if (!readVec2D(fs, g.ZhiJingNeiBianJie.dw)) return false; + if (!readVec1D(fs, g.ZhiJingNeiBianJie.rw)) return false; + if (!readIVec2D(fs, g.ZhiJingNeiBianJie.inwell)) return false; + + if (!readInt(fs, g.LieFengJingNeiBianJie.n)) return false; + if (!readVec2D(fs, g.LieFengJingNeiBianJie.XiLinf)) return false; + if (!readVec2D(fs, g.LieFengJingNeiBianJie.lf)) return false; + if (!readVec2D(fs, g.LieFengJingNeiBianJie.df)) return false; + if (!readVec1D(fs, g.LieFengJingNeiBianJie.xf)) return false; + if (!readIVec2D(fs, g.LieFengJingNeiBianJie.infra)) return false; + + if (!readInt(fs, g.DuoJiYaLieShuiPingJingNeiBianJie.n)) return false; + if (!readVec2D(fs, g.DuoJiYaLieShuiPingJingNeiBianJie.XiLinh)) return false; + if (!readVec2D(fs, g.DuoJiYaLieShuiPingJingNeiBianJie.lh)) return false; + if (!readVec2D(fs, g.DuoJiYaLieShuiPingJingNeiBianJie.dh)) return false; + if (!readVec2D(fs, g.DuoJiYaLieShuiPingJingNeiBianJie.dsxf)) return false; + if (!readIVec2D(fs, g.DuoJiYaLieShuiPingJingNeiBianJie.inhor)) return false; + if (!readIVec1D(fs, g.DuoJiYaLieShuiPingJingNeiBianJie.nhor)) return false; + + if (!readInt(fs, g.WaiBianJie.n)) return false; + if (!readVec2D(fs, g.WaiBianJie.WaiBianh)) return false; + if (!readVec2D(fs, g.WaiBianJie.WaiBianl)) return false; + if (!readVec2D(fs, g.WaiBianJie.WaiBiand)) return false; + + if (!readInt(fs, g.NeiBuDuanCeng.n)) return false; + if (!readVec2D(fs, g.NeiBuDuanCeng.faultb1)) return false; + if (!readVec2D(fs, g.NeiBuDuanCeng.faultb2)) return false; + if (!readVec2D(fs, g.NeiBuDuanCeng.faultl1)) return false; + if (!readVec2D(fs, g.NeiBuDuanCeng.faultd1)) return false; + + if (!readIVec1D(fs, g.YuChuLiJuZhen.ia)) return false; + if (!readIVec1D(fs, g.YuChuLiJuZhen.ja)) return false; + if (!readIVec2D(fs, g.YuChuLiJuZhen.nzeros)) return false; + if (!readInt(fs, g.YuChuLiJuZhen.numk)) return false; + + return true; +} + +// ---------------- FNV-1a 64 hash ---------------- + +unsigned long long DatasetIO::fnv1a64_init() +{ + return 14695981039346656037ULL; +} + +void DatasetIO::fnv1a64_update(unsigned long long& h, const void* data, size_t len) +{ + const unsigned char* p = (const unsigned char*)data; + const unsigned long long prime = 1099511628211ULL; + for (size_t i = 0; i < len; ++i) { + h ^= (unsigned long long)p[i]; + h *= prime; + } +} + +void DatasetIO::hashInt(unsigned long long& h, int v) { fnv1a64_update(h, &v, sizeof(int)); } +void DatasetIO::hashUInt(unsigned long long& h, unsigned int v) { fnv1a64_update(h, &v, sizeof(unsigned int)); } +void DatasetIO::hashDouble(unsigned long long& h, double v) { fnv1a64_update(h, &v, sizeof(double)); } + +void DatasetIO::hashVec1D(unsigned long long& h, const dVec1& v) +{ + unsigned int sz = (unsigned int)v.size(); + hashUInt(h, sz); + for (size_t i = 0; i < v.size(); ++i) hashDouble(h, v[i]); +} + +void DatasetIO::hashIVec1D(unsigned long long& h, const iVec1& v) +{ + unsigned int sz = (unsigned int)v.size(); + hashUInt(h, sz); + for (size_t i = 0; i < v.size(); ++i) hashInt(h, v[i]); +} + +void DatasetIO::hashVec2D(unsigned long long& h, const dVec2& v) +{ + unsigned int rows = (unsigned int)v.size(); + hashUInt(h, rows); + for (size_t i = 0; i < v.size(); ++i) { + unsigned int cols = (unsigned int)v[i].size(); + hashUInt(h, cols); + for (size_t j = 0; j < v[i].size(); ++j) hashDouble(h, v[i][j]); + } +} + +void DatasetIO::hashIVec2D(unsigned long long& h, const iVec2& v) +{ + unsigned int rows = (unsigned int)v.size(); + hashUInt(h, rows); + for (size_t i = 0; i < v.size(); ++i) { + unsigned int cols = (unsigned int)v[i].size(); + hashUInt(h, cols); + for (size_t j = 0; j < v[i].size(); ++j) hashInt(h, v[i][j]); + } +} + +void DatasetIO::hashVec3D(unsigned long long& h, const dVec3& v) +{ + unsigned int depth = (unsigned int)v.size(); + hashUInt(h, depth); + for (size_t i = 0; i < v.size(); ++i) { + hashVec2D(h, v[i]); + } +} + +unsigned long long DatasetIO::computeSceneKey64(const HX_NWTM_GRID_INPUT& in) +{ + unsigned long long h = fnv1a64_init(); + hashInt(h, in.D); + hashDouble(h, in.GridControl); + + hashVec2D(h, in.Boundary); + hashVec2D(h, in.VerticalWell); + hashVec2D(h, in.HorizontalWell); + hashVec2D(h, in.FractureVerticalWell); + hashVec3D(h, in.MultistageFracturedHorizontalWell); + hashVec2D(h, in.InclinedWell); + hashVec2D(h, in.Fault); + return h; +} \ No newline at end of file diff --git a/ML/Training/Training/DatasetIO.h b/ML/Training/Training/DatasetIO.h new file mode 100644 index 0000000..55f6eca --- /dev/null +++ b/ML/Training/Training/DatasetIO.h @@ -0,0 +1,127 @@ +#pragma once +#ifndef DATASET_IO_H +#define DATASET_IO_H + +#include +#include +#include +#include + +#include "pch.h" +#include "SceneIO.h" // 为了用 PebiScene + +// DatasetIO：把 (scene + gridOutput2) 打包成一个 dataset.bin，供后续直接加载 +class DatasetIO +{ +public: + // 保存 dataset（包含 scene + gridOutput2 + header(sceneKey)） + static bool saveDataset(const std::string& filename, + const PebiScene& scene, + const HX_NWTM_GRID_INPUT& gridInput, + const HX_NWTM_GRID_OUTPUT2& grid); + + // 加载 dataset（读出 scene + gridOutput2，且校验 header） + // strictVerify=true：会用读出来的 scene 重新计算 sceneKey 与 header 对比（防止文件损坏/写错） + static bool loadDataset(const std::string& filename, + PebiScene& outScene, + HX_NWTM_GRID_OUTPUT2& outGrid, + bool strictVerify = true); + + // 只读 header 并返回 sceneKey（用于与外部 scene.bin 比对串场景） + static bool readDatasetSceneKey(const std::string& filename, + unsigned long long& outSceneKey, + unsigned int& outCells, + unsigned int& outWells, + int& outSolverType); + + static bool fileExists(const std::string& filename); + + // 计算 sceneKey：与 GridCacheIO 相同算法（FNV-1a 64） + static unsigned long long computeSceneKey64(const HX_NWTM_GRID_INPUT& in); + +private: + // ---------------- 文件头 ---------------- + struct DatasetHeader + { + unsigned int magic; // 'DTHX' 等 + unsigned int version; // 格式版本 + unsigned int reserved; // 预留 + unsigned long long sceneKey; // 串场景指纹 + unsigned int nCells; // cell 数 + unsigned int nWells; // 井数（从 Rate.t.size() 或 scene.wellName.size() 推） + int solverType; // 求解器类型 + }; + + enum { MAGIC = 0x58485444 }; // 'DTHX' little-endian + enum { VERSION = 1 }; + + // 从 scene 生成 gridInput（用于 strictVerify 重新算 sceneKey） + static HX_NWTM_GRID_INPUT buildGridInputFromScene(const PebiScene& scene); + + // ---------------- 基础 IO（安全尺寸） ---------------- + static bool writeUInt(std::ofstream& fs, unsigned int v); + static bool writeULL(std::ofstream& fs, unsigned long long v); + static bool writeInt(std::ofstream& fs, int v); + static bool writeDouble(std::ofstream& fs, double v); + + static bool readUInt(std::ifstream& fs, unsigned int& v); + static bool readULL(std::ifstream& fs, unsigned long long& v); + static bool readInt(std::ifstream& fs, int& v); + static bool readDouble(std::ifstream& fs, double& v); + + static bool writeSafeSize(std::ofstream& fs, size_t sz); + static bool readSafeSize(std::ifstream& fs, size_t& sz); + + // ---------------- string/vector/scene IO ---------------- + static bool writeString(std::ofstream& fs, const std::string& s); + static bool readString(std::ifstream& fs, std::string& s); + + static bool writeStdVecD(std::ofstream& fs, const std::vector& v); + static bool readStdVecD(std::ifstream& fs, std::vector& v); + + static bool writeStdVecI(std::ofstream& fs, const std::vector& v); + static bool readStdVecI(std::ifstream& fs, std::vector& v); + + static bool writeStdVec2D(std::ofstream& fs, const std::vector >& v); + static bool readStdVec2D(std::ifstream& fs, std::vector >& v); + + static bool writeStdVec3D(std::ofstream& fs, const std::vector > >& v); + static bool readStdVec3D(std::ifstream& fs, std::vector > >& v); + + static bool writeScene(std::ofstream& fs, const PebiScene& scene); + static bool readScene(std::ifstream& fs, PebiScene& scene); + + // ---------------- gridOutput2 IO ---------------- + static bool writeVec1D(std::ofstream& fs, const dVec1& v); + static bool writeIVec1D(std::ofstream& fs, const iVec1& v); + static bool writeVec2D(std::ofstream& fs, const dVec2& v); + static bool writeIVec2D(std::ofstream& fs, const iVec2& v); + static bool writeVec3D(std::ofstream& fs, const dVec3& v); + + static bool readVec1D(std::ifstream& fs, dVec1& v); + static bool readIVec1D(std::ifstream& fs, iVec1& v); + static bool readVec2D(std::ifstream& fs, dVec2& v); + static bool readIVec2D(std::ifstream& fs, iVec2& v); + static bool readVec3D(std::ifstream& fs, dVec3& v); + + static bool writeGridOutput2(std::ofstream& fs, const HX_NWTM_GRID_OUTPUT2& g); + static bool readGridOutput2(std::ifstream& fs, HX_NWTM_GRID_OUTPUT2& g); + + static bool writeHeader(std::ofstream& fs, const DatasetHeader& hdr); + static bool readHeader(std::ifstream& fs, DatasetHeader& hdr); + + // ---------------- FNV-1a 64 hash ---------------- + static unsigned long long fnv1a64_init(); + static void fnv1a64_update(unsigned long long& h, const void* data, size_t len); + + static void hashInt(unsigned long long& h, int v); + static void hashUInt(unsigned long long& h, unsigned int v); + static void hashDouble(unsigned long long& h, double v); + static void hashVec1D(unsigned long long& h, const dVec1& v); + static void hashIVec1D(unsigned long long& h, const iVec1& v); + static void hashVec2D(unsigned long long& h, const dVec2& v); + static void hashIVec2D(unsigned long long& h, const iVec2& v); + static void hashVec3D(unsigned long long& h, const dVec3& v); +}; + +#endif // DATASET_IO_H diff --git a/ML/Training/Training/GridCacheIO.cpp b/ML/Training/Training/GridCacheIO.cpp new file mode 100644 index 0000000..84f473e --- /dev/null +++ b/ML/Training/Training/GridCacheIO.cpp @@ -0,0 +1,516 @@ +#include "GridCacheIO.h" + +// ---------------- public ---------------- + +bool GridCacheIO::fileExists(const std::string& filename) +{ + std::ifstream fs(filename.c_str(), std::ios::binary); + return fs.good(); +} + +bool GridCacheIO::saveGrid(const std::string& filename, + HX_NWTM_GRID_OUTPUT2 grid, + const HX_NWTM_GRID_INPUT& gridInput) +{ + // 写入前先清洗，避免把 DLL 未初始化的 n 写进缓存 + sanitizeGridOutput2(grid); + + std::ofstream fs(filename.c_str(), std::ios::binary); + if (!fs) { + std::cerr << "GridCacheIO::saveGrid: 无法打开 " << filename << std::endl; + return false; + } + + GridCacheHeader hdr; + hdr.magic = (unsigned int)MAGIC; + hdr.version = (unsigned int)VERSION; + hdr.reserved = 0; + hdr.sceneKey = computeSceneKey64(gridInput); + hdr.nCells = (unsigned int)grid.Trinodexy.size(); + + if (!writeHeader(fs, hdr)) { + std::cerr << "GridCacheIO::saveGrid: 写 header 失败\n"; + return false; + } + + if (!writeGridOutput2(fs, grid)) { + std::cerr << "GridCacheIO::saveGrid: 写 gridOutput2 失败\n"; + return false; + } + + return fs.good(); +} + +bool GridCacheIO::isCacheMatch(const std::string& filename, + const HX_NWTM_GRID_INPUT& expectedGridInput) +{ + std::ifstream fs(filename.c_str(), std::ios::binary); + if (!fs) return false; + + GridCacheHeader hdr; + if (!readHeader(fs, hdr)) return false; + if (hdr.magic != (unsigned int)MAGIC) return false; + if (hdr.version != (unsigned int)VERSION) return false; + + unsigned long long key = computeSceneKey64(expectedGridInput); + return (hdr.sceneKey == key); +} + +bool GridCacheIO::loadGrid(const std::string& filename, + HX_NWTM_GRID_OUTPUT2& grid, + const HX_NWTM_GRID_INPUT& expectedGridInput) +{ + std::ifstream fs(filename.c_str(), std::ios::binary); + if (!fs) { + std::cerr << "GridCacheIO::loadGrid: 无法打开 " << filename << std::endl; + return false; + } + + GridCacheHeader hdr; + if (!readHeader(fs, hdr)) { + std::cerr << "GridCacheIO::loadGrid: 读 header 失败\n"; + return false; + } + + if (hdr.magic != (unsigned int)MAGIC) { + std::cerr << "GridCacheIO::loadGrid: magic 错误\n"; + return false; + } + if (hdr.version != (unsigned int)VERSION) { + std::cerr << "GridCacheIO::loadGrid: version 错误\n"; + return false; + } + + unsigned long long expectedKey = computeSceneKey64(expectedGridInput); + if (hdr.sceneKey != expectedKey) { + std::cerr << "GridCacheIO::loadGrid: sceneKey 不匹配\n"; + return false; + } + + if (!readGridOutput2(fs, grid)) { + std::cerr << "GridCacheIO::loadGrid: 读 gridOutput2 失败\n"; + return false; + } + + // 读取后清洗，保证 n 字段稳定 + sanitizeGridOutput2(grid); + + // 快速一致性检查 + if (hdr.nCells != (unsigned int)grid.Trinodexy.size()) { + std::cerr << "GridCacheIO::loadGrid: nCells 不匹配，缓存可能损坏\n"; + return false; + } + + return fs.good(); +} + +// ---------------- header io ---------------- + +bool GridCacheIO::writeHeader(std::ofstream& fs, const GridCacheHeader& hdr) +{ + if (!writeUInt(fs, hdr.magic)) return false; + if (!writeUInt(fs, hdr.version)) return false; + if (!writeUInt(fs, hdr.reserved)) return false; + if (!writeULL(fs, hdr.sceneKey)) return false; + if (!writeUInt(fs, hdr.nCells)) return false; + return true; +} + +bool GridCacheIO::readHeader(std::ifstream& fs, GridCacheHeader& hdr) +{ + if (!readUInt(fs, hdr.magic)) return false; + if (!readUInt(fs, hdr.version)) return false; + if (!readUInt(fs, hdr.reserved)) return false; + if (!readULL(fs, hdr.sceneKey)) return false; + if (!readUInt(fs, hdr.nCells)) return false; + return true; +} + +// ---------------- sanitize ---------------- + +void GridCacheIO::sanitizeGridOutput2(HX_NWTM_GRID_OUTPUT2& g) +{ + // 这些 n 字段在日志里出现了“随机大数”，非常像 DLL 未初始化/非确定性字段 + // 最安全策略：统一用 vector size 来覆盖它们，避免缓存写入/读出后污染流程。 + + g.ZhiJingNeiBianJie.n = (int)g.ZhiJingNeiBianJie.XiLinw.size(); + g.LieFengJingNeiBianJie.n = (int)g.LieFengJingNeiBianJie.XiLinf.size(); + g.DuoJiYaLieShuiPingJingNeiBianJie.n = (int)g.DuoJiYaLieShuiPingJingNeiBianJie.XiLinh.size(); + g.WaiBianJie.n = (int)g.WaiBianJie.WaiBianh.size(); + + // 内部断层：用 faultb1 的行数作为 n（你对比时 vector 内容一致，而 n 不一致） + g.NeiBuDuanCeng.n = (int)g.NeiBuDuanCeng.faultb1.size(); + + // 预处理矩阵：numk 如果为 0 或异常，可用 ia 大小兜底（这里不强制覆盖，只做兜底） + if (g.YuChuLiJuZhen.numk < 0) g.YuChuLiJuZhen.numk = 0; +} + +// ---------------- FNV1a hash ---------------- + +unsigned long long GridCacheIO::fnv1a64_init() +{ + return 14695981039346656037ULL; // offset basis +} + +void GridCacheIO::fnv1a64_update(unsigned long long& h, const void* data, size_t len) +{ + const unsigned char* p = (const unsigned char*)data; + const unsigned long long prime = 1099511628211ULL; + for (size_t i = 0; i < len; ++i) { + h ^= (unsigned long long)p[i]; + h *= prime; + } +} + +void GridCacheIO::hashInt(unsigned long long& h, int v) +{ + fnv1a64_update(h, &v, sizeof(int)); +} + +void GridCacheIO::hashUInt(unsigned long long& h, unsigned int v) +{ + fnv1a64_update(h, &v, sizeof(unsigned int)); +} + +void GridCacheIO::hashDouble(unsigned long long& h, double v) +{ + fnv1a64_update(h, &v, sizeof(double)); +} + +void GridCacheIO::hashVec1D(unsigned long long& h, const dVec1& v) +{ + unsigned int sz = (unsigned int)v.size(); + hashUInt(h, sz); + for (size_t i = 0; i < v.size(); ++i) hashDouble(h, v[i]); +} + +void GridCacheIO::hashIVec1D(unsigned long long& h, const iVec1& v) +{ + unsigned int sz = (unsigned int)v.size(); + hashUInt(h, sz); + for (size_t i = 0; i < v.size(); ++i) hashInt(h, v[i]); +} + +void GridCacheIO::hashVec2D(unsigned long long& h, const dVec2& v) +{ + unsigned int rows = (unsigned int)v.size(); + hashUInt(h, rows); + for (size_t i = 0; i < v.size(); ++i) { + unsigned int cols = (unsigned int)v[i].size(); + hashUInt(h, cols); + for (size_t j = 0; j < v[i].size(); ++j) hashDouble(h, v[i][j]); + } +} + +void GridCacheIO::hashIVec2D(unsigned long long& h, const iVec2& v) +{ + unsigned int rows = (unsigned int)v.size(); + hashUInt(h, rows); + for (size_t i = 0; i < v.size(); ++i) { + unsigned int cols = (unsigned int)v[i].size(); + hashUInt(h, cols); + for (size_t j = 0; j < v[i].size(); ++j) hashInt(h, v[i][j]); + } +} + +void GridCacheIO::hashVec3D(unsigned long long& h, const dVec3& v) +{ + unsigned int depth = (unsigned int)v.size(); + hashUInt(h, depth); + for (size_t i = 0; i < v.size(); ++i) { + hashVec2D(h, v[i]); + } +} + +unsigned long long GridCacheIO::computeSceneKey64(const HX_NWTM_GRID_INPUT& in) +{ + unsigned long long h = fnv1a64_init(); + + // 只 hash 会影响网格生成的关键输入（串场景保护） + hashInt(h, in.D); + hashDouble(h, in.GridControl); + + hashVec2D(h, in.Boundary); + hashVec2D(h, in.VerticalWell); + hashVec2D(h, in.HorizontalWell); + hashVec2D(h, in.FractureVerticalWell); + hashVec3D(h, in.MultistageFracturedHorizontalWell); + hashVec2D(h, in.InclinedWell); + hashVec2D(h, in.Fault); + + return h; +} + +// ---------------- primitive io ---------------- + +bool GridCacheIO::writeUInt(std::ofstream& fs, unsigned int v) +{ + fs.write((const char*)&v, sizeof(unsigned int)); + return fs.good(); +} + +bool GridCacheIO::writeULL(std::ofstream& fs, unsigned long long v) +{ + fs.write((const char*)&v, sizeof(unsigned long long)); + return fs.good(); +} + +bool GridCacheIO::writeInt(std::ofstream& fs, int v) +{ + fs.write((const char*)&v, sizeof(int)); + return fs.good(); +} + +bool GridCacheIO::writeDouble(std::ofstream& fs, double v) +{ + fs.write((const char*)&v, sizeof(double)); + return fs.good(); +} + +bool GridCacheIO::readUInt(std::ifstream& fs, unsigned int& v) +{ + fs.read((char*)&v, sizeof(unsigned int)); + return fs.good(); +} + +bool GridCacheIO::readULL(std::ifstream& fs, unsigned long long& v) +{ + fs.read((char*)&v, sizeof(unsigned long long)); + return fs.good(); +} + +bool GridCacheIO::readInt(std::ifstream& fs, int& v) +{ + fs.read((char*)&v, sizeof(int)); + return fs.good(); +} + +bool GridCacheIO::readDouble(std::ifstream& fs, double& v) +{ + fs.read((char*)&v, sizeof(double)); + return fs.good(); +} + +bool GridCacheIO::writeSafeSize(std::ofstream& fs, size_t sz) +{ + unsigned int u = (unsigned int)sz; // 存为 uint32 + return writeUInt(fs, u); +} + +bool GridCacheIO::readSafeSize(std::ifstream& fs, size_t& sz) +{ + unsigned int u = 0; + if (!readUInt(fs, u)) return false; + + // 防御：防止读到异常大值导致内存爆炸 + const unsigned int MAX_SIZE = 100000000U; + if (u > MAX_SIZE) return false; + + sz = (size_t)u; + return true; +} + +// ---------------- vector io ---------------- + +bool GridCacheIO::writeVec1D(std::ofstream& fs, const dVec1& v) +{ + if (!writeSafeSize(fs, v.size())) return false; + for (size_t i = 0; i < v.size(); ++i) { + if (!writeDouble(fs, v[i])) return false; + } + return true; +} + +bool GridCacheIO::writeIVec1D(std::ofstream& fs, const iVec1& v) +{ + if (!writeSafeSize(fs, v.size())) return false; + for (size_t i = 0; i < v.size(); ++i) { + if (!writeInt(fs, v[i])) return false; + } + return true; +} + +bool GridCacheIO::writeVec2D(std::ofstream& fs, const dVec2& v) +{ + if (!writeSafeSize(fs, v.size())) return false; + for (size_t i = 0; i < v.size(); ++i) { + if (!writeVec1D(fs, v[i])) return false; + } + return true; +} + +bool GridCacheIO::writeIVec2D(std::ofstream& fs, const iVec2& v) +{ + if (!writeSafeSize(fs, v.size())) return false; + for (size_t i = 0; i < v.size(); ++i) { + if (!writeIVec1D(fs, v[i])) return false; + } + return true; +} + +bool GridCacheIO::writeVec3D(std::ofstream& fs, const dVec3& v) +{ + if (!writeSafeSize(fs, v.size())) return false; + for (size_t i = 0; i < v.size(); ++i) { + if (!writeVec2D(fs, v[i])) return false; + } + return true; +} + +bool GridCacheIO::readVec1D(std::ifstream& fs, dVec1& v) +{ + size_t sz = 0; + if (!readSafeSize(fs, sz)) return false; + v.resize(sz); + for (size_t i = 0; i < sz; ++i) { + if (!readDouble(fs, v[i])) return false; + } + return true; +} + +bool GridCacheIO::readIVec1D(std::ifstream& fs, iVec1& v) +{ + size_t sz = 0; + if (!readSafeSize(fs, sz)) return false; + v.resize(sz); + for (size_t i = 0; i < sz; ++i) { + if (!readInt(fs, v[i])) return false; + } + return true; +} + +bool GridCacheIO::readVec2D(std::ifstream& fs, dVec2& v) +{ + size_t rows = 0; + if (!readSafeSize(fs, rows)) return false; + v.resize(rows); + for (size_t i = 0; i < rows; ++i) { + if (!readVec1D(fs, v[i])) return false; + } + return true; +} + +bool GridCacheIO::readIVec2D(std::ifstream& fs, iVec2& v) +{ + size_t rows = 0; + if (!readSafeSize(fs, rows)) return false; + v.resize(rows); + for (size_t i = 0; i < rows; ++i) { + if (!readIVec1D(fs, v[i])) return false; + } + return true; +} + +bool GridCacheIO::readVec3D(std::ifstream& fs, dVec3& v) +{ + size_t depth = 0; + if (!readSafeSize(fs, depth)) return false; + v.resize(depth); + for (size_t i = 0; i < depth; ++i) { + if (!readVec2D(fs, v[i])) return false; + } + return true; +} + +// ---------------- gridOutput2 io ---------------- + +bool GridCacheIO::writeGridOutput2(std::ofstream& fs, const HX_NWTM_GRID_OUTPUT2& g) +{ + // 顶层 + if (!writeVec2D(fs, g.Trinodexy)) return false; + if (!writeVec1D(fs, g.Area)) return false; + if (!writeVec2D(fs, g.D)) return false; + + // ZhiJingNeiBianJie + if (!writeInt(fs, g.ZhiJingNeiBianJie.n)) return false; + if (!writeVec2D(fs, g.ZhiJingNeiBianJie.XiLinw)) return false; + if (!writeVec2D(fs, g.ZhiJingNeiBianJie.lw)) return false; + if (!writeVec2D(fs, g.ZhiJingNeiBianJie.dw)) return false; + if (!writeVec1D(fs, g.ZhiJingNeiBianJie.rw)) return false; + if (!writeIVec2D(fs, g.ZhiJingNeiBianJie.inwell)) return false; + + // LieFengJingNeiBianJie + if (!writeInt(fs, g.LieFengJingNeiBianJie.n)) return false; + if (!writeVec2D(fs, g.LieFengJingNeiBianJie.XiLinf)) return false; + if (!writeVec2D(fs, g.LieFengJingNeiBianJie.lf)) return false; + if (!writeVec2D(fs, g.LieFengJingNeiBianJie.df)) return false; + if (!writeVec1D(fs, g.LieFengJingNeiBianJie.xf)) return false; + if (!writeIVec2D(fs, g.LieFengJingNeiBianJie.infra)) return false; + + // DuoJiYaLieShuiPingJingNeiBianJie + if (!writeInt(fs, g.DuoJiYaLieShuiPingJingNeiBianJie.n)) return false; + if (!writeVec2D(fs, g.DuoJiYaLieShuiPingJingNeiBianJie.XiLinh)) return false; + if (!writeVec2D(fs, g.DuoJiYaLieShuiPingJingNeiBianJie.lh)) return false; + if (!writeVec2D(fs, g.DuoJiYaLieShuiPingJingNeiBianJie.dh)) return false; + if (!writeVec2D(fs, g.DuoJiYaLieShuiPingJingNeiBianJie.dsxf)) return false; + if (!writeIVec2D(fs, g.DuoJiYaLieShuiPingJingNeiBianJie.inhor)) return false; + if (!writeIVec1D(fs, g.DuoJiYaLieShuiPingJingNeiBianJie.nhor)) return false; + + // WaiBianJie + if (!writeInt(fs, g.WaiBianJie.n)) return false; + if (!writeVec2D(fs, g.WaiBianJie.WaiBianh)) return false; + if (!writeVec2D(fs, g.WaiBianJie.WaiBianl)) return false; + if (!writeVec2D(fs, g.WaiBianJie.WaiBiand)) return false; + + // NeiBuDuanCeng + if (!writeInt(fs, g.NeiBuDuanCeng.n)) return false; + if (!writeVec2D(fs, g.NeiBuDuanCeng.faultb1)) return false; + if (!writeVec2D(fs, g.NeiBuDuanCeng.faultb2)) return false; + if (!writeVec2D(fs, g.NeiBuDuanCeng.faultl1)) return false; + if (!writeVec2D(fs, g.NeiBuDuanCeng.faultd1)) return false; + + // YuChuLiJuZhen + if (!writeIVec1D(fs, g.YuChuLiJuZhen.ia)) return false; + if (!writeIVec1D(fs, g.YuChuLiJuZhen.ja)) return false; + if (!writeIVec2D(fs, g.YuChuLiJuZhen.nzeros)) return false; + if (!writeInt(fs, g.YuChuLiJuZhen.numk)) return false; + + return true; +} + +bool GridCacheIO::readGridOutput2(std::ifstream& fs, HX_NWTM_GRID_OUTPUT2& g) +{ + if (!readVec2D(fs, g.Trinodexy)) return false; + if (!readVec1D(fs, g.Area)) return false; + if (!readVec2D(fs, g.D)) return false; + + if (!readInt(fs, g.ZhiJingNeiBianJie.n)) return false; + if (!readVec2D(fs, g.ZhiJingNeiBianJie.XiLinw)) return false; + if (!readVec2D(fs, g.ZhiJingNeiBianJie.lw)) return false; + if (!readVec2D(fs, g.ZhiJingNeiBianJie.dw)) return false; + if (!readVec1D(fs, g.ZhiJingNeiBianJie.rw)) return false; + if (!readIVec2D(fs, g.ZhiJingNeiBianJie.inwell)) return false; + + if (!readInt(fs, g.LieFengJingNeiBianJie.n)) return false; + if (!readVec2D(fs, g.LieFengJingNeiBianJie.XiLinf)) return false; + if (!readVec2D(fs, g.LieFengJingNeiBianJie.lf)) return false; + if (!readVec2D(fs, g.LieFengJingNeiBianJie.df)) return false; + if (!readVec1D(fs, g.LieFengJingNeiBianJie.xf)) return false; + if (!readIVec2D(fs, g.LieFengJingNeiBianJie.infra)) return false; + + if (!readInt(fs, g.DuoJiYaLieShuiPingJingNeiBianJie.n)) return false; + if (!readVec2D(fs, g.DuoJiYaLieShuiPingJingNeiBianJie.XiLinh)) return false; + if (!readVec2D(fs, g.DuoJiYaLieShuiPingJingNeiBianJie.lh)) return false; + if (!readVec2D(fs, g.DuoJiYaLieShuiPingJingNeiBianJie.dh)) return false; + if (!readVec2D(fs, g.DuoJiYaLieShuiPingJingNeiBianJie.dsxf)) return false; + if (!readIVec2D(fs, g.DuoJiYaLieShuiPingJingNeiBianJie.inhor)) return false; + if (!readIVec1D(fs, g.DuoJiYaLieShuiPingJingNeiBianJie.nhor)) return false; + + if (!readInt(fs, g.WaiBianJie.n)) return false; + if (!readVec2D(fs, g.WaiBianJie.WaiBianh)) return false; + if (!readVec2D(fs, g.WaiBianJie.WaiBianl)) return false; + if (!readVec2D(fs, g.WaiBianJie.WaiBiand)) return false; + + if (!readInt(fs, g.NeiBuDuanCeng.n)) return false; + if (!readVec2D(fs, g.NeiBuDuanCeng.faultb1)) return false; + if (!readVec2D(fs, g.NeiBuDuanCeng.faultb2)) return false; + if (!readVec2D(fs, g.NeiBuDuanCeng.faultl1)) return false; + if (!readVec2D(fs, g.NeiBuDuanCeng.faultd1)) return false; + + if (!readIVec1D(fs, g.YuChuLiJuZhen.ia)) return false; + if (!readIVec1D(fs, g.YuChuLiJuZhen.ja)) return false; + if (!readIVec2D(fs, g.YuChuLiJuZhen.nzeros)) return false; + if (!readInt(fs, g.YuChuLiJuZhen.numk)) return false; + + return true; +} \ No newline at end of file diff --git a/ML/Training/Training/GridCacheIO.h b/ML/Training/Training/GridCacheIO.h new file mode 100644 index 0000000..e67f530 --- /dev/null +++ b/ML/Training/Training/GridCacheIO.h @@ -0,0 +1,101 @@ +#pragma once +#ifndef GRID_CACHE_IO_H +#define GRID_CACHE_IO_H + +#include +#include +#include + +#include "pch.h" + +class GridCacheIO +{ +public: + // 保存网格（写入 header + gridOutput2） + // 注意：会自动“清洗”一些不稳定的 n 字段，避免写入垃圾值 + static bool saveGrid(const std::string& filename, + HX_NWTM_GRID_OUTPUT2 grid, // 这里用值传递：方便内部清洗 + const HX_NWTM_GRID_INPUT& gridInput); + + // 读取网格（读取后会检查 sceneKey 是否匹配） + // 注意：读取后也会做一次“清洗”，保证 n 字段稳定 + static bool loadGrid(const std::string& filename, + HX_NWTM_GRID_OUTPUT2& grid, + const HX_NWTM_GRID_INPUT& expectedGridInput); + + // 只检查 header 的 sceneKey 是否匹配（不读完整网格） + static bool isCacheMatch(const std::string& filename, + const HX_NWTM_GRID_INPUT& expectedGridInput); + + static bool fileExists(const std::string& filename); + +private: + // ---------------- 文件头 ---------------- + struct GridCacheHeader + { + unsigned int magic; // 'GCHX' + unsigned int version; // 格式版本 + unsigned int reserved; // 保留 + unsigned long long sceneKey; // 串场景指纹 + unsigned int nCells; // 快速检查 + }; + + enum { MAGIC = 0x58484347 }; // 'GCHX' little-endian + enum { VERSION = 1 }; + +private: + // 计算串场景 key + static unsigned long long computeSceneKey64(const HX_NWTM_GRID_INPUT& in); + + // FNV-1a 64 + static unsigned long long fnv1a64_init(); + static void fnv1a64_update(unsigned long long& h, const void* data, size_t len); + + // hash 辅助 + static void hashInt(unsigned long long& h, int v); + static void hashUInt(unsigned long long& h, unsigned int v); + static void hashDouble(unsigned long long& h, double v); + static void hashVec1D(unsigned long long& h, const dVec1& v); + static void hashIVec1D(unsigned long long& h, const iVec1& v); + static void hashVec2D(unsigned long long& h, const dVec2& v); + static void hashIVec2D(unsigned long long& h, const iVec2& v); + static void hashVec3D(unsigned long long& h, const dVec3& v); + + // IO 基元 + static bool writeUInt(std::ofstream& fs, unsigned int v); + static bool writeULL(std::ofstream& fs, unsigned long long v); + static bool writeInt(std::ofstream& fs, int v); + static bool writeDouble(std::ofstream& fs, double v); + + static bool readUInt(std::ifstream& fs, unsigned int& v); + static bool readULL(std::ifstream& fs, unsigned long long& v); + static bool readInt(std::ifstream& fs, int& v); + static bool readDouble(std::ifstream& fs, double& v); + + static bool writeSafeSize(std::ofstream& fs, size_t sz); + static bool readSafeSize(std::ifstream& fs, size_t& sz); + + static bool writeVec1D(std::ofstream& fs, const dVec1& v); + static bool writeIVec1D(std::ofstream& fs, const iVec1& v); + static bool writeVec2D(std::ofstream& fs, const dVec2& v); + static bool writeIVec2D(std::ofstream& fs, const iVec2& v); + static bool writeVec3D(std::ofstream& fs, const dVec3& v); + + static bool readVec1D(std::ifstream& fs, dVec1& v); + static bool readIVec1D(std::ifstream& fs, iVec1& v); + static bool readVec2D(std::ifstream& fs, dVec2& v); + static bool readIVec2D(std::ifstream& fs, iVec2& v); + static bool readVec3D(std::ifstream& fs, dVec3& v); + + // gridOutput2 序列化 + static bool writeGridOutput2(std::ofstream& fs, const HX_NWTM_GRID_OUTPUT2& g); + static bool readGridOutput2(std::ifstream& fs, HX_NWTM_GRID_OUTPUT2& g); + + static bool readHeader(std::ifstream& fs, GridCacheHeader& hdr); + static bool writeHeader(std::ofstream& fs, const GridCacheHeader& hdr); + + // “清洗”不稳定字段：用 vector size 统一覆盖 n，避免 DLL 未初始化导致的随机值 + static void sanitizeGridOutput2(HX_NWTM_GRID_OUTPUT2& g); +}; + +#endif // GRID_CACHE_IO_H diff --git a/ML/Training/Training/SceneIO.cpp b/ML/Training/Training/SceneIO.cpp new file mode 100644 index 0000000..e7b8a7d --- /dev/null +++ b/ML/Training/Training/SceneIO.cpp @@ -0,0 +1,226 @@ +#include "SceneIO.h" +#include +#include + +bool SceneIO::readInt(std::ifstream& fs, int& value) +{ + fs.read(reinterpret_cast(&value), sizeof(int)); + return fs.good(); +} + +bool SceneIO::readDouble(std::ifstream& fs, double& value) +{ + fs.read(reinterpret_cast(&value), sizeof(double)); + return fs.good(); +} + +bool SceneIO::readSafeSize(std::ifstream& fs, size_t& size) +{ + uint32_t sz = 0; + fs.read(reinterpret_cast(&sz), sizeof(uint32_t)); + if (!fs.good()) return false; + + const uint32_t MAX_SIZE = 100000000; + if (sz > MAX_SIZE) return false; + + size = static_cast(sz); + return true; +} + +bool SceneIO::readString(std::ifstream& fs, std::string& value) +{ + uint32_t byteLen = 0; + fs.read(reinterpret_cast(&byteLen), sizeof(uint32_t)); + if (!fs.good()) return false; + + // 0xFFFFFFFF 表示 null 字符串 + if (byteLen == 0xFFFFFFFF) { + value.clear(); + return true; + } + + // UTF-16 字节数 -> 字符数 + size_t charCount = (size_t)(byteLen / 2); + + std::vector utf16(charCount); + if (charCount > 0) { + fs.read(reinterpret_cast(&utf16[0]), byteLen); + if (!fs.good()) return false; + } + + // 简化转换：只取低字节（你的数据目前可用） + value.clear(); + value.reserve(charCount); + for (size_t i = 0; i < charCount; ++i) { + value += static_cast(utf16[i] & 0xFF); + } + + return true; +} + +bool SceneIO::readVector1D(std::ifstream& fs, std::vector& vec) +{ + size_t sz = 0; + if (!readSafeSize(fs, sz)) return false; + + vec.resize(sz); + if (sz > 0) { + fs.read(reinterpret_cast(&vec[0]), sz * sizeof(double)); + if (!fs.good()) return false; + } + return true; +} + +bool SceneIO::readVector2D(std::ifstream& fs, std::vector >& vec) +{ + size_t rows = 0; + if (!readSafeSize(fs, rows)) return false; + + vec.resize(rows); + for (size_t i = 0; i < rows; ++i) { + if (!readVector1D(fs, vec[i])) return false; + } + return true; +} + +bool SceneIO::readVector3D(std::ifstream& fs, std::vector > >& vec) +{ + size_t depth = 0; + if (!readSafeSize(fs, depth)) return false; + + vec.resize(depth); + for (size_t i = 0; i < depth; ++i) { + if (!readVector2D(fs, vec[i])) return false; + } + return true; +} + +bool SceneIO::readStdVecI(std::ifstream& fs, std::vector& vec) +{ + size_t sz = 0; + if (!readSafeSize(fs, sz)) return false; + + vec.resize(sz); + if (sz > 0) { + for (size_t i = 0; i < sz; ++i) { + if (!readInt(fs, vec[i])) return false; + } + } + return true; +} + +bool SceneIO::loadScene(const std::string& filename, PebiScene& scene) +{ + std::ifstream fs(filename.c_str(), std::ios::binary); + if (!fs) { + std::cerr << "无法打开 scene 文件: " << filename << std::endl; + return false; + } + + // magic + uint32_t magic = 0; + fs.read(reinterpret_cast(&magic), sizeof(uint32_t)); + if (!fs.good() || magic != 0x4E4D5343) { // 'CSMN' + std::cerr << "scene.bin magic 错误" << std::endl; + return false; + } + + if (!readInt(fs, scene.version)) return false; + + // 网格基础 + if (!readInt(fs, scene.D)) return false; + if (!readDouble(fs, scene.GridControl)) return false; + + // 几何 + if (!readVector2D(fs, scene.Boundary)) return false; + if (!readVector2D(fs, scene.VerticalWell)) return false; + if (!readVector2D(fs, scene.HorizontalWell)) return false; + if (!readVector2D(fs, scene.FractureVerticalWell)) return false; + if (!readVector3D(fs, scene.MultistageFracturedHorizontalWell)) return false; + if (!readVector2D(fs, scene.InclinedWell)) return false; + if (!readVector2D(fs, scene.Fault)) return false; + + // wellType + size_t nType = 0; + if (!readSafeSize(fs, nType)) return false; + scene.wellType.resize(nType); + for (size_t i = 0; i < nType; ++i) { + if (!readInt(fs, scene.wellType[i])) return false; + } + + // wellName + size_t nName = 0; + if (!readSafeSize(fs, nName)) return false; + scene.wellName.resize(nName); + for (size_t i = 0; i < nName; ++i) { + if (!readString(fs, scene.wellName[i])) return false; + } + + // solverType + if (!readInt(fs, scene.solverType)) return false; + + // Rate + if (!readVector2D(fs, scene.Rate.t)) return false; + if (!readVector2D(fs, scene.Rate.qo)) return false; + if (!readVector2D(fs, scene.Rate.qg)) return false; + if (!readVector2D(fs, scene.Rate.qw)) return false; + + // CS + if (!readVector1D(fs, scene.CS.C)) return false; + if (!readVector1D(fs, scene.CS.S)) return false; + + // WellFlowSectionIndex + if (!readStdVecI(fs, scene.wellFlowSectionIndex)) return false; + + // PVT + if (!readVector1D(fs, scene.PVT.p)) return false; + if (!readDouble(fs, scene.PVT.pb)) return false; + + if (!readVector1D(fs, scene.PVT.Rso)) return false; + if (!readVector1D(fs, scene.PVT.Bo)) return false; + if (!readVector1D(fs, scene.PVT.Co)) return false; + if (!readVector1D(fs, scene.PVT.miuo)) return false; + if (!readVector1D(fs, scene.PVT.rouo)) return false; + + if (!readVector1D(fs, scene.PVT.Rv)) return false; + if (!readVector1D(fs, scene.PVT.Bg)) return false; + if (!readVector1D(fs, scene.PVT.Cg)) return false; + if (!readVector1D(fs, scene.PVT.miug)) return false; + if (!readVector1D(fs, scene.PVT.roug)) return false; + if (!readVector1D(fs, scene.PVT.Z)) return false; + + if (!readVector1D(fs, scene.PVT.Rsw)) return false; + if (!readVector1D(fs, scene.PVT.Bw)) return false; + if (!readVector1D(fs, scene.PVT.Cw)) return false; + if (!readVector1D(fs, scene.PVT.miuw)) return false; + if (!readVector1D(fs, scene.PVT.rouw)) return false; + + if (!readVector1D(fs, scene.PVT.V)) return false; + if (!readVector1D(fs, scene.PVT.k_kinitial)) return false; + if (!readVector1D(fs, scene.PVT.Cf_Cfinitial)) return false; + + if (!readVector1D(fs, scene.PVT.So)) return false; + if (!readVector1D(fs, scene.PVT.Kro)) return false; + if (!readVector1D(fs, scene.PVT.Sg)) return false; + if (!readVector1D(fs, scene.PVT.Krg)) return false; + if (!readVector1D(fs, scene.PVT.Sw)) return false; + if (!readVector1D(fs, scene.PVT.Krw)) return false; + + // Base + if (!readDouble(fs, scene.Base.Pi)) return false; + if (!readDouble(fs, scene.Base.Cti)) return false; + if (!readDouble(fs, scene.Base.Cf)) return false; + if (!readDouble(fs, scene.Base.Soi)) return false; + if (!readDouble(fs, scene.Base.Sgi)) return false; + if (!readDouble(fs, scene.Base.Swi)) return false; + + if (!readDouble(fs, scene.Base.d)) return false; + if (!readDouble(fs, scene.Base.dt_Min)) return false; + if (!readDouble(fs, scene.Base.dt_Max)) return false; + + if (!readDouble(fs, scene.Base.k_ref)) return false; + if (!readDouble(fs, scene.Base.phi_ref)) return false; + if (!readDouble(fs, scene.Base.h_ref)) return false; + + return true; +} \ No newline at end of file diff --git a/ML/Training/Training/SceneIO.h b/ML/Training/Training/SceneIO.h new file mode 100644 index 0000000..35fdbd0 --- /dev/null +++ b/ML/Training/Training/SceneIO.h @@ -0,0 +1,84 @@ +#pragma once +#include +#include +#include + +// 场景结构体 +struct PebiScene { + int version; + + // 网格几何 + int D; + double GridControl; + std::vector > Boundary; + std::vector > VerticalWell; + std::vector > HorizontalWell; + std::vector > FractureVerticalWell; + std::vector > > MultistageFracturedHorizontalWell; + std::vector > InclinedWell; + std::vector > Fault; + + // 井顺序 + std::vector wellType; + std::vector wellName; + + // 求解器参数 + int solverType; + + // Rate 流量数据 + struct { + std::vector > t; + std::vector > qo; + std::vector > qg; + std::vector > qw; + } Rate; + + // CS 井筒参数 + struct { + std::vector C; + std::vector S; + } CS; + + // 每口井的流量段索引 + std::vector wellFlowSectionIndex; + + // PVT 数据 + struct { + std::vector p; + double pb; + std::vector Rso, Bo, Co, miuo, rouo; + std::vector Rv, Bg, Cg, miug, roug, Z; + std::vector Rsw, Bw, Cw, miuw, rouw; + std::vector V, k_kinitial, Cf_Cfinitial; + std::vector So, Kro, Sg, Krg, Sw, Krw; + } PVT; + + // Base 储层参数 + struct { + double Pi, Cti, Cf, Soi, Sgi, Swi; + double d, dt_Min, dt_Max; + double k_ref, phi_ref, h_ref; + } Base; +}; + +// scene.bin 读取器 +class SceneIO { +public: + static bool loadScene(const std::string& filename, PebiScene& scene); + +private: + static bool readInt(std::ifstream& fs, int& value); + static bool readDouble(std::ifstream& fs, double& value); + + // 读取 UTF-16 字符串（scene.bin 的格式） + static bool readString(std::ifstream& fs, std::string& value); + + static bool readSafeSize(std::ifstream& fs, size_t& size); + + static bool readVector1D(std::ifstream& fs, std::vector& vec); + static bool readVector2D(std::ifstream& fs, std::vector >& vec); + static bool readVector3D(std::ifstream& fs, std::vector > >& vec); + + // 读取整数向量 + static bool readStdVecI(std::ifstream& fs, std::vector& vec); +}; diff --git a/ML/Training/Training/Training.cpp b/ML/Training/Training/Training.cpp new file mode 100644 index 0000000..59ee8b8 --- /dev/null +++ b/ML/Training/Training/Training.cpp @@ -0,0 +1,318 @@ +#include +#include +#include +#include + +#include "pch.h" +#include "SceneIO.h" +#include "GridCacheIO.h" +#include "DatasetIO.h" + +typedef void (*HX_NWTM_GRID_Func)(HX_NWTM_GRID_OUTPUT1&, HX_NWTM_GRID_OUTPUT2&, + const HX_NWTM_GRID_INPUT&, std::string); +typedef void (*HX_NWTM_MODEL_Func)(HX_NWTM_MODEL_OUTPUT&, const HX_NWTM_MODEL_INPUT&, std::string); + +// 获取目录 +static std::string getExeDir() +{ + char buf[MAX_PATH] = {0}; + DWORD len = GetModuleFileNameA(NULL, buf, MAX_PATH); + if (len == 0) { + return "."; + } + std::string path(buf, len); + size_t pos = path.find_last_of("\\/"); + if (pos != std::string::npos) { + path.resize(pos); + } + return path; +} + +// 从 scene 重建 GRID 输入（与之前一致） +static HX_NWTM_GRID_INPUT rebuildGridInput(const PebiScene& scene) +{ + HX_NWTM_GRID_INPUT input; + input.D = scene.D; + input.GridControl = scene.GridControl; + input.Boundary = scene.Boundary; + input.VerticalWell = scene.VerticalWell; + input.HorizontalWell = scene.HorizontalWell; + input.FractureVerticalWell = scene.FractureVerticalWell; + input.MultistageFracturedHorizontalWell = scene.MultistageFracturedHorizontalWell; + input.InclinedWell = scene.InclinedWell; + input.Fault = scene.Fault; + return input; +} + +// 从 scene + gridOutput2 重建 MODEL 输入（与之前一致） +static HX_NWTM_MODEL_INPUT rebuildModelInput(const PebiScene& scene, const HX_NWTM_GRID_OUTPUT2& gridOutput) +{ + HX_NWTM_MODEL_INPUT input(gridOutput); + + input.T = scene.solverType; + + input.Rate.t = scene.Rate.t; + input.Rate.qo = scene.Rate.qo; + input.Rate.qg = scene.Rate.qg; + input.Rate.qw = scene.Rate.qw; + + input.CS.C = scene.CS.C; + input.CS.S = scene.CS.S; + + input.PVT.p = scene.PVT.p; + input.PVT.pb = scene.PVT.pb; + input.PVT.Rso = scene.PVT.Rso; + input.PVT.Bo = scene.PVT.Bo; + input.PVT.Co = scene.PVT.Co; + input.PVT.miuo = scene.PVT.miuo; + input.PVT.rouo = scene.PVT.rouo; + input.PVT.Rv = scene.PVT.Rv; + input.PVT.Bg = scene.PVT.Bg; + input.PVT.Cg = scene.PVT.Cg; + input.PVT.miug = scene.PVT.miug; + input.PVT.roug = scene.PVT.roug; + input.PVT.Z = scene.PVT.Z; + input.PVT.Rsw = scene.PVT.Rsw; + input.PVT.Bw = scene.PVT.Bw; + input.PVT.Cw = scene.PVT.Cw; + input.PVT.miuw = scene.PVT.miuw; + input.PVT.rouw = scene.PVT.rouw; + input.PVT.V = scene.PVT.V; + input.PVT.k_kinitial = scene.PVT.k_kinitial; + input.PVT.Cf_Cfinitial = scene.PVT.Cf_Cfinitial; + input.PVT.So = scene.PVT.So; + input.PVT.Kro = scene.PVT.Kro; + input.PVT.Sg = scene.PVT.Sg; + input.PVT.Krg = scene.PVT.Krg; + input.PVT.Sw = scene.PVT.Sw; + input.PVT.Krw = scene.PVT.Krw; + + input.Base.Pi = scene.Base.Pi; + input.Base.Cti = scene.Base.Cti; + input.Base.Cf = scene.Base.Cf; + input.Base.Soi = scene.Base.Soi; + input.Base.Sgi = scene.Base.Sgi; + input.Base.Swi = scene.Base.Swi; + input.Base.d = scene.Base.d; + input.Base.dt_Min = scene.Base.dt_Min; + input.Base.dt_Max = scene.Base.dt_Max; + + // 用参考值填满每个网格单元（最小可用） + size_t nCells = gridOutput.Trinodexy.size(); + input.Base.k = dVec1(nCells, scene.Base.k_ref); + input.Base.phi = dVec1(nCells, scene.Base.phi_ref); + input.Base.h = dVec1(nCells, scene.Base.h_ref); + + return input; +} + +// 导出验证样本（CSV） +static void exportValidationSample(const HX_NWTM_MODEL_OUTPUT& output, const std::string& filename) +{ + std::ofstream f(filename.c_str()); + f << "t,pw\n"; + if (!output.pw.empty()) { + for (size_t i = 0; i < output.t.size() && i < output.pw[0].size(); ++i) { + f << output.t[i] << "," << output.pw[0][i] << "\n"; + } + } + f.close(); + std::cout << "已导出: " << filename << std::endl; +} + +// 打印输入摘要：用于快速确认“为什么 steps 变了” +static void printInputSummary(const PebiScene& scene) +{ + size_t wells = scene.Rate.t.size(); + size_t tlen0 = (wells > 0 ? scene.Rate.t[0].size() : 0); + + std::cout << "输入摘要：wells=" << (unsigned int)wells + << ", solver=" << scene.solverType + << ", Rate.t[0].len=" << (unsigned int)tlen0 + << ", PVT点=" << (unsigned int)scene.PVT.p.size() + << std::endl; +} + +int main() +{ + std::cout << "=== PEBI Scene Processor (Dataset First) ===" << std::endl; + + // 按需开关：发布版保持 true（因为 dataset 是后续批量采样的核心） + const bool ENABLE_DATASET_FIRST = true; + const bool ENABLE_DATASET_WRITE = true; // 当走 scene+grid 生成时，是否写 dataset.bin + + // 路径 + std::string exeDir = getExeDir(); + + // license 仍然在工程根目录下的 Bin\Res\license + std::string lic = exeDir + "\\..\\..\\..\\Bin\\Res\\license\\HXNWTM_license.dat"; + + // 所有数据统一放到 ML/nmWTAI-ML/data/temp + //std::string dataDir = exeDir + "\\..\\..\\Data"; + //std::string scenePath = dataDir + "\\scene.bin"; + //std::string datasetPath = dataDir + "\\dataset.bin"; + //std::string gridCachePath= dataDir + "\\grid_cache.bin"; + std::string dataDir = exeDir + "\\..\\..\\nmWTAI-ML\\data\\temp"; + std::string scenePath = dataDir + "\\scene.bin"; + std::string datasetPath = dataDir + "\\dataset.bin"; + std::string gridCachePath= dataDir + "\\grid_cache.bin"; + + // 1) 加载 DLL + HMODULE hx = LoadLibraryW(L"HX_NWTM.dll"); + if (!hx) { + std::cerr << "LoadLibrary failed, err=" << GetLastError() << std::endl; + return 1; + } + + HX_NWTM_GRID_Func HX_NWTM_GRID = (HX_NWTM_GRID_Func)GetProcAddress(hx, "HX_NWTM_GRID"); + HX_NWTM_MODEL_Func HX_NWTM_MODEL = (HX_NWTM_MODEL_Func)GetProcAddress(hx, "HX_NWTM_MODEL"); + if (!HX_NWTM_GRID || !HX_NWTM_MODEL) { + std::cerr << "GetProcAddress failed, err=" << GetLastError() << std::endl; + FreeLibrary(hx); + return 2; + } + + // 2) 检查 license + DWORD attr = GetFileAttributesA(lic.c_str()); + if (attr == INVALID_FILE_ATTRIBUTES) { + std::cerr << "License not found: " << lic << std::endl; + FreeLibrary(hx); + return 3; + } + std::cout << "License OK: " << lic << std::endl; + + // 3) dataset 优先尝试（但如果 scene.bin 存在，会做 sceneKey 串场景校验） + PebiScene scene; + HX_NWTM_GRID_OUTPUT2 gridOutput2; + bool gotFromDataset = false; + + if (ENABLE_DATASET_FIRST && DatasetIO::fileExists(datasetPath)) { + std::cout << "\n发现 dataset: " << datasetPath << "，尝试直接加载..." << std::endl; + + // 先读 dataset header + unsigned long long dsKey = 0; + unsigned int dsCells = 0; + unsigned int dsWells = 0; + int dsSolver = 0; + + bool headerOK = DatasetIO::readDatasetSceneKey(datasetPath, dsKey, dsCells, dsWells, dsSolver); + + // 如果 scene.bin 存在，我们就顺便做一次“串场景保护”：sceneKey 不匹配就不用 dataset + bool sceneKeyMatch = true; + + if (headerOK && DatasetIO::fileExists(scenePath)) { + std::cout << "\n加载 scene.bin（仅用于串场景校验）..." << std::endl; + if (SceneIO::loadScene(scenePath, scene)) { + HX_NWTM_GRID_INPUT gi = rebuildGridInput(scene); + unsigned long long expected = DatasetIO::computeSceneKey64(gi); + if (expected != dsKey) { + sceneKeyMatch = false; + std::cout << "dataset sceneKey 不匹配（串场景保护触发），忽略 dataset，准备重建..." << std::endl; + } + } else { + std::cout << "scene.bin 加载失败，跳过串场景校验，直接尝试 dataset 自检..." << std::endl; + sceneKeyMatch = true; + } + } + + if (sceneKeyMatch) { + // 真正加载 dataset（strictVerify=true 会用 dataset 内部 scene 重新算 key） + PebiScene dsScene; + HX_NWTM_GRID_OUTPUT2 dsGrid; + if (DatasetIO::loadDataset(datasetPath, dsScene, dsGrid, true)) { + scene = dsScene; + gridOutput2 = dsGrid; + gotFromDataset = true; + + std::cout << "dataset 加载成功：cells=" << (unsigned int)gridOutput2.Trinodexy.size() + << ", wells=" << (unsigned int)(scene.Rate.t.size()) + << ", solver=" << scene.solverType << std::endl; + } else { + std::cout << "dataset 加载失败，准备走 scene+grid_cache..." << std::endl; + } + } + } + + // 4) 如果没有拿到 dataset，就走 scene.bin + grid_cache.bin（带串场景保护） + if (!gotFromDataset) { + std::cout << "\n加载 scene.bin..." << std::endl; + if (!SceneIO::loadScene(scenePath, scene)) { + std::cerr << "scene.bin 读取失败: " << scenePath << std::endl; + FreeLibrary(hx); + return 4; + } + + // 输入摘要：帮助你定位 “为什么 steps 变了” + std::cout << "scene.bin 读取成功：Wells=" << (unsigned int)scene.wellName.size() + << ", PVT点=" << (unsigned int)scene.PVT.p.size() + << ", Solver=" << scene.solverType << std::endl; + printInputSummary(scene); + + HX_NWTM_GRID_INPUT gridInput = rebuildGridInput(scene); + + HX_NWTM_GRID_OUTPUT1 gridOutput1; + bool loadedFromCache = false; + + if (GridCacheIO::fileExists(gridCachePath)) { + std::cout << "\n读取网格缓存: " << gridCachePath << std::endl; + if (GridCacheIO::loadGrid(gridCachePath, gridOutput2, gridInput)) { + loadedFromCache = true; + } else { + std::cout << "缓存无效或串场景不匹配，准备重新生成..." << std::endl; + } + } + + if (!loadedFromCache) { + std::cout << "\n生成网格..." << std::endl; + try { + HX_NWTM_GRID(gridOutput1, gridOutput2, gridInput, lic); + } catch (...) { + std::cerr << "HX_NWTM_GRID 异常" << std::endl; + FreeLibrary(hx); + return 5; + } + + std::cout << "保存网格缓存: " << gridCachePath << std::endl; + if (!GridCacheIO::saveGrid(gridCachePath, gridOutput2, gridInput)) { + std::cout << "警告：保存网格缓存失败（继续执行）" << std::endl; + } + } + + std::cout << "网格就绪: " << (unsigned int)gridOutput2.Trinodexy.size() + << " cells" << (loadedFromCache ? " (from cache)" : " (generated)") << std::endl; + + // 生成/更新 dataset.bin（可开关） + if (ENABLE_DATASET_WRITE) { + std::cout << "\n写入 dataset: " << datasetPath << std::endl; + if (!DatasetIO::saveDataset(datasetPath, scene, gridInput, gridOutput2)) { + std::cout << "警告：dataset 写入失败（继续执行）" << std::endl; + } else { + std::cout << "dataset 写入成功" << std::endl; + } + } + } + + //// 5) 求解模型（这一步不是“多余验证”，而是 dataset 生成/采样最终都要用到的核心步骤） + //std::cout << "\n开始求解模型..." << std::endl; + + //HX_NWTM_MODEL_INPUT modelInput = rebuildModelInput(scene, gridOutput2); + //HX_NWTM_MODEL_OUTPUT modelOutput; + + //try { + // HX_NWTM_MODEL(modelOutput, modelInput, lic); + //} catch (...) { + // std::cerr << "HX_NWTM_MODEL 异常" << std::endl; + // FreeLibrary(hx); + // return 6; + //} + + //std::cout << "===================================\n"; + //std::cout << "求解完成: t=" << (unsigned int)modelOutput.t.size() + // << " steps, wells=" << (unsigned int)modelOutput.pw.size() << std::endl; + + //exportValidationSample(modelOutput, dataDir + "\\validation.csv"); + + FreeLibrary(hx); + std::cout << "\nDone!" << std::endl; + return 0; +} diff --git a/ML/Training/include/framework.h b/ML/Training/include/framework.h new file mode 100644 index 0000000..2529f7a --- /dev/null +++ b/ML/Training/include/framework.h @@ -0,0 +1,4 @@ +#pragma once +#define WIN32_LEAN_AND_MEAN // 从 Windows 头文件中排除极少使用的内容 +// Windows 头文件 +#include \ No newline at end of file diff --git a/ML/Training/include/netgenmesh.h b/ML/Training/include/netgenmesh.h new file mode 100644 index 0000000..0df975d --- /dev/null +++ b/ML/Training/include/netgenmesh.h @@ -0,0 +1,167 @@ +// netgen网格生成库的C++接口，传入边界、井筒、裂缝、断层和复合区，生成二维三角形网格。 +// 作者：万义钊 +// 日期：2025年1月28日，除夕 +// 版本：1.0 +#pragma once +namespace nglib +{ +#include "nglib.h" +} +using namespace nglib; +#include +#include +#include +enum BOUND_TYPE +{ + CONST_PRESSURE, + CONST_RATE, + GEO_LIMIT +}; +enum BOUND_SHAPE +{ + CIRCLE, + RECTANGLE, + POLYGON +}; +enum WELL_TYPE +{ + VERTICAL, + VERTICAL_FRACTURED, + MULTI_FRACED_HORIZONTAL +}; +struct Point +{ + double x, y, z; + std::vector 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; + } +}; + +struct Cell +{ + int type; + std::vector pointIndices; + std::vector cellData; + Cell(const Cell &other) + { + type = other.type; + pointIndices = other.pointIndices; + cellData = other.cellData; + } + Cell() {} +}; +class Line +{ + +public: + Point p1; + Point p2; +}; + +// 裂缝线 +class Frac : public Line +{ +public: + Frac(Point p1, Point p2) + { + this->p1 = p1; + this->p2 = p2; + } + Frac() {} + // 交点 + int nIntersect; // 交点数量 + std::vector vector_interPoint; // 交点坐标 +}; + +// 定义外边界的线 +class Segment : public Line +{ +public: + Segment(Point p1, Point p2) + { + this->p1 = p1; + this->p2 = p2; + } + Segment() {} + BOUND_TYPE boundType; // 边界线的边界条件类型 + int nNumNodes; // 边界线上的网格点数 +}; + +// 外边界的基类 +class CBoundShape +{ +public: + BOUND_SHAPE boundShape; // 定义外边界类型，0为圆形，1为四边形，2为多边形 + + // 表示多边形（包括四边形）的参数 + int nNumSegs; // 边的数量 + std::vector vecSegments; // 边 + + // 表示圆形 + Point cCenter; + double dRadius; + int nNodeNum; + BOUND_TYPE boundType; +}; + +// 多边形,内部的复合区域 +class CBoundPolygon : public CBoundShape +{ +public: + int nNumSegs; // 边的数量 + std::vector vecSegments; // 边 +}; + +// 表示井筒圆形边界 +class CBoundWell +{ +public: + CBoundWell(){wellType=VERTICAL;}//默认井为直井 +public: + + WELL_TYPE wellType; + + //直井的几何参数，wellType=VERTICAL + Point cCenter; + double dRadius; + int nNodeNum; + BOUND_TYPE boundType; + + //直井压裂井独有的参数，wellType=VERTICAL_FRACTURED + //直井的中心和半径仍然使用cCenter和dRaduis表示 + double dHf;//裂缝半长 + double dTheata;//裂缝角度(相对于x轴) + + //多段压裂水平井的参数，cCenter表示水平井中心点坐标，wellType=MULTI_FRACED_HORIZONTAL + double dLength;//长度 + double dBeta;//角度 + int nFracNum;//压裂裂缝的数量 + std::vector vecFracs;//压裂裂缝 + +}; + +#ifdef __cplusplus +extern "C" +{ +#endif + __declspec(dllexport) Ng_Mesh *NgMeshGen2D(std::string strFileIn2D); + __declspec(dllexport) Ng_Mesh *NgMeshGen2DByIn2d(std::vector wells, std::vector limits, std::vector faults, std::vector fracs, CBoundShape shape); + __declspec(dllexport) void Ng2VTK(Ng_Mesh *mesh); + +#ifdef __cplusplus +} +#endif \ No newline at end of file diff --git a/ML/Training/include/nginterface.h b/ML/Training/include/nginterface.h new file mode 100644 index 0000000..0477035 --- /dev/null +++ b/ML/Training/include/nginterface.h @@ -0,0 +1,476 @@ +#ifndef NGINTERFACE +#define NGINTERFACE + + + + + +/**************************************************************************/ +/* File: nginterface.h */ +/* Author: Joachim Schoeberl */ +/* Date: 20. Nov. 99 */ +/**************************************************************************/ + +/* + Application program interface to Netgen + +*/ + +#ifdef WIN32 + #if NGINTERFACE_EXPORTS || NGLIB_EXPORTS || nglib_EXPORTS + #define DLL_HEADER __declspec(dllexport) + #else + #define DLL_HEADER __declspec(dllimport) + #endif +#else + #define DLL_HEADER +#endif + + +// max number of nodes per element +#define NG_ELEMENT_MAXPOINTS 12 + +// max number of nodes per surface element +#define NG_SURFACE_ELEMENT_MAXPOINTS 8 + + + +// implemented element types: +enum NG_ELEMENT_TYPE { + NG_PNT = 0, + NG_SEGM = 1, NG_SEGM3 = 2, + NG_TRIG = 10, NG_QUAD=11, NG_TRIG6 = 12, NG_QUAD6 = 13, + NG_TET = 20, NG_TET10 = 21, + NG_PYRAMID = 22, NG_PRISM = 23, NG_PRISM12 = 24, + NG_HEX = 25 +}; + +typedef double NG_POINT[3]; // coordinates +typedef int NG_EDGE[2]; // initial point, end point +typedef int NG_FACE[4]; // points, last one is 0 for trig + + +#ifdef __cplusplus +extern "C" { +#endif + + // load geomtry from file + DLL_HEADER void Ng_LoadGeometry (const char * filename); + + // load netgen mesh + DLL_HEADER void Ng_LoadMesh (const char * filename); + + // load netgen mesh + DLL_HEADER void Ng_LoadMeshFromString (const char * mesh_as_string); + + // space dimension (2 or 3) + DLL_HEADER int Ng_GetDimension (); + + // number of mesh points + DLL_HEADER int Ng_GetNP (); + + // number of mesh vertices (differs from GetNP for 2nd order elements) + DLL_HEADER int Ng_GetNV (); + + // number of mesh elements + DLL_HEADER int Ng_GetNE (); + + // number of surface triangles + DLL_HEADER int Ng_GetNSE (); + + // Get Point coordintes, index from 1 .. np + DLL_HEADER void Ng_GetPoint (int pi, double * p); + + // Get Element Points + DLL_HEADER NG_ELEMENT_TYPE Ng_GetElement (int ei, int * epi, int * np = 0); + + // Get Element Type + DLL_HEADER NG_ELEMENT_TYPE Ng_GetElementType (int ei); + + // Get sub-domain of element ei + DLL_HEADER int Ng_GetElementIndex (int ei); + + DLL_HEADER void Ng_SetElementIndex(const int ei, const int index); + + // Get Material of element ei + DLL_HEADER char * Ng_GetElementMaterial (int ei); + + // Get Material of domain dom + DLL_HEADER char * Ng_GetDomainMaterial (int dom); + + // Get User Data + DLL_HEADER int Ng_GetUserDataSize (char * id); + DLL_HEADER void Ng_GetUserData (char * id, double * data); + + // Get Surface Element Points + DLL_HEADER NG_ELEMENT_TYPE Ng_GetSurfaceElement (int ei, int * epi, int * np = 0); + + // Get Surface Element Type + DLL_HEADER NG_ELEMENT_TYPE Ng_GetSurfaceElementType (int ei); + + // Get Surface Element Index + DLL_HEADER int Ng_GetSurfaceElementIndex (int ei); + + // Get Surface Element Surface Number + DLL_HEADER int Ng_GetSurfaceElementSurfaceNumber (int ei); + + // Get Surface Element Number + DLL_HEADER int Ng_GetSurfaceElementFDNumber (int ei); + + // Get BCName for Surface Element + DLL_HEADER char * Ng_GetSurfaceElementBCName (int ei); + //void Ng_GetSurfaceElementBCName (int ei, char * name); + + // Get BCName for bc-number + DLL_HEADER char * Ng_GetBCNumBCName (int bcnr); + //void Ng_GetBCNumBCName (int bcnr, char * name); + + // Get normal vector of surface element node + DLL_HEADER void Ng_GetNormalVector (int sei, int locpi, double * nv); + + + DLL_HEADER void Ng_SetPointSearchStartElement(int el); + + // Find element of point, returns local coordinates + DLL_HEADER int Ng_FindElementOfPoint (double * p, double * lami, + int build_searchtrees = 0, + const int * const indices = NULL, const int numind = 0); + + // Find surface element of point, returns local coordinates + DLL_HEADER int Ng_FindSurfaceElementOfPoint (double * p, double * lami, + int build_searchtrees = 0, + const int * const indices = NULL, const int numind = 0); + + + // is elment ei curved ? + DLL_HEADER int Ng_IsElementCurved (int ei); + // is elment sei curved ? + DLL_HEADER int Ng_IsSurfaceElementCurved (int sei); + + /// Curved Elemens: + /// xi..local coordinates + /// x ..global coordinates + /// dxdxi...D x D Jacobian matrix (row major storage) + DLL_HEADER void Ng_GetElementTransformation (int ei, const double * xi, + double * x, double * dxdxi); + + + /// buffer must be at least 100 doubles, alignment of double + DLL_HEADER void Ng_GetBufferedElementTransformation (int ei, const double * xi, + double * x, double * dxdxi, + void * buffer, int buffervalid); + + + + /// Curved Elemens: + /// xi..local coordinates + /// x ..global coordinates + /// dxdxi...D x D-1 Jacobian matrix (row major storage) + /// curved ...is element curved ? + DLL_HEADER void Ng_GetSurfaceElementTransformation (int sei, const double * xi, + double * x, double * dxdxi); + + /// Curved Elemens: + /// xi..local coordinates + /// sxi..step xi + /// x ..global coordinates + /// dxdxi...D x D Jacobian matrix (row major storage) + DLL_HEADER void Ng_GetMultiElementTransformation (int ei, int n, + const double * xi, size_t sxi, + double * x, size_t sx, + double * dxdxi, size_t sdxdxi); + + + + DLL_HEADER int Ng_GetSegmentIndex (int elnr); + DLL_HEADER NG_ELEMENT_TYPE Ng_GetSegment (int elnr, int * epi, int * np = 0); + + + + + // Mark element for refinement + DLL_HEADER void Ng_SetRefinementFlag (int ei, int flag); + DLL_HEADER void Ng_SetSurfaceRefinementFlag (int sei, int flag); + + // Do local refinement + enum NG_REFINEMENT_TYPE { NG_REFINE_H = 0, NG_REFINE_P = 1, NG_REFINE_HP = 2 }; + DLL_HEADER void Ng_Refine (NG_REFINEMENT_TYPE reftype); + + // Use second order elements + DLL_HEADER void Ng_SecondOrder (); + DLL_HEADER void Ng_HighOrder (int order, bool rational = false); + //void Ng_HPRefinement (int levels, double parameter = 0.125); + DLL_HEADER void Ng_HPRefinement (int levels, double parameter = 0.125, + bool setorders = true,bool ref_level = false); + // void Ng_HPRefinement (int levels); + // void Ng_HPRefinement (int levels, double parameter); + + + // Topology and coordinate information of master element: + + DLL_HEADER int Ng_ME_GetNVertices (NG_ELEMENT_TYPE et); + DLL_HEADER int Ng_ME_GetNEdges (NG_ELEMENT_TYPE et); + DLL_HEADER int Ng_ME_GetNFaces (NG_ELEMENT_TYPE et); + + DLL_HEADER const NG_POINT * Ng_ME_GetVertices (NG_ELEMENT_TYPE et); + DLL_HEADER const NG_EDGE * Ng_ME_GetEdges (NG_ELEMENT_TYPE et); + DLL_HEADER const NG_FACE * Ng_ME_GetFaces (NG_ELEMENT_TYPE et); + + DLL_HEADER void Ng_UpdateTopology(); + + DLL_HEADER int Ng_GetNEdges(); + DLL_HEADER int Ng_GetNFaces(); + + + DLL_HEADER int Ng_GetElement_Edges (int elnr, int * edges, int * orient = 0); + DLL_HEADER int Ng_GetElement_Faces (int elnr, int * faces, int * orient = 0); + + DLL_HEADER int Ng_GetSurfaceElement_Edges (int selnr, int * edges, int * orient = 0); + DLL_HEADER int Ng_GetSurfaceElement_Face (int selnr, int * orient = 0); + + DLL_HEADER void Ng_GetSurfaceElementNeighbouringDomains(const int selnr, int & in, int & out); + + DLL_HEADER int Ng_GetFace_Vertices (int fnr, int * vert); + DLL_HEADER void Ng_GetEdge_Vertices (int ednr, int * vert); + DLL_HEADER int Ng_GetFace_Edges (int fnr, int * edge); + + DLL_HEADER int Ng_GetNVertexElements (int vnr); + DLL_HEADER void Ng_GetVertexElements (int vnr, int * els); + + DLL_HEADER int Ng_GetElementOrder (int enr); + DLL_HEADER void Ng_GetElementOrders (int enr, int * ox, int * oy, int * oz); + + DLL_HEADER void Ng_SetElementOrder (int enr, int order); + DLL_HEADER void Ng_SetElementOrders (int enr, int ox, int oy, int oz); + + DLL_HEADER int Ng_GetSurfaceElementOrder (int enr); + DLL_HEADER void Ng_GetSurfaceElementOrders (int enr, int * ox, int * oy); + + DLL_HEADER void Ng_SetSurfaceElementOrder (int enr, int order); + DLL_HEADER void Ng_SetSurfaceElementOrders (int enr, int ox, int oy); + + // Multilevel functions: + + // number of levels: + DLL_HEADER int Ng_GetNLevels (); + // get two parent nodes (indeed vertices !) of node ni + DLL_HEADER void Ng_GetParentNodes (int ni, int * parents); + + // get parent element (first child has always same number) + DLL_HEADER int Ng_GetParentElement (int ei); + + // get parent surface element (first child has always same number) + DLL_HEADER int Ng_GetParentSElement (int ei); + + // representant of anisotropic cluster + DLL_HEADER int Ng_GetClusterRepVertex (int vi); + DLL_HEADER int Ng_GetClusterRepEdge (int edi); + DLL_HEADER int Ng_GetClusterRepFace (int fai); + DLL_HEADER int Ng_GetClusterRepElement (int eli); + + + void Ng_SurfaceElementTransformation (int eli, double x, double y, + double * p3d, double * jacobian); + +#ifdef PARALLEL + + // the folling functions are 0-base !! + + // number on distant processor + // returns pairs (dist_proc, num_on_dist_proc) + int NgPar_GetDistantNodeNums ( int nodetype, int locnum, int * pnums ); + int NgPar_GetNDistantNodeNums ( int nodetype, int locnum ); + + int NgPar_GetGlobalNodeNum (int nodetype, int locnum); + +#endif + + namespace netgen { + // #include "../visualization/soldata.hpp" + class SolutionData; + class MouseEventHandler; + class UserVisualizationObject; + } + + enum Ng_SolutionType + { NG_SOLUTION_NODAL = 1, + NG_SOLUTION_ELEMENT = 2, + NG_SOLUTION_SURFACE_ELEMENT = 3, + NG_SOLUTION_NONCONTINUOUS = 4, + NG_SOLUTION_SURFACE_NONCONTINUOUS = 5, + NG_SOLUTION_VIRTUAL_FUNCTION = 6, + NG_SOLUTION_MARKED_ELEMENTS = 10, + NG_SOLUTION_ELEMENT_ORDER = 11 + }; + + struct Ng_SolutionData + { + const char * name; // name of gridfunction + double * data; // solution values + int components; // relevant (double) components in solution vector + int dist; // # doubles per entry alignment! + int iscomplex; // complex vector ? + bool draw_surface; + bool draw_volume; + int order; // order of elements, only partially supported + Ng_SolutionType soltype; // type of solution function + netgen::SolutionData * solclass; + }; + + // initialize solution data with default arguments + DLL_HEADER void Ng_InitSolutionData (Ng_SolutionData * soldata); + // set solution data + DLL_HEADER void Ng_SetSolutionData (Ng_SolutionData * soldata); + /// delete gridfunctions + DLL_HEADER void Ng_ClearSolutionData(); + // redraw + DLL_HEADER void Ng_Redraw(); + /// + DLL_HEADER void Ng_SetMouseEventHandler (netgen::MouseEventHandler * handler); + /// + DLL_HEADER void Ng_SetUserVisualizationObject (netgen::UserVisualizationObject * vis); + // + DLL_HEADER void Ng_SetVisualizationParameter (const char * name, + const char * value); + + + // number of periodic vertices + DLL_HEADER int Ng_GetNPeriodicVertices (int idnr); + // pairs should be an integer array of 2*npairs + DLL_HEADER void Ng_GetPeriodicVertices (int idnr, int * pairs); + + // number of periodic edges + DLL_HEADER int Ng_GetNPeriodicEdges (int idnr); + // pairs should be an integer array of 2*npairs + DLL_HEADER void Ng_GetPeriodicEdges (int idnr, int * pairs); + + DLL_HEADER void RunParallel ( void * (*fun)(void *), void * in); + + DLL_HEADER void Ng_PushStatus (const char * str); + DLL_HEADER void Ng_PopStatus (); + DLL_HEADER void Ng_SetThreadPercentage (double percent); + DLL_HEADER void Ng_GetStatus (char ** str, double & percent); + + DLL_HEADER void Ng_SetTerminate(void); + DLL_HEADER void Ng_UnSetTerminate(void); + DLL_HEADER int Ng_ShouldTerminate(void); + DLL_HEADER void Ng_SetRunning(int flag); + DLL_HEADER int Ng_IsRunning(); + + //// added by Roman Stainko .... + DLL_HEADER int Ng_GetVertex_Elements( int vnr, int* elems); + DLL_HEADER int Ng_GetVertex_SurfaceElements( int vnr, int* elems ); + DLL_HEADER int Ng_GetVertex_NElements( int vnr ); + DLL_HEADER int Ng_GetVertex_NSurfaceElements( int vnr ); + + +#ifdef SOCKETS + int Ng_SocketClientOpen( const int port, const char * host ); + void Ng_SocketClientWrite( const char * write, char ** reply); + void Ng_SocketClientClose ( void ); + void Ng_SocketClientGetServerHost ( const int number, char ** host ); + void Ng_SocketClientGetServerPort ( const int number, int * port ); + void Ng_SocketClientGetServerClientID ( const int number, int * id ); +#endif + + DLL_HEADER void Ng_InitPointCurve(double red, double green, double blue); + DLL_HEADER void Ng_AddPointCurvePoint(const double * point); + + +#ifdef PARALLEL + void Ng_SetElementPartition ( int elnr, int part ); + int Ng_GetElementPartition ( int elnr ); +#endif + + DLL_HEADER void Ng_SaveMesh ( const char * meshfile ); + DLL_HEADER void Ng_Bisect ( const char * refinementfile ); + + // if qualityloss is not equal to NULL at input, a (1-based) list of qualitylosses (due to projection) + // is saved in *qualityloss, its size is the return value + DLL_HEADER int Ng_Bisect_WithInfo ( const char * refinementfile, double ** qualityloss); + + typedef void * Ng_Mesh; + DLL_HEADER Ng_Mesh Ng_SelectMesh (Ng_Mesh mesh); + + DLL_HEADER void Ng_GetArgs (int & argc, char ** &argv); + + +#ifdef __cplusplus +} +#endif + +#endif + + + + + + +/* + The new node interface ... + it is 0-based ! + */ + +extern "C" { + + /* + number of nodes of type nt + nt = 0 is Vertex + nt = 1 is Edge + nt = 2 is Face + nt = 3 is Cell + */ + DLL_HEADER int Ng_GetNNodes (int nt); + + /* + closure nodes of node (nt, nodenr): + nodeset is bit-coded, bit 0 includes Vertices, bit 1 edges, etc + E.g., nodeset = 6 includes edge and face nodes + nodes consists of pairs of integers (nodetype, nodenr) + return value is number of nodes + */ + DLL_HEADER int Ng_GetClosureNodes (int nt, int nodenr, int nodeset, int * nodes); + + + /* + number of dim-dimensional elements + dim = 3 ... volume elements + dim = 2 ... surface elements + dim = 1 ... segments + dim = 0 ... not available + */ + DLL_HEADER int Ng_GetNElements (int dim); + + /* + closure nodes of dim-dimensional element elmentnr: + nodeset is bit-coded, bit 0 includes Vertices, bit 1 edges, etc + E.g., nodeset = 6 includes edge and face nodes + nodes consists of pairs of integers (nodetype, nodenr) + return value is number of nodes + */ + DLL_HEADER int Ng_GetElementClosureNodes (int dim, int elementnr, int nodeset, int * nodes); + + + struct Ng_Tcl_Interp; + typedef int (Ng_Tcl_CmdProc) (Ng_Tcl_Interp *interp, int argc, const char *argv[]); + + DLL_HEADER void Ng_Tcl_CreateCommand (Ng_Tcl_Interp * interp, + const char * cmdName, Ng_Tcl_CmdProc * proc); + + void Ng_Tcl_SetResult (Ng_Tcl_Interp * interp, const char * result); +} + + + + + +#ifdef __cplusplus +#include +namespace netgen +{ + DLL_HEADER extern std::ostream * testout; + DLL_HEADER extern int printmessage_importance; +} + +#endif + diff --git a/ML/Training/include/nglib.h b/ML/Training/include/nglib.h new file mode 100644 index 0000000..57e0296 --- /dev/null +++ b/ML/Training/include/nglib.h @@ -0,0 +1,747 @@ +#ifndef NGLIB +#define NGLIB + +/**************************************************************************/ +/* File: nglib.h */ +/* Author: Joachim Schoeberl */ +/* Date: 7. May. 2000 */ +/**************************************************************************/ + +/*! + \file nglib.h + \brief Library interface to the netgen meshing kernel + \author Joachim Schoeberl + \date 7. May 2000 + + This header file provides access to the core functionality of the Netgen + Mesher via a library interface, without an interactive User Interface. + + The intention of providing these set of functions is to allow system + developers to integrate Netgen into top-level code, to act as the low + level mesh generation / optimisation kernel. +*/ + +// Philippose - 14.02.2009 +// Modifications for creating a DLL in Windows +#ifdef WIN32 + #ifdef NGLIB_EXPORTS || nglib_EXPORTS + #define DLL_HEADER __declspec(dllexport) + #else + #define DLL_HEADER __declspec(dllimport) + #endif +#else + #define DLL_HEADER +#endif + + + +// ** Constants used within Netgen ********************* +/// Maximum allowed number of nodes per volume element +#define NG_VOLUME_ELEMENT_MAXPOINTS 10 + +/// Maximum allowed number of nodes per surface element +#define NG_SURFACE_ELEMENT_MAXPOINTS 8 + + + +// *** Data-types for accessing Netgen functionality *** +/// Data type for NETGEN mesh +typedef void * Ng_Mesh; + +/// Data type for NETGEN CSG geometry +typedef void * Ng_CSG_Geometry; + +/// Data type for NETGEN 2D geometry +typedef void * Ng_Geometry_2D; + +/// Data type for NETGEN STL geometry +typedef void * Ng_STL_Geometry; + +#ifdef OCCGEOMETRY +/// Data type for NETGEN OpenCascade geometry +typedef void * Ng_OCC_Geometry; +typedef void * Ng_OCC_TopTools_IndexedMapOfShape; +#endif + + +// *** Special Enum types used within Netgen *********** +/// Currently implemented surface element types +enum Ng_Surface_Element_Type + { NG_TRIG = 1, NG_QUAD = 2, NG_TRIG6 = 3, NG_QUAD6 = 4, NG_QUAD8 = 5 }; + +/// Currently implemented volume element types +enum Ng_Volume_Element_Type + { NG_TET = 1, NG_PYRAMID = 2, NG_PRISM = 3, NG_TET10 = 4 }; + +/// Values returned by Netgen functions +enum Ng_Result + { + NG_ERROR = -1, + NG_OK = 0, + NG_SURFACE_INPUT_ERROR = 1, + NG_VOLUME_FAILURE = 2, + NG_STL_INPUT_ERROR = 3, + NG_SURFACE_FAILURE = 4, + NG_FILE_NOT_FOUND = 5 + }; + + + +// *** Classes required for use within Netgen ********** +/// Netgen Meshing Parameters class +class Ng_Meshing_Parameters +{ +public: + int uselocalh; //!< Switch to enable / disable usage of local mesh size modifiers + + double maxh; //!< Maximum global mesh size allowed + double minh; //!< Minimum global mesh size allowed + + double fineness; //!< Mesh density: 0...1 (0 => coarse; 1 => fine) + double grading; //!< Mesh grading: 0...1 (0 => uniform mesh; 1 => aggressive local grading) + + double elementsperedge; //!< Number of elements to generate per edge of the geometry + double elementspercurve; //!< Elements to generate per curvature radius + + int closeedgeenable; //!< Enable / Disable mesh refinement at close edges + double closeedgefact; //!< Factor to use for refinement at close edges (larger => finer) + + int minedgelenenable; //!< Enable / Disable user defined minimum edge length for edge subdivision + double minedgelen; //!< Minimum edge length to use while subdividing the edges (default = 1e-4) + + int second_order; //!< Generate second-order surface and volume elements + int quad_dominated; //!< Creates a Quad-dominated mesh + + char * meshsize_filename; //!< Optional external mesh size file + + int optsurfmeshenable; //!< Enable / Disable automatic surface mesh optimization + int optvolmeshenable; //!< Enable / Disable automatic volume mesh optimization + + int optsteps_3d; //!< Number of optimize steps to use for 3-D mesh optimization + int optsteps_2d; //!< Number of optimize steps to use for 2-D mesh optimization + + // Philippose - 13/09/2010 + // Added a couple more parameters into the meshing parameters list + // from Netgen into Nglib + int invert_tets; //!< Invert all the volume elements + int invert_trigs; //!< Invert all the surface triangle elements + + int check_overlap; //!< Check for overlapping surfaces during Surface meshing + int check_overlapping_boundary; //!< Check for overlapping surface elements before volume meshing + + + /*! + Default constructor for the Mesh Parameters class + + Note: This constructor initialises the variables in the + class with the following default values + - #uselocalh: 1 + - #maxh: 1000.0 + - #fineness: 0.5 + - #grading: 0.3 + - #elementsperedge: 2.0 + - #elementspercurve: 2.0 + - #closeedgeenable: 0 + - #closeedgefact: 2.0 + - #secondorder: 0 + - #meshsize_filename: null + - #quad_dominated: 0 + - #optsurfmeshenable: 1 + - #optvolmeshenable: 1 + - #optsteps_2d: 3 + - #optsteps_3d: 3 + - #invert_tets: 0 + - #invert_trigs:0 + - #check_overlap: 1 + - #check_overlapping_boundary: 1 + */ + DLL_HEADER Ng_Meshing_Parameters(); + + + + /*! + Reset the meshing parameters to their defaults + + This member function resets all the meshing parameters + of the object to the default values + */ + DLL_HEADER void Reset_Parameters(); + + + + /*! + Transfer local meshing parameters to internal meshing parameters + + This member function transfers all the meshing parameters + defined in the local meshing parameters structure of nglib into + the internal meshing parameters structure used by the Netgen core + */ + DLL_HEADER void Transfer_Parameters(); +}; + + + + +// *** Functions Exported by this Library ************* + +// ------------------------------------------------------------------ +// Netgen library initialisation / destruction functions + +/*! \brief Initialise the Netgen library and prepare for use + + This function needs to be called by the third-party + program before beginning to use the other Netgen + specific functions. +*/ +DLL_HEADER void Ng_Init (); + + +/*! \brief Exit the Netgen meshing kernel in a clean manner + + Use this function to exit the meshing sub-system in + a clean and orderly manner. +*/ +DLL_HEADER void Ng_Exit (); + + +/*! \brief Create a new (and empty) Netgen Mesh Structure + + This function creates a new Netgen Mesh, initialises + it, and returns a pointer to the created mesh structure. + + Use the returned pointer for subsequent operations + which involve mesh operations. + + \return Ng_Mesh Pointer to a Netgen Mesh type #Ng_Mesh +*/ +DLL_HEADER Ng_Mesh * Ng_NewMesh (); + + +/*! \brief Delete an existing Netgen Mesh Structure + + Use this function to delete an existing Netgen mesh + structure and release the used memory. + + \param mesh Pointer to an existing Netgen Mesh structure + of type #Ng_Mesh +*/ +DLL_HEADER void Ng_DeleteMesh (Ng_Mesh * mesh); + + +/*! \brief Save a Netgen Mesh to disk + + This function allows a generated mesh structure to be saved + to disk. + + A Mesh saved using this function, will be written to disk + in the Netgen VOL file format. + + \param mesh Pointer to an existing Netgen Mesh structure + of type #Ng_Mesh + \param filename Pointer to a character array containing the + name of the file to which the mesh should + be saved +*/ +DLL_HEADER void Ng_SaveMesh(Ng_Mesh * mesh, const char* filename); + + +/*! \brief Load a Netgen VOL Mesh from disk into memory + + A Netgen mesh saved in the internal VOL format can be loaded + into a Netgen Mesh structure using this function. + + \param filename Pointer to a character array containing the + name of the file to load + \return Ng_Mesh Pointer to a Netgen Mesh type #Ng_Mesh containing + the mesh loaded from disk +*/ +DLL_HEADER Ng_Mesh * Ng_LoadMesh(const char* filename); + + +/*! \brief Merge a Netgen VOL Mesh from disk into an existing mesh in memory + + A Netgen mesh saved in the internal VOL format can be merged + into an existing Netgen Mesh structure using this function. + + \param mesh Name of the Mesh structure already existent in memory + \param filename Pointer to a character array containing the + name of the file to load + \return Ng_Result Status of the merge operation +*/ +DLL_HEADER Ng_Result Ng_MergeMesh(Ng_Mesh * mesh, const char* filename); + + +/*! \brief Merge one Netgen Mesh into another Netgen Mesh in the case + when both are already in memory + + (NOTE: FUNCTION STILL WORK IN PROGRESS!!!) + + This function can be used to merge two Netgen meshes already present + in memory. + + \param mesh1 Parent Mesh structure into which the second mesh + will be merged + \param mesh2 Child mesh structure which will get merged into + the parent mesh + \return Ng_Result Status of the merge operation +*/ +DLL_HEADER Ng_Result Ng_MergeMesh(Ng_Mesh * mesh1, Ng_Mesh * mesh2); +// ------------------------------------------------------------------ + + + +// ------------------------------------------------------------------ +// Basic Meshing functions for manually adding points, surface elements +// and volume elements to a Netgen Mesh structure + +/*! \brief Add a point to a given Netgen Mesh Structure + + This function allows points to be directly added to a Netgen + mesh structure by providing the co-ordinates. + + Each call to the function allows only one point to be added. + + \param mesh Pointer to an existing Netgen Mesh structure of + type #Ng_Mesh + \param x Pointer to an array of type double containing the co-ordinates + of the point to be added in the form: \n + - x[0] = X co-ordinate + - x[1] = Y co-ordinate + - x[2] = Z co-ordinate +*/ +DLL_HEADER void Ng_AddPoint (Ng_Mesh * mesh, double * x); + + +/*! \brief Add a surface element to a given Netgen Mesh Structure + + This function allows the top-level code to directly add individual + Surface Elements to a Netgen Mesh Structure by providing the type of + element to be added and the indices of the points which constitute the + element. + + Note: + - The points referred to by the surface elements must have been + added prior to calling this function. + - Currently only triangular elements are supported, and the Surface Element + Type argument is not used. + + \param mesh Pointer to an existing Netgen Mesh structure of + type #Ng_Mesh + \param et Surface Element type provided via the enumerated type + #Ng_Surface_Element_Type + \param pi Pointer to an array of integers containing the indices of the + points which constitute the surface element being added +*/ +DLL_HEADER void Ng_AddSurfaceElement (Ng_Mesh * mesh, Ng_Surface_Element_Type et, int * pi); + + +/*! \brief Add a volume element to a given Netgen Mesh Structure + + This function allows the top-level code to directly add individual + Volume Elements to a Netgen Mesh Structure by providing the type of + element to be added and the indices of the points which constitute the + element. + + Note: + - The points referred to by the volume elements must have been + added prior to calling this function. + - Currently only tetrahedral elements are supported, and the Volume Element + Type argument is not used. + + \param mesh Pointer to an existing Netgen Mesh structure of + type #Ng_Mesh + \param et Volume Element type provided via the enumerated type + #Ng_Volume_Element_Type + \param pi Pointer to an array of integers containing the indices of the + points which constitute the volume element being added + +*/ +DLL_HEADER void Ng_AddVolumeElement (Ng_Mesh * mesh, Ng_Volume_Element_Type et, int * pi); + +// ------------------------------------------------------------------ + + + +// ------------------------------------------------------------------ +// Local Mesh Size restriction / limiting utilities + +/*! \brief Apply a global restriction on mesh element size + + This utility allows the user to apply a global mesh element + size limitation. + + During mesh creation, in the absence of an explicit local + size restriction around the neighbourhood of a point within + the meshing domain, this global size restriction will be + utilised. + + Note: This function only limits the Maximum + size of an element within the mesh. + + \param mesh Pointer to an existing Netgen Mesh structure of + type #Ng_Mesh + \param h Variable of type double, specifying the maximum + allowable mesh size +*/ +DLL_HEADER void Ng_RestrictMeshSizeGlobal (Ng_Mesh * mesh, double h); + + +/*! \brief Locally restrict the mesh element size at the given point + + Unlike the function #Ng_RestrictMeshSizeGlobal, this function + allows the user to locally restrict the maximum allowable mesh + size at a given point. + + The point is specified via its three cartesian co-ordinates. + + Note: This function only limits the Maximum size + of the elements around the specified point. + + \param mesh Pointer to an existing Netgen Mesh structure of + type #Ng_Mesh + \param p Pointer to an Array of type double, containing + the three co-ordinates of the point in the form: \n + - p[0] = X co-ordinate + - p[1] = Y co-ordinate + - p[2] = Z co-ordinate + \param h Variable of type double, specifying the maximum + allowable mesh size at that point +*/ +DLL_HEADER void Ng_RestrictMeshSizePoint (Ng_Mesh * mesh, double * p, double h); + + +/*! \brief Locally restrict the mesh element size within a specified box + + Similar to the function #Ng_RestrictMeshSizePoint, this function + allows the size of elements within a mesh to be locally limited. + + However, rather than limit the mesh size at a single point, this + utility restricts the local mesh size within a 3D Box region, specified + via the co-ordinates of the two diagonally opposite points of a cuboid. + + Note: This function only limits the Maximum size + of the elements within the specified region. + + \param mesh Pointer to an existing Netgen Mesh structure of + type #Ng_Mesh + \param pmin Pointer to an Array of type double, containing + the three co-ordinates of the first point of the cuboid: \n + - pmin[0] = X co-ordinate + - pmin[1] = Y co-ordinate + - pmin[2] = Z co-ordinate + \param pmax Pointer to an Array of type double, containing + the three co-ordinates of the opposite point of the + cuboid: \n + - pmax[0] = X co-ordinate + - pmax[1] = Y co-ordinate + - pmax[2] = Z co-ordinate + \param h Variable of type double, specifying the maximum + allowable mesh size at that point +*/ +DLL_HEADER void Ng_RestrictMeshSizeBox (Ng_Mesh * mesh, double * pmin, double * pmax, double h); + +// ------------------------------------------------------------------ + + + +// ------------------------------------------------------------------ +// 3D Mesh Generation functions + +/*! \brief Create a 3D Volume Mesh given a Surface Mesh + + After creating a surface mesh, this function can be utilised + to automatically generate the corresponding 3D Volume Mesh. + + Mesh generation parameters (such as grading, maximum element size, + etc.) are specified via the meshing parameters class which also + needs to be passed to this function. + + Note: Currently, Netgen generates pure tetrahedral volume + meshes. + + \param mesh Pointer to an existing Netgen Mesh structure of + type #Ng_Mesh + \param mp Pointer to a copy of the Meshing Parameters class + (#Ng_Meshing_Parameters), filled up with the + required values + + \return Ng_Result Status of the Mesh Generation routine. More + details regarding the return value can be + found in the description of #Ng_Result +*/ +DLL_HEADER Ng_Result Ng_GenerateVolumeMesh (Ng_Mesh * mesh, Ng_Meshing_Parameters * mp); + +// ------------------------------------------------------------------ + + + +// ------------------------------------------------------------------ +// Basic Mesh information functions + +/*! \brief Returns the Number of Points present in the specified Mesh + + Given an already existent Netgen Mesh Structure, this function + returns the number of points currently present within the Mesh. + + \param mesh Pointer to an existing Netgen Mesh structure of + type #Ng_Mesh + \return + Integer Data-type with the number of points in the Mesh +*/ +DLL_HEADER int Ng_GetNP (Ng_Mesh * mesh); + + +/*! \brief Returns the Number of Surface Elements present in the specified Mesh + + Given an already existent Netgen Mesh Structure, this function + returns the number of surface elements currently present within + the Mesh. + + \param mesh Pointer to an existing Netgen Mesh structure of + type #Ng_Mesh + \return + Integer Data-type with the number of surface elements in the Mesh +*/ +DLL_HEADER int Ng_GetNSE (Ng_Mesh * mesh); + + +/*! \brief Returns the Number of Volume Elements present in the specified Mesh + + Given an already existent Netgen Mesh Structure, this function + returns the number of volume elements currently present within + the Mesh. + + \param mesh Pointer to an existing Netgen Mesh structure of + type #Ng_Mesh + \return + Integer Data-type with the number of volume elements in the Mesh +*/ +DLL_HEADER int Ng_GetNE (Ng_Mesh * mesh); + +// ------------------------------------------------------------------ + + + +// ------------------------------------------------------------------ +// Mesh Topology functions +// Use these functions to extract points, surface / volume elements, +// perform topological searches, etc..etc... + +// Return the Point Coordinates of a specified Point +// The x, y and z co-ordinates are returned in the array pointer as +// x[0] = x ; x[1] = y ; x[2] = z +DLL_HEADER void Ng_GetPoint (Ng_Mesh * mesh, int num, double * x); + + + +// return surface and volume element in pi +DLL_HEADER Ng_Surface_Element_Type +Ng_GetSurfaceElement (Ng_Mesh * mesh, int num, int * pi); + +DLL_HEADER Ng_Volume_Element_Type +Ng_GetVolumeElement (Ng_Mesh * mesh, int num, int * pi); + +// ------------------------------------------------------------------ + + + + +// ********************************************************** +// ** 2D Meshing ** +// ********************************************************** + + +// feeds points and boundary to mesh + +DLL_HEADER void Ng_AddPoint_2D (Ng_Mesh * mesh, double * x); +DLL_HEADER void Ng_AddBoundarySeg_2D (Ng_Mesh * mesh, int pi1, int pi2); + +// ask for number of points, elements and boundary segments +DLL_HEADER int Ng_GetNP_2D (Ng_Mesh * mesh); +DLL_HEADER int Ng_GetNE_2D (Ng_Mesh * mesh); +DLL_HEADER int Ng_GetNSeg_2D (Ng_Mesh * mesh); + +// return point coordinates +DLL_HEADER void Ng_GetPoint_2D (Ng_Mesh * mesh, int num, double * x); + +// return 2d elements +DLL_HEADER Ng_Surface_Element_Type +Ng_GetElement_2D (Ng_Mesh * mesh, int num, int * pi, int * matnum = NULL); + +// return 2d boundary segment +DLL_HEADER void Ng_GetSegment_2D (Ng_Mesh * mesh, int num, int * pi, int * matnum = NULL); + + +// load 2d netgen spline geometry +DLL_HEADER Ng_Geometry_2D * Ng_LoadGeometry_2D (const char * filename); + +// generate 2d mesh, mesh is allocated by function +DLL_HEADER Ng_Result Ng_GenerateMesh_2D (Ng_Geometry_2D * geom, + Ng_Mesh ** mesh, + Ng_Meshing_Parameters * mp); + +DLL_HEADER void Ng_HP_Refinement (Ng_Geometry_2D * geom, + Ng_Mesh * mesh, + int levels); + + + + + +// ********************************************************** +// ** STL Meshing ** +// ********************************************************** + + +// loads geometry from STL file +DLL_HEADER Ng_STL_Geometry * Ng_STL_LoadGeometry (const char * filename, int binary = 0); + + +// generate new STL Geometry +DLL_HEADER Ng_STL_Geometry * Ng_STL_NewGeometry (); + + +// fills STL Geometry +// positive orientation +// normal vector may be null-pointer +DLL_HEADER void Ng_STL_AddTriangle (Ng_STL_Geometry * geom, + double * p1, double * p2, double * p3, + double * nv = NULL); + +// add (optional) edges : +DLL_HEADER void Ng_STL_AddEdge (Ng_STL_Geometry * geom, + double * p1, double * p2); + +// after adding triangles (and edges) initialize +DLL_HEADER Ng_Result Ng_STL_InitSTLGeometry (Ng_STL_Geometry * geom); + +// automatically generates edges: +DLL_HEADER Ng_Result Ng_STL_MakeEdges (Ng_STL_Geometry * geom, + Ng_Mesh* mesh, + Ng_Meshing_Parameters * mp); + + +// generates mesh, empty mesh must be already created. +DLL_HEADER Ng_Result Ng_STL_GenerateSurfaceMesh (Ng_STL_Geometry * geom, + Ng_Mesh * mesh, + Ng_Meshing_Parameters * mp); + + +#ifdef ACIS + +// ********************************************************** +// ** ACIS Meshing ** +// ********************************************************** + +/// Data type for NETGEN STL geomty +typedef void * Ng_ACIS_Geometry; + +// loads geometry from STL file +DLL_HEADER Ng_ACIS_Geometry * Ng_ACIS_LoadGeometry (const char * filename); + +// generates mesh, empty mesh must be already created. +DLL_HEADER Ng_Result Ng_ACIS_GenerateSurfaceMesh (Ng_ACIS_Geometry * geom, + Ng_Mesh * mesh, + Ng_Meshing_Parameters * mp); + + +#endif + + + +#ifdef OCCGEOMETRY + +// ********************************************************** +// ** OpenCascade Geometry / Meshing Utilities ** +// ********************************************************** + +// Create new OCC Geometry Object +DLL_HEADER Ng_OCC_Geometry * Ng_OCC_NewGeometry (); + +// Delete an OCC Geometry Object +DLL_HEADER Ng_Result Ng_OCC_DeleteGeometry (Ng_OCC_Geometry * geom); + +// Loads geometry from STEP file +DLL_HEADER Ng_OCC_Geometry * Ng_OCC_Load_STEP (const char * filename); + +// Loads geometry from IGES file +DLL_HEADER Ng_OCC_Geometry * Ng_OCC_Load_IGES (const char * filename); + +// Loads geometry from BREP file +DLL_HEADER Ng_OCC_Geometry * Ng_OCC_Load_BREP (const char * filename); + +// Set the local mesh size based on geometry / topology +DLL_HEADER Ng_Result Ng_OCC_SetLocalMeshSize (Ng_OCC_Geometry * geom, + Ng_Mesh * mesh, + Ng_Meshing_Parameters * mp); + +// Mesh the edges and add Face descriptors to prepare for surface meshing +DLL_HEADER Ng_Result Ng_OCC_GenerateEdgeMesh (Ng_OCC_Geometry * geom, + Ng_Mesh * mesh, + Ng_Meshing_Parameters * mp); + +// Mesh the surfaces of an OCC geometry +DLL_HEADER Ng_Result Ng_OCC_GenerateSurfaceMesh (Ng_OCC_Geometry * geom, + Ng_Mesh * mesh, + Ng_Meshing_Parameters * mp); + +// Get the face map of an already loaded OCC geometry +DLL_HEADER Ng_Result Ng_OCC_GetFMap(Ng_OCC_Geometry * geom, + Ng_OCC_TopTools_IndexedMapOfShape * FMap); + +#endif // OCCGEOMETRY + + + +// ********************************************************** +// ** Mesh refinement algorithms ** +// ********************************************************** + +// uniform mesh refinement +DLL_HEADER void Ng_Uniform_Refinement (Ng_Mesh * mesh); + + +// uniform mesh refinement with geometry adaption: + +DLL_HEADER void Ng_2D_Uniform_Refinement (Ng_Geometry_2D * geom, + Ng_Mesh * mesh); + +DLL_HEADER void Ng_STL_Uniform_Refinement (Ng_STL_Geometry * geom, + Ng_Mesh * mesh); + +DLL_HEADER void Ng_CSG_Uniform_Refinement (Ng_CSG_Geometry * geom, + Ng_Mesh * mesh); + +#ifdef OCCGEOMETRY +DLL_HEADER void Ng_OCC_Uniform_Refinement (Ng_OCC_Geometry * geom, + Ng_Mesh * mesh); +#endif + + + +// ********************************************************** +// ** Second Order mesh algorithms ** +// ********************************************************** + +// convert mesh to second order +DLL_HEADER void Ng_Generate_SecondOrder (Ng_Mesh * mesh); + + +// convert mesh to second order with geometry adaption: + +DLL_HEADER void Ng_2D_Generate_SecondOrder (Ng_Geometry_2D * geom, + Ng_Mesh * mesh); + +DLL_HEADER void Ng_STL_Generate_SecondOrder (Ng_STL_Geometry * geom, + Ng_Mesh * mesh); + +DLL_HEADER void Ng_CSG_Generate_SecondOrder (Ng_CSG_Geometry * geom, + Ng_Mesh * mesh); + +#ifdef OCCGEOMETRY +DLL_HEADER void Ng_OCC_Generate_SecondOrder (Ng_OCC_Geometry * geom, + Ng_Mesh * mesh); +#endif + + +#endif // NGLIB diff --git a/ML/Training/include/pch.h b/ML/Training/include/pch.h new file mode 100644 index 0000000..1fa98bc --- /dev/null +++ b/ML/Training/include/pch.h @@ -0,0 +1,570 @@ +#pragma once +#ifndef PCH_H +#define PCH_H +#include "framework.h" +#endif //PCH_H + +#define HX_API extern "C" _declspec(dllexport) +#include +#include +#include +#include +#include +#include +#include +#include +#include + +//const double M_PI = acos(-1.0); +typedef std::vector>>dVec3; //三维数组:double +typedef std::vector>dVec2; //二维数组:double +typedef std::vector>iVec2; //二维数组:int +typedef std::vectordVec1; //一维数组:double +typedef std::vectoriVec1; //一维数组:int + +template void HX_copy(std::vector>>& p1, const std::vector>>& p0) +{ + int m = p0.size(); p1.resize(m); + for (int i = 0; i < m; ++i) + { + int n = p0[i].size(); p1[i].resize(n); + for (int j = 0; j < n; ++j) + { + int l = p0[i][j].size(); p1[i][j].resize(l); + for (int k = 0; k < l; ++k) + { + p1[i][j][k] = p0[i][j][k]; + } + } + } +} +template void HX_copy(std::vector>& p1, const std::vector>& p0) +{ + int m = p0.size(); p1.resize(m); + for (int i = 0; i < m; ++i) + { + int n = p0[i].size(); p1[i].resize(n); + for (int j = 0; j < n; ++j) + { + p1[i][j] = p0[i][j]; + } + } +} +template void HX_copy(std::vector& p1, const std::vector& p0) +{ + int m = p0.size(); p1.resize(m); + for (int i = 0; i < m; ++i) + { + p1[i] = p0[i]; + } +} +template void HX_copy(std::vector>>& p1, T*** p0, int m, int* n, int l) +{ + p1.resize(m); + for (int i = 0; i < m; ++i) + { + p1[i].resize(n[i]); + for (int j = 0; j < n[i]; ++j) + { + p1[i][j].resize(l); + for (int k = 0; k < l; ++k) + { + p1[i][j][k] = p0[i][j][k]; + } + } + } +} +template void HX_copy(std::vector>& p1, T** p0, int m, int n) +{ + p1.resize(m); + for (int i = 0; i < m; ++i) + { + p1[i].resize(n); + for (int j = 0; j < n; ++j) + { + p1[i][j] = p0[i][j]; + } + } +} +template void HX_copy(std::vector& p1, T* p0, int m) +{ + p1.resize(m); + for (int i = 0; i < m; ++i) + { + p1[i] = p0[i]; + } +} + +//点结构体 +struct point +{ + //点结构体 + double x; double y; //点坐标 + point() { x = 0; y = 0; } + ~point() {} + void set(const double& x_ = 0, const double& y_ = 0) { x = x_; y = y_; } + void set(const point& p) { x = p.x; y = p.y; } +}; +struct point3 +{ + double x; + double y; + double z; + point3() { x = 0; y = 0; z = 0; } + ~point3() {} + point3(const dVec1&p) { x = p[0]; y = p[1]; z = p[2]; } + point3(const double& x_ = 0, const double& y_ = 0, const double&z_ = 0) { x = x_; y = y_; z = z_;} + void set(const double& x_ = 0, const double& y_ = 0, const double& z_ = 0) { x = x_; y = y_; z = z_; } + void set(const point3&p) { x = p.x; y = p.y; z = p.z; } +}; +//网格结构体 +struct cell +{ + //网格单元结构体 + std::vector p; + iVec2 pindex; + iVec1 isplot; + cell() {} + ~cell() {} +}; + + +//网格算法输入参数结构体 +struct HX_NWTM_GRID_INPUT +{ + // 网格划分算法输入参数结构体 + dVec2 Boundary; //3D:{x0, y0, z0, x1, y1, z1} //2D:{x0, y0, x1, y1} 边界数据 + dVec2 VerticalWell; //3D:{x0, y0, z0, x1, y1, z1, rw} //2D:{x0, y0, x1, y1, rw} 直井数据 + dVec2 HorizontalWell; //3D:{x0, y0, z0, x1, y1, z1, rw} 水平井数据 + dVec2 FractureVerticalWell; //3D:{x0, y0, z0, x1, y1, z1, wf} //2D:{x0, y0, x1, y1, wf} 压裂直井数据 + dVec3 MultistageFracturedHorizontalWell; //3D:{x0, y0, z0, x1, y1, z1, wf} //2D:{x0, y0, x1, y1, wf} 多级压裂水平井数据 + dVec2 InclinedWell; //3D:{x0, y0, z0, x1, y1, z1, rw} 斜井数据 + dVec2 Fault; //3D:{x0, y0, z0, x1, y1, z1} //2D:{x0, y0, x1, y1} 断层数据 + double GridControl; // 网格大小控制参数 + int D; // 维数 + + //默认初始化 + HX_NWTM_GRID_INPUT() + { + dVec1 a(3), b(4), c(5); + Boundary.resize(4); + b[0] = -1500.0; b[1] = -1500.0; b[2] = -1500.0; b[3] = 1500.0; Boundary[0] = b; + b[0] = -1500.0; b[1] = 1500.0; b[2] = 1500.0; b[3] = 1500.0; Boundary[1] = b; + b[0] = 1500.0; b[1] = 1500.0; b[2] = 1500.0; b[3] = -1500.0; Boundary[2] = b; + b[0] = 1500.0; b[1] = -1500.0; b[2] = -1500.0; b[3] = -1500.0; Boundary[3] = b; + + /*VerticalWell.resize(3); + a[0] = 0; a[1] = 0; a[2] = 0.1; VerticalWell[0] = a; + a[0] = 1000; a[1] = 1000; a[2] = 0.1; VerticalWell[1] = a; + a[0] = -1000; a[1] = -1000; a[2] = 0.1; VerticalWell[2] = a; + + HorizontalWell.resize(0); + + FractureVerticalWell.resize(1); + c[0] = -200; c[1] = -200; c[2] = 200; c[3] = -200; c[4] = 0.05; FractureVerticalWell[0] = c; + + MultistageFracturedHorizontalWell.resize(1); + MultistageFracturedHorizontalWell[0].resize(3, dVec1(5)); + c[0] = -600; c[1] = 600; c[2] = -400; c[3] = 600; c[4] = 0.1; MultistageFracturedHorizontalWell[0][0] = c; + c[0] = -600; c[1] = 400; c[2] = -400; c[3] = 400; c[4] = 0.1; MultistageFracturedHorizontalWell[0][1] = c; + c[0] = -600; c[1] = 200; c[2] = -400; c[3] = 200; c[4] = 0.1; MultistageFracturedHorizontalWell[0][2] = c; + + InclinedWell.resize(0); + + Fault.resize(1); + c[0] = -500; c[1] = 1000; c[2] = 500; c[3] = 500; Fault[0] = c;*/ + //单一直井 + VerticalWell.resize(1); + a[0] = 0; a[1] = 0; a[2] = 0.1; VerticalWell[0] = a; + + HorizontalWell.resize(0); + + FractureVerticalWell.resize(0); + + MultistageFracturedHorizontalWell.resize(0); + + InclinedWell.resize(0); + + Fault.resize(0); + //单一压裂直井 + /*VerticalWell.resize(0); + + HorizontalWell.resize(0); + + FractureVerticalWell.resize(1); + c[0] = -200; c[1] = 0; c[2] = 200; c[3] = 0; c[4] = 0.05; FractureVerticalWell[0] = c; + + MultistageFracturedHorizontalWell.resize(0); + + InclinedWell.resize(0); + + Fault.resize(0);*/ + //单一多段压裂水平井 + /*VerticalWell.resize(0); + + HorizontalWell.resize(0); + + FractureVerticalWell.resize(0); + + MultistageFracturedHorizontalWell.resize(1); + MultistageFracturedHorizontalWell[0].resize(7, dVec1(5)); + c[0] = -600; c[1] = -200; c[2] = -600; c[3] = 200; c[4] = 0.05; MultistageFracturedHorizontalWell[0][0] = c; + c[0] = -400; c[1] = -200; c[2] = -400; c[3] = 200; c[4] = 0.05; MultistageFracturedHorizontalWell[0][1] = c; + c[0] = -200; c[1] = -200; c[2] = -200; c[3] = 200; c[4] = 0.05; MultistageFracturedHorizontalWell[0][2] = c; + c[0] = 0; c[1] = -200; c[2] =0; c[3] = 200; c[4] = 0.05; MultistageFracturedHorizontalWell[0][3] = c; + c[0] = 200; c[1] = -200; c[2] = 200; c[3] = 200; c[4] = 0.05; MultistageFracturedHorizontalWell[0][4] = c; + c[0] = 400; c[1] = -200; c[2] = 400; c[3] = 200; c[4] = 0.05; MultistageFracturedHorizontalWell[0][5] = c; + c[0] = 600; c[1] = -200; c[2] = 600; c[3] = 200; c[4] = 0.05; MultistageFracturedHorizontalWell[0][6] = c; + InclinedWell.resize(0); + + Fault.resize(0);*/ + + //单一直井+断层 + /*VerticalWell.resize(1); + a[0] = 0; a[1] = 0; a[2] = 0.1; VerticalWell[0] = a; + + HorizontalWell.resize(0); + + FractureVerticalWell.resize(0); + + MultistageFracturedHorizontalWell.resize(0); + + InclinedWell.resize(0); + + Fault.resize(1); + c[0] = -50; c[1] = -100; c[2] = -50; c[3] = 100; Fault[0] = c;*/ + //单一压裂直井+断层 + /*VerticalWell.resize(0); + + HorizontalWell.resize(0); + + FractureVerticalWell.resize(1); + c[0] = -200; c[1] = 0; c[2] = 200; c[3] = 0; c[4] = 0.05; FractureVerticalWell[0] = c; + + MultistageFracturedHorizontalWell.resize(0); + + InclinedWell.resize(0); + + Fault.resize(1); + c[0] = -100; c[1] = 100; c[2] = 100; c[3] = 100; Fault[0] = c;*/ + //单一多段压裂水平井+断层 + /*VerticalWell.resize(0); + + HorizontalWell.resize(0); + + FractureVerticalWell.resize(0); + + MultistageFracturedHorizontalWell.resize(1); + MultistageFracturedHorizontalWell[0].resize(7, dVec1(5)); + c[0] = -600; c[1] = -200; c[2] = -600; c[3] = 200; c[4] = 0.05; MultistageFracturedHorizontalWell[0][0] = c; + c[0] = -400; c[1] = -200; c[2] = -400; c[3] = 200; c[4] = 0.05; MultistageFracturedHorizontalWell[0][1] = c; + c[0] = -200; c[1] = -200; c[2] = -200; c[3] = 200; c[4] = 0.05; MultistageFracturedHorizontalWell[0][2] = c; + c[0] = 0; c[1] = -200; c[2] =0; c[3] = 200; c[4] = 0.05; MultistageFracturedHorizontalWell[0][3] = c; + c[0] = 200; c[1] = -200; c[2] = 200; c[3] = 200; c[4] = 0.05; MultistageFracturedHorizontalWell[0][4] = c; + c[0] = 400; c[1] = -200; c[2] = 400; c[3] = 200; c[4] = 0.05; MultistageFracturedHorizontalWell[0][5] = c; + c[0] = 600; c[1] = -200; c[2] = 600; c[3] = 200; c[4] = 0.05; MultistageFracturedHorizontalWell[0][6] = c; + InclinedWell.resize(0); + + Fault.resize(1); + c[0] = -300; c[1] = 300; c[2] = 300; c[3] = 300; Fault[0] = c;*/ + + + GridControl = 150.0; + + D = 2; + } + ~HX_NWTM_GRID_INPUT() {} +}; + +//网格算法输出参数结构体(绘图用) +struct HX_NWTM_GRID_OUTPUT1 +{ + cell TRI_cell; //三角形网格 + cell PEBI_cell; //PEBI网格 + + HX_NWTM_GRID_OUTPUT1() {} + ~HX_NWTM_GRID_OUTPUT1() {} +}; + +//网格算法输出参数结构体(模型用) +struct HX_NWTM_GRID_OUTPUT2 +{ + dVec2 Trinodexy; + dVec1 Area; + dVec2 D; + struct { + int n; + dVec2 XiLinw; + dVec2 lw; + dVec2 dw; + dVec1 rw; + iVec2 inwell; + } ZhiJingNeiBianJie; + struct { + int n; + dVec2 XiLinf; + dVec2 lf; + dVec2 df; + dVec1 xf; + iVec2 infra; + } LieFengJingNeiBianJie; + struct { + int n; + dVec2 XiLinh; + dVec2 lh; + dVec2 dh; + dVec2 dsxf; + iVec2 inhor; + iVec1 nhor; + } DuoJiYaLieShuiPingJingNeiBianJie; + struct { + int n; + dVec2 WaiBianh; + dVec2 WaiBianl; + dVec2 WaiBiand; + } WaiBianJie; + struct { + int n; + dVec2 faultb1; + dVec2 faultb2; + dVec2 faultl1; + dVec2 faultd1; + } NeiBuDuanCeng; + struct { + iVec1 ia; + iVec1 ja; + iVec2 nzeros; + int numk; + } YuChuLiJuZhen; + HX_NWTM_GRID_OUTPUT2() {} + ~HX_NWTM_GRID_OUTPUT2() {} +}; + + +//KRINGING插值输入参数结构体 +struct HX_KRING_INPUT +{ + double nugget; //块金值:表示空间点在零距离处的变异程度，即测量误差和小于采样尺度的随机变异之和 + double sill; //基台值:表示变差函数随距离增加而趋于稳定的极限值，反映区域化变量的总变异程度 + double range; //变程:表示空间相关性的有效距离。当两点间距离超过 range 时，它们之间不再具有空间相关性 + double model; //变差函数模型:指定变差函数的数学形式，描述空间相关性随距离的变化规律{, , } + //高斯模型SPHERICAL(0)：相关性随距离增加呈指数衰减，适用于连续性较强的变量 + //指数模型EXPONENTIAL(1)：相关性快速衰减，适用于局部变异性较大的变量 + //球状模型GAUSSIAN(2)：在变程内呈抛物线变化，超过变程后相关性为零 + dVec2 p; //插值点 + dVec2 v; // + ~HX_KRING_INPUT() {} + HX_KRING_INPUT(double nugget0, double sill0, double range0, double model0 , const dVec2& p0, const dVec2& v0) + { + nugget = nugget0; sill = sill0; range = range0; model = model0; + p = p0; + v = v0; + } +}; + +//KRINGING插值输出参数结构体 +struct HX_KRING_OUTPUT +{ + dVec1 v; + HX_KRING_OUTPUT() {} + ~HX_KRING_OUTPUT() {} +}; + +//数值试井模型求解器输入参数结构体 +struct HX_NWTM_MODEL_INPUT +{ + int T; //1:油单相常数pvt; 2:油单相变化pvt; 3:水单相常数pvt; 4:水单相变化pvt; 5:气单相变化pvt; 6:气单相拟压力; 7:油气两相; 8:油水两相; 9:气水两相; 10:油气水三相 + HX_NWTM_GRID_OUTPUT2 GRID; + struct Rate //流量数据 + { + dVec2 t; //时间, h [一口井一组数] + dVec2 qo; //油流量,m^3/d [一口井一组数] + dVec2 qg; //气流量,m^3/d [一口井一组数] + dVec2 qw; //水流量,m^3/d [一口井一组数] + }Rate; + struct Pressure //压力数据 + { + dVec2 t; //时间, h [一口井一组数] + dVec2 p; //压力, MPa [一口井一组数] + }Pressure; + struct CS //井储表皮数据 + { + dVec1 C; //井储, m^3/MPa [一口井一个数] + dVec1 S; //表皮, [一口井一个数] + }CS; + struct PVT //流体性质数据 + { + dVec1 p; //压力, MPa + double pb; //饱和压力, MPa + dVec1 Rso; //溶解气油比, m^3/m^3 + dVec1 Bo; //油体积系数, m^3/m^3 + dVec1 Co; //油压缩系数, 1/MPa + dVec1 miuo; //油粘度, mPa·s + dVec1 rouo; //油密度, kg/m^3 + dVec1 Rv; //凝析油气比, m^3/m^3 + dVec1 Bg; //气体积系数, m^3/m^3 + dVec1 Cg; //气压缩系数, 1/MPa + dVec1 miug; //气粘度, mPa·s + dVec1 roug; //气密度, kg/m^3 + dVec1 Z; //气偏差因子, 1 + dVec1 Rsw; //溶解气水比, m^3/m^3 + dVec1 Bw; //水体积系数, m^3/m^3 + dVec1 Cw; //水压缩系数, 1/MPa + dVec1 miuw; //水粘度, mPa·s + dVec1 rouw; //水密度, kg/m^3 + dVec1 V; //吸附气量, m^3/kg + dVec1 k_kinitial; //渗透率比, 1 + dVec1 Cf_Cfinitial; //岩石压缩系数比, 1 + dVec1 So; //油饱和度 + dVec1 Kro; //油相对渗透率 + dVec1 Sg; //气饱和度 + dVec1 Krg; //气相对渗透率 + dVec1 Sw; //水饱和度 + dVec1 Krw; //水相对渗透率 + }PVT; + struct Base //基础数据 + { + double Pi; //初始压力, MPa + double Cti; //综合压缩系数, 1/MPa + double Cf; //岩石压缩系数, 1/MPa + double Soi; //初始含油饱和度 + double Sgi; //初始含气饱和度 + double Swi; //初始含水饱和度 + dVec1 k; //渗透率, D [一个网格单元一个值] + dVec1 phi; //孔隙度, 1 [一个网格单元一个值] + dVec1 h; //储层厚度, m [一个网格单元一个值] + double d; //时间增长指数 + double dt_Min; //最小时间间隔, h + double dt_Max; //最大时间间隔, h + }Base; + + //初始化 + HX_NWTM_MODEL_INPUT() {} + ~HX_NWTM_MODEL_INPUT() {} + HX_NWTM_MODEL_INPUT(const HX_NWTM_GRID_OUTPUT2& p0) + { + T =1; + GRID = p0; + + /*Rate.t.resize(5); + Rate.qo.resize(5); + Rate.qw.resize(5); + Rate.qg.resize(5); + + Rate.t[0].resize(2); Rate.t[0][0] = 2000; Rate.t[0][1] = 500; + Rate.qo[0].resize(2); Rate.qo[0][0] = 10; Rate.qo[0][1] = 0; + Rate.qw[0].resize(2); Rate.qw[0][0] = 2; Rate.qw[0][1] = 0; + Rate.qg[0].resize(2); Rate.qg[0][0] = 20000; Rate.qg[0][1] = 0; + + Rate.t[1].resize(0); + Rate.qo[1].resize(0); + Rate.qw[1].resize(0); + Rate.qg[1].resize(0); + + Rate.t[2].resize(3); Rate.t[2][0] = 1000; Rate.t[2][1] = 1000; Rate.t[2][2] = 500; + Rate.qo[2].resize(3); Rate.qo[2][0] = 30; Rate.qo[2][1] = 40; Rate.qo[2][2] = 20; + Rate.qw[2].resize(3); Rate.qw[2][0] = 3; Rate.qw[2][1] = 4; Rate.qw[2][2] = 2; + Rate.qg[2].resize(3); Rate.qg[2][0] = 30000; Rate.qg[2][1] = 40000; Rate.qg[2][2] = 20000; + + Rate.t[3].resize(2); Rate.t[3][0] = 1500; Rate.t[3][1] = 1000; + Rate.qo[3].resize(2); Rate.qo[3][0] = 30; Rate.qo[3][1] = 20; + Rate.qw[3].resize(2); Rate.qw[3][0] = 5; Rate.qw[3][1] = 2; + Rate.qg[3].resize(2); Rate.qg[3][0] = 50000; Rate.qg[3][1] = 20000; + Rate.t[4].resize(2); Rate.t[4][0] = 1000; Rate.t[4][1] = 1500; + Rate.qo[4].resize(2); Rate.qo[4][0] = -50; Rate.qo[4][1] = -60; + Rate.qw[4].resize(2); Rate.qw[4][0] = -2; Rate.qw[4][1] = -5; + Rate.qg[4].resize(2); Rate.qg[4][0] = -20000; Rate.qg[4][1] = -50000; + + Pressure.t.resize(0); + Pressure.p.resize(0); + + CS.C.resize(5); + CS.C[0] = 0.1; CS.C[1] = 0.1; CS.C[2] = 0.1; CS.C[3] = 0.1; CS.C[4] = 0.1; + CS.S.resize(5); + CS.S[0] = 0.1; CS.S[1] = 0.1; CS.S[2] = 0.1; CS.S[3] = 0.1; CS.S[4] = 0.1;*/ + + + Rate.t.resize(1); + Rate.qo.resize(1); + Rate.qw.resize(1); + Rate.qg.resize(1); + + Rate.t[0].resize(2); Rate.t[0][0] = 2000; Rate.t[0][1] = 500; + Rate.qo[0].resize(2); Rate.qo[0][0] =10; Rate.qo[0][1] = 0; + Rate.qw[0].resize(2); Rate.qw[0][0] =-10; Rate.qw[0][1] = 0; + Rate.qg[0].resize(2); Rate.qg[0][0] = 50000; Rate.qg[0][1] = 0; + Pressure.t.resize(0); + Pressure.p.resize(0); + + CS.C.resize(1); + CS.C[0] = 0.1; + CS.S.resize(1); + CS.S[0] = 0.1; + + PVT.p = dVec1(200, 0); for (int i = 0; i < 200; ++i) { PVT.p[i] = (i + 1.0); } + PVT.pb = 40.0; + PVT.Rso = dVec1(200, 0); + PVT.Bo = dVec1(200, 1.2); + PVT.Co = dVec1(200, 5e-4); + PVT.miuo = dVec1(200, 0.5); + PVT.rouo = dVec1(200, 800); + PVT.Rv = dVec1(200, 0); + PVT.Bg = dVec1(200, 5e-3); + PVT.Cg = dVec1(200, 2e-2); + PVT.miug = dVec1(200, 2e-2); + PVT.roug = dVec1(200, 200); + PVT.Z = dVec1(200, 1); + PVT.Rsw = dVec1(200, 0); + PVT.Bw = dVec1(200, 1.05); + PVT.Cw = dVec1(200, 1e-4); + PVT.miuw = dVec1(200, 0.8); + PVT.rouw = dVec1(200, 1000); + PVT.V = dVec1(200, 0); + + PVT.k_kinitial = dVec1(200, 1); + PVT.Cf_Cfinitial = dVec1(200, 1); + + PVT.So = dVec1(81, 0); for (int i = 0; i < 81; ++i) { PVT.So[i] = (0.1+i*0.01); } + PVT.Kro = dVec1(81, 0); for (int i = 0; i < 81; ++i) { PVT.Kro[i] = (i*0.0125); } + PVT.Krw = dVec1(81, 0); for (int i = 0; i < 81; ++i) { PVT.Krw[i] = (1 - i * 0.0125); } + PVT.Sg = dVec1(100, 0); + PVT.Krg = dVec1(100, 0); + PVT.Sw = dVec1(100, 0); + + + + Base.Pi = 40.0; + Base.Cti = 1e-3; + Base.Cf = 1e-4; + Base.Soi = 0.8; + Base.Sgi = 0.0; + Base.Swi = 0.2; + Base.k = dVec1(p0.Trinodexy.size(), 0.001); + Base.phi = dVec1(p0.Trinodexy.size(), 0.1); + Base.h = dVec1(p0.Trinodexy.size(), 10); + Base.d = 1.05; + Base.dt_Min = 0.0025; + Base.dt_Max = 12.5; + } +}; + +//数值试井模型求解器输出参数结构体 +struct HX_NWTM_MODEL_OUTPUT +{ + dVec1 t; //时间, h + dVec2 pw; //井底压力, MPa [一口井一组数] + dVec2 p; //压力分布, MPa [一个时间一组数] + dVec2 So; //油饱和度分布 [一个时间一组数] + dVec2 Sg; //气饱和度分布 [一个时间一组数] + dVec2 Sw; //水饱和度分布 [一个时间一组数] + dVec2 k; //渗透率分布,mD [一个时间一组数] + HX_NWTM_MODEL_OUTPUT() {} + ~HX_NWTM_MODEL_OUTPUT() {} +}; + +HX_API void HX_NWTM_GRID(HX_NWTM_GRID_OUTPUT1& p1, HX_NWTM_GRID_OUTPUT2& p2, const HX_NWTM_GRID_INPUT& p0, std::string LIC); //数值试井网格接口 +HX_API void HX_NWTM_KRINGING(HX_KRING_OUTPUT& p1, const HX_KRING_INPUT p0, std::string LIC); //数值试井非均质性计算接口 +HX_API void HX_NWTM_MODEL(HX_NWTM_MODEL_OUTPUT& p1, const HX_NWTM_MODEL_INPUT& p0, std::string LIC); //数值试井模型求解器接口 + + + diff --git a/ML/Training/include/singlePhaseSolver.h b/ML/Training/include/singlePhaseSolver.h new file mode 100644 index 0000000..ca9be3b --- /dev/null +++ b/ML/Training/include/singlePhaseSolver.h @@ -0,0 +1,396 @@ +#pragma once +#include +#include +#include +#include +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 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 pointIndices; + std::vector 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 a_sorted(pointIndices); + std::vector 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 vector_interPoint; // 交点坐标 + + //裂缝上的线单元 + std::vector fracCells; + + double dE;//裂缝宽度 + double dKr;//裂缝渗透率比,Kf/K + double dW;//裂缝储能比,(phi*Ct)f/(phi*Ct) + double dFc;//裂缝导流能力,dFc=dKr*dE; + + //记录netgen的in2d文件中几何点的编号 + //一个裂缝只需要2个节点即可记录 + std::vector 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 vecNodes;//边界线上的网格节点编号 +}; + +// 外边界的基类 +class CBoundShape { + public: + BOUND_SHAPE boundShape; // 定义外边界类型，0为圆形，1为四边形，2为多边形 + + // 表示多边形（包括四边形）的参数 + int nNumSegs; // 边的数量 + std::vector vecSegments; // 边 + + // 表示圆形 + Point cCenter; + double dRadius; + int nNodeNum; + BOUND_TYPE boundType; + double dConstPressure;//圆形定压边界时的压力值 + std::vector vecNodes;//圆形边界上的网格点编号 + + //记录netgen的in2d文件中几何点的编号 + std::vector ngPointIndex; +}; + +// 多边形,内部的复合区域 +class CBoundPolygon : public CBoundShape { + public: + CBoundPolygon() { + bAnchor = false; + } + + CBoundPolygon(std::vector 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 vecSegments; // 边 + + double dComKr;//复合区的渗透率比 + double dComW;//复合区的储能比 + bool bAnchor;//是否作为复合区进行材料属性的分区赋值，false表示不进行复合区的材料编号赋值，即其材料编号与背景值相同 + + //记录netgen的in2d文件中几何点的编号 + std::vector 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 m_NodeIndex;//all the nodes index of the wellbore + std::vector 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 vecFracs;//压裂裂缝 + double dWidth;//2d模型中水平井宽度 + CBoundPolygon wellPoly;//水平井的边界形状 + + //记录netgen的in2d文件中几何点的编号 + std::vector 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& wells, std::vector limits, + std::vector& faults, std::vector& fracs, CBoundShape &shape); +__declspec(dllexport) void Ng2VTK(Ng_Mesh *mesh); + + +//singlePhase solver API +__declspec(dllexport) int singlePhaseSolver(double* pBaseData, std::vector points, std::vector cells, + int m_nWellNum, CBoundWell* pwells, std::vector>& vecBP, + std::vector>>& vecNodePre); +__declspec(dllexport) int singlePhaseSolverNgmesh(const TimeDis time, double* m_pBaseData, nglib::Ng_Mesh* mesh, + std::vector& m_wWell, const std::vector& comPolygon, const std::vector& faultPolygon, + std::vector& fracs, CBoundShape &outBound, std::vector>>& 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& m_wWell, + const std::vector& comPolygon, const std::vector& faultPolygon, std::vector& fracs, + CBoundShape &outBound, std::vector>>& vecNodePre, int& timeStep, int& totalSteps); +__declspec(dllexport) bool logLogPre(const std::vector& pressure, const int& nSection, double* pdTimeQ, + double* pdQ, int nTimeQ, std::vector& logPre); +#ifdef __cplusplus +} + +#endif diff --git a/ML/Training/include/soldata.hpp b/ML/Training/include/soldata.hpp new file mode 100644 index 0000000..1a7dfd3 --- /dev/null +++ b/ML/Training/include/soldata.hpp @@ -0,0 +1,130 @@ +#ifndef FILE_SOLDATA +#define FILE_SOLDATA + + +namespace netgen +{ + + using namespace std; + + class SolutionData + { + protected: + + string name; + int components; + bool iscomplex; + + int multidimcomponent; + + public: + SolutionData (const string & aname, + int acomponents = 1, bool aiscomplex = 0) + : name(aname), components(acomponents), iscomplex(aiscomplex) + { ; } + + virtual ~SolutionData () + { ; } + + int GetComponents() + { + return components; + } + + bool IsComplex() + { + return iscomplex; + } + + virtual bool GetValue (int /* elnr */, + double /* lam1 */, double /* lam2 */, double /* lam3 */, + double * /* values */) + { + return false; + } + + virtual bool GetValue (int selnr, + const double xref[], const double x[], const double dxdxref[], + double * values) + { + return GetValue (selnr, xref[0], xref[1], xref[2], values); + } + + virtual bool GetMultiValue (int elnr, int facetnr, int npts, + const double * xref, int sxref, + const double * x, int sx, + const double * dxdxref, int sdxdxref, + double * values, int svalues) + { + bool res = false; + for (int i = 0; i < npts; i++) + res = GetValue (elnr, &xref[i*sxref], &x[i*sx], &dxdxref[i*sdxdxref], &values[i*svalues]); + return res; + } + + + + virtual bool GetSurfValue (int /* selnr */, int facetnr, + double /* lam1 */, double /* lam2 */, + double * /* values */) + { + return false; + } + + + virtual bool GetSurfValue (int selnr, int facetnr, + const double xref[], const double x[], const double dxdxref[], + double * values) + { + return GetSurfValue (selnr, facetnr, xref[0], xref[1], values); + } + + + virtual bool GetMultiSurfValue (int selnr, int facetnr, int npts, + const double * xref, int sxref, + const double * x, int sx, + const double * dxdxref, int sdxdxref, + double * values, int svalues) + { + bool res = false; + for (int i = 0; i < npts; i++) + res = GetSurfValue (selnr, facetnr, &xref[i*sxref], &x[i*sx], &dxdxref[i*sdxdxref], &values[i*svalues]); + return res; + } + + virtual bool GetSegmentValue (int segnr, double xref, double * values) + { return false; } + + + virtual int GetNumMultiDimComponents () + { + return 1; + } + + void SetMultiDimComponent (int mc) + { + if (mc >= GetNumMultiDimComponents()) mc = GetNumMultiDimComponents()-1; + if (mc < 0) mc = 0; + multidimcomponent = mc; + } + }; + + + class DLL_HEADER MouseEventHandler + { + public: + virtual void DblClick (int elnr, double x, double y, double z) = 0; + }; + + + class DLL_HEADER UserVisualizationObject + { + public: + virtual void Draw () = 0; + }; + + +} + +#endif +