chore:忽略ML/Training下的Debug和Release文件夹

main
lh 3 weeks ago
parent 5b948c52d2
commit 99e74efa8c

3
.gitignore vendored

@ -119,4 +119,5 @@ ML/nmWTAI-ML/data
ML/nmWTAI-ML/models
ML/nmWTAI-ML/results
__pycache__
ML/Training
ML/Training/Debug
ML/Training/Release

@ -0,0 +1,196 @@
#include "RunnerIO.h"
#include <fstream>
#include <iostream>
#include <stdint.h>
#include <algorithm>
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<char> 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();
}
}

@ -0,0 +1,48 @@
#pragma once
#include <string>
#include <vector>
#include <cstdint>
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<double> timeQ; // hours
std::vector<double> 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<double> t;
std::vector< std::vector<double> > pw;
std::vector< std::vector<double> > loglog_t;
std::vector< std::vector<double> > loglog_p;
std::vector< std::vector<double> > 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);
}

@ -0,0 +1,444 @@
#include "singlePhaseSolver.h"
#include <Windows.h>
#include <iostream>
#include <string>
#include <vector>
#include <fstream>
#include <sstream>
#include <algorithm>
#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<Point>&, const int&, double*, double*, int, std::vector<Point>&);
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 / wellboreCfill 更快)
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/Cfresize + 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<double>());
rr.loglog_p.assign(rr.nWells, std::vector<double>());
rr.loglog_deriv.assign(rr.nWells, std::vector<double>());
if (!preLogFun) {
return true; // 没 dllloglog 为空,但结构正确
}
for (unsigned int w = 0; w < rr.nWells; ++w)
{
std::vector<Point> 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<double> timeQ;
std::vector<double> 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<Point> 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<double>& lt = rr.loglog_t[w];
std::vector<double>& lp = rr.loglog_p[w];
std::vector<double>& 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 <dataset.bin> <HX_NWTM.dll> <license.dat>\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;
}

@ -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 bytesUTF-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<double>& 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<double>& 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<int>& 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<int>& 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<std::vector<double> >& 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<std::vector<double> >& 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<std::vector<std::vector<double> > >& 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<std::vector<std::vector<double> > >& 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<vector<double>>
// 为了兼容,这里存原始数据,不做转换
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;
}

@ -0,0 +1,127 @@
#pragma once
#ifndef DATASET_IO_H
#define DATASET_IO_H
#include <string>
#include <vector>
#include <fstream>
#include <iostream>
#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<double>& v);
static bool readStdVecD(std::ifstream& fs, std::vector<double>& v);
static bool writeStdVecI(std::ofstream& fs, const std::vector<int>& v);
static bool readStdVecI(std::ifstream& fs, std::vector<int>& v);
static bool writeStdVec2D(std::ofstream& fs, const std::vector<std::vector<double> >& v);
static bool readStdVec2D(std::ifstream& fs, std::vector<std::vector<double> >& v);
static bool writeStdVec3D(std::ofstream& fs, const std::vector<std::vector<std::vector<double> > >& v);
static bool readStdVec3D(std::ifstream& fs, std::vector<std::vector<std::vector<double> > >& 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

@ -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;
}

@ -0,0 +1,101 @@
#pragma once
#ifndef GRID_CACHE_IO_H
#define GRID_CACHE_IO_H
#include <string>
#include <fstream>
#include <iostream>
#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

@ -0,0 +1,226 @@
#include "SceneIO.h"
#include <iostream>
#include <stdint.h>
bool SceneIO::readInt(std::ifstream& fs, int& value)
{
fs.read(reinterpret_cast<char*>(&value), sizeof(int));
return fs.good();
}
bool SceneIO::readDouble(std::ifstream& fs, double& value)
{
fs.read(reinterpret_cast<char*>(&value), sizeof(double));
return fs.good();
}
bool SceneIO::readSafeSize(std::ifstream& fs, size_t& size)
{
uint32_t sz = 0;
fs.read(reinterpret_cast<char*>(&sz), sizeof(uint32_t));
if (!fs.good()) return false;
const uint32_t MAX_SIZE = 100000000;
if (sz > MAX_SIZE) return false;
size = static_cast<size_t>(sz);
return true;
}
bool SceneIO::readString(std::ifstream& fs, std::string& value)
{
uint32_t byteLen = 0;
fs.read(reinterpret_cast<char*>(&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<uint16_t> utf16(charCount);
if (charCount > 0) {
fs.read(reinterpret_cast<char*>(&utf16[0]), byteLen);
if (!fs.good()) return false;
}
// 简化转换:只取低字节(你的数据目前可用)
value.clear();
value.reserve(charCount);
for (size_t i = 0; i < charCount; ++i) {
value += static_cast<char>(utf16[i] & 0xFF);
}
return true;
}
bool SceneIO::readVector1D(std::ifstream& fs, std::vector<double>& vec)
{
size_t sz = 0;
if (!readSafeSize(fs, sz)) return false;
vec.resize(sz);
if (sz > 0) {
fs.read(reinterpret_cast<char*>(&vec[0]), sz * sizeof(double));
if (!fs.good()) return false;
}
return true;
}
bool SceneIO::readVector2D(std::ifstream& fs, std::vector<std::vector<double> >& 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<std::vector<std::vector<double> > >& 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<int>& 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<char*>(&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;
}

@ -0,0 +1,84 @@
#pragma once
#include <vector>
#include <string>
#include <fstream>
// 场景结构体
struct PebiScene {
int version;
// 网格几何
int D;
double GridControl;
std::vector<std::vector<double> > Boundary;
std::vector<std::vector<double> > VerticalWell;
std::vector<std::vector<double> > HorizontalWell;
std::vector<std::vector<double> > FractureVerticalWell;
std::vector<std::vector<std::vector<double> > > MultistageFracturedHorizontalWell;
std::vector<std::vector<double> > InclinedWell;
std::vector<std::vector<double> > Fault;
// 井顺序
std::vector<int> wellType;
std::vector<std::string> wellName;
// 求解器参数
int solverType;
// Rate 流量数据
struct {
std::vector<std::vector<double> > t;
std::vector<std::vector<double> > qo;
std::vector<std::vector<double> > qg;
std::vector<std::vector<double> > qw;
} Rate;
// CS 井筒参数
struct {
std::vector<double> C;
std::vector<double> S;
} CS;
// 每口井的流量段索引
std::vector<int> wellFlowSectionIndex;
// PVT 数据
struct {
std::vector<double> p;
double pb;
std::vector<double> Rso, Bo, Co, miuo, rouo;
std::vector<double> Rv, Bg, Cg, miug, roug, Z;
std::vector<double> Rsw, Bw, Cw, miuw, rouw;
std::vector<double> V, k_kinitial, Cf_Cfinitial;
std::vector<double> 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<double>& vec);
static bool readVector2D(std::ifstream& fs, std::vector<std::vector<double> >& vec);
static bool readVector3D(std::ifstream& fs, std::vector<std::vector<std::vector<double> > >& vec);
// 读取整数向量
static bool readStdVecI(std::ifstream& fs, std::vector<int>& vec);
};

@ -0,0 +1,318 @@
#include <Windows.h>
#include <iostream>
#include <fstream>
#include <string>
#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) {
// 真正加载 datasetstrictVerify=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;
}

@ -0,0 +1,4 @@
#pragma once
#define WIN32_LEAN_AND_MEAN // 从 Windows 头文件中排除极少使用的内容
// Windows 头文件
#include <windows.h>

@ -0,0 +1,167 @@
// netgen网格生成库的C++接口,传入边界、井筒、裂缝、断层和复合区,生成二维三角形网格。
// 作者:万义钊
// 日期2025年1月28日除夕
// 版本1.0
#pragma once
namespace nglib
{
#include "nglib.h"
}
using namespace nglib;
#include <iostream>
#include <string>
#include <vector>
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<double> pointData;
Point(double x0, double y0)
{
x = x0;
y = y0;
}
Point() {};
// 拷贝构造函数
Point(const Point &other)
{
// std::cout << "Copy constructor called" << std::endl;
x = other.x; // 复制数据
y = other.y;
z = other.z;
pointData = other.pointData;
}
};
struct Cell
{
int type;
std::vector<int> pointIndices;
std::vector<double> 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<Point> 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<Segment> vecSegments; // 边
// 表示圆形
Point cCenter;
double dRadius;
int nNodeNum;
BOUND_TYPE boundType;
};
// 多边形,内部的复合区域
class CBoundPolygon : public CBoundShape
{
public:
int nNumSegs; // 边的数量
std::vector<Segment> 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<Frac> vecFracs;//压裂裂缝
};
#ifdef __cplusplus
extern "C"
{
#endif
__declspec(dllexport) Ng_Mesh *NgMeshGen2D(std::string strFileIn2D);
__declspec(dllexport) Ng_Mesh *NgMeshGen2DByIn2d(std::vector<CBoundWell> wells, std::vector<CBoundPolygon> limits, std::vector<CBoundPolygon> faults, std::vector<Frac> fracs, CBoundShape shape);
__declspec(dllexport) void Ng2VTK(Ng_Mesh *mesh);
#ifdef __cplusplus
}
#endif

@ -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 <iostream>
namespace netgen
{
DLL_HEADER extern std::ostream * testout;
DLL_HEADER extern int printmessage_importance;
}
#endif

@ -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.
<i>Note:</i>
- 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.
<i>Note:</i>
- 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.
<b>Note</b>: This function only limits the <b>Maximum</b>
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 <i>double</i>, 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.
<b>Note</b>: This function only limits the <b>Maximum</b> 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 <i>double</i>, 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 <i>double</i>, 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.
<b>Note</b>: This function only limits the <b>Maximum</b> 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 <i>double</i>, 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 <i>double</i>, 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 <i>double</i>, 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.
<b>Note</b>: 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

@ -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 <iostream>
#include <vector>
#include <cmath>
#include <algorithm>
#include <fstream>
#include <ctime>
#include <iomanip>
#include <unordered_set>
#include <Windows.h>
//const double M_PI = acos(-1.0);
typedef std::vector<std::vector<std::vector<double>>>dVec3; //三维数组:double
typedef std::vector<std::vector<double>>dVec2; //二维数组:double
typedef std::vector<std::vector<int>>iVec2; //二维数组:int
typedef std::vector<double>dVec1; //一维数组:double
typedef std::vector<int>iVec1; //一维数组:int
template<class T> void HX_copy(std::vector<std::vector<std::vector<T>>>& p1, const std::vector<std::vector<std::vector<T>>>& 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<class T> void HX_copy(std::vector<std::vector<T>>& p1, const std::vector<std::vector<T>>& 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<class T> void HX_copy(std::vector<T>& p1, const std::vector<T>& p0)
{
int m = p0.size(); p1.resize(m);
for (int i = 0; i < m; ++i)
{
p1[i] = p0[i];
}
}
template<class T> void HX_copy(std::vector<std::vector<std::vector<T>>>& 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<class T> void HX_copy(std::vector<std::vector<T>>& 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<class T> void HX_copy(std::vector<T>& 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<point> 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); //数值试井模型求解器接口

@ -0,0 +1,396 @@
#pragma once
#include <string>
#include <vector>
#include <unordered_set>
#include <algorithm>
namespace nglib {
#include "nglib.h"
}
using namespace nglib;
enum BOUND_TYPE {
CLOSED,
CONST_PRESSURE,
CONST_RATE,
};
enum BOUND_SHAPE {
CIRCLE,
RECTANGLE,
POLYGON
};
enum WELL_TYPE {
VERTICAL,
VERTICAL_FRACTURED,
MULTI_FRACED_HORIZONTAL
};
struct Point {
double x, y, z;
std::vector<double> pointData;
Point(double x0, double y0) {
x = x0;
y = y0;
}
Point() {};
// 拷贝构造函数
Point(const Point& other) {
//std::cout << "Copy constructor called" << std::endl;
x = other.x; // 复制数据
y = other.y;
z = other.z;
pointData = other.pointData;
}
bool operator==(const Point& other) const {
const double epsilon = 1e-4; // 误差范围
return fabs(x - other.x) < epsilon &&
fabs(y - other.y) < epsilon;
}
};
struct TimeDis {
int nLogNum;//每个对数周期的时间点数
double dTB;//起始时间
TimeDis() {
nLogNum = 10; //默认为每个对数周期10个点
dTB = 1.0e-4; //默认初始时刻为1e-4,如果计算的双对数曲线的井储段未显示,则需要减小该值
}
};
class Cell {
public:
int type;
std::vector<int> pointIndices;
std::vector<double> cellData;
Cell(const Cell& other) {
type = other.type;
pointIndices = other.pointIndices;
cellData = other.cellData;
}
Cell() {}
bool operator==(const Cell&other)const {
if (type != other.type) {
return false;
}
if (pointIndices.size() != other.pointIndices.size()) {
return false;
}
std::vector<int> a_sorted(pointIndices);
std::vector<int> b_sorted(other.pointIndices);
std::sort(a_sorted.begin(), a_sorted.end());
std::sort(b_sorted.begin(), b_sorted.end());
return a_sorted == b_sorted;
}
};
class hash_fun {
public:
int operator()(const Cell &A) const {
return A.type;
}
};
/*
class CBoundLine
{
// DECLARE_SERIAL(CBoundLine)
public:
CBoundLine() {
nNodeNum = 8;
nboundtype = 0;
ddl = 100;
bisSel = false;
};
virtual ~CBoundLine() {};
Point pPre, pOri;
int nboundtype;
int nNodeNum;
double ddl;
int nNodeMin;
int nNodeMax;
bool bisSel;
};*/
class Line {
public:
Point p1;
Point p2;
};
// 裂缝线
class Frac : public Line {
public:
Frac(Point p1, Point p2) {
this->p1 = p1;
this->p2 = p2;
dE = 0.01;
dW = 1;
dFc = 1000;
dKr = dFc / dE;
}
Frac() {
dE = 0.01;
dW = 1;
dFc = 1000;
dKr = dFc / dE;
}
// 交点
int nIntersect; // 交点数量
std::vector<Point> vector_interPoint; // 交点坐标
//裂缝上的线单元
std::vector<Cell> fracCells;
double dE;//裂缝宽度
double dKr;//裂缝渗透率比,Kf/K
double dW;//裂缝储能比,(phi*Ct)f/(phi*Ct)
double dFc;//裂缝导流能力,dFc=dKr*dE;
//记录netgen的in2d文件中几何点的编号
//一个裂缝只需要2个节点即可记录
std::vector<int> ngPointIndex;
};
// 定义外边界的线
class Segment : public Line {
public:
Segment(Point p1, Point p2) {
this->p1 = p1;
this->p2 = p2;
}
Segment() {}
BOUND_TYPE boundType; // 边界线的边界条件类型
double dConstPressure;//定义边界时的压力值
int nNumNodes; // 边界线上的网格点数
std::vector<int> vecNodes;//边界线上的网格节点编号
};
// 外边界的基类
class CBoundShape {
public:
BOUND_SHAPE boundShape; // 定义外边界类型0为圆形1为四边形2为多边形
// 表示多边形(包括四边形)的参数
int nNumSegs; // 边的数量
std::vector<Segment> vecSegments; // 边
// 表示圆形
Point cCenter;
double dRadius;
int nNodeNum;
BOUND_TYPE boundType;
double dConstPressure;//圆形定压边界时的压力值
std::vector<int> vecNodes;//圆形边界上的网格点编号
//记录netgen的in2d文件中几何点的编号
std::vector<int> ngPointIndex;
};
// 多边形,内部的复合区域
class CBoundPolygon : public CBoundShape {
public:
CBoundPolygon() {
bAnchor = false;
}
CBoundPolygon(std::vector<Point> points) { //按拟时针给定的点来构建多边形
for (int i = 0; i < points.size() - 1; ++i) {
vecSegments.push_back(Segment(points[i], points[i + 1]));
}
vecSegments.push_back(Segment(points[points.size() - 1], points[0]));
nNumSegs = points.size();
bAnchor = false;
}
int nNumSegs; // 边的数量
std::vector<Segment> vecSegments; // 边
double dComKr;//复合区的渗透率比
double dComW;//复合区的储能比
bool bAnchor;//是否作为复合区进行材料属性的分区赋值false表示不进行复合区的材料编号赋值即其材料编号与背景值相同
//记录netgen的in2d文件中几何点的编号
std::vector<int> ngPointIndex;
};
class CBoundWell {
// DECLARE_SERIAL(CWell)
public:
CBoundWell() {
nNodeNum = 12;
// dRc = 0.1;
nFlowType = 0;
nTimeNumQ = 0;
nTimeNumP = 0;
dSkin = 0;
dC = 0.01;
dMinPress = 0.101325;
bisSel = false;
nDenseNum = 10;
dDenseRadius[0] = 1;
nDenseNodeNUm[0] = 18;
dDenseRadius[1] = 10;
nDenseNodeNUm[1] = 23;
dDenseRadius[2] = 50;
nDenseNodeNUm[2] = 22;
dDenseRadius[3] = 100;
nDenseNodeNUm[3] = 25;
dDenseRadius[4] = 500;
nDenseNodeNUm[4] = 26;
dDenseRadius[5] = 1000;
nDenseNodeNUm[5] = 27;
dDenseRadius[6] = 2000;
nDenseNodeNUm[6] = 32;
dDenseRadius[7] = 3000;
nDenseNodeNUm[7] = 34;
dDenseRadius[8] = 4000;
nDenseNodeNUm[8] = 35;
dDenseRadius[9] = 5000;
nDenseNodeNUm[9] = 47;
dDenseRadius[10] = 10000;
nDenseNodeNUm[10] = 48;
dDenseRadius[11] = 20000;
nDenseNodeNUm[11] = 48;
bisAutoChoose = false;
bisMultFlow = true;
nInComNote = 0;
wellType = VERTICAL;
dWidth = 0.1;
};
virtual ~CBoundWell() {};
void SetDenseValue();
public:
//void Serialize(CArchive& ar);
std::string sWellname;
Point pCenter;
int nFlowType;
int nNodeNum;
//double dRc;
double dSkin;
double dC;
double dMinPress;
int nTimeNumQ;
int nTimeNumP;
double pdTimeQ[100];
double pdTimeP[100];
double pdQ[100];
double pdP[100];
int nNodeMin;
int nNodeMax;
bool bisSel;
int nDenseNum;
double dDenseRadius[20];
int nDenseNodeNUm[20];
double dS;
bool bisAutoChoose;
bool bisMultFlow;
int nInComNote;
std::vector<int> m_NodeIndex;//all the nodes index of the wellbore
std::vector<Point> m_PressureData;//计算的井筒压力数据
WELL_TYPE wellType;
//直井的几何参数wellType=VERTICAL
Point cCenter;
double dRadius;
//int nNodeNum;
BOUND_TYPE boundType;
//直井压裂井独有的参数wellType=VERTICAL_FRACTURED
//直井的中心和半径仍然使用cCenter和dRaduis表示
double dHf;//裂缝半长
double dTheata;//裂缝角度(相对于x轴)
Frac wFrac[2];//通过dHf和dTheata计算的裂缝存储在该变量中用于后续的计算
//多段压裂水平井的参数cCenter表示水平井中心点坐标wellType=MULTI_FRACED_HORIZONTAL
double dLength;//长度
double dBeta;//角度
int nFracNum;//压裂裂缝的数量
std::vector<Frac> vecFracs;//压裂裂缝
double dWidth;//2d模型中水平井宽度
CBoundPolygon wellPoly;//水平井的边界形状
//记录netgen的in2d文件中几何点的编号
std::vector<int> ngPointIndex;
};
/*
class CBoundCir
{
// DECLARE_SERIAL(CBoundCir)
public:
CBoundCir() {
nNodeNum = 24;
nboundtype = 0;
bisSel = false;
dComKr = 10;
dComW = 3;
};
virtual ~CBoundCir() {};
public:
Point pCenter;
double dRc;
int nboundtype;
int nNodeNum;
bool bisSel;
int nNodeMin;
int nNodeMax;
int nMeshMin;
int nMeshMax;
double dComKr;
double dComW;
};*/
#ifdef __cplusplus
extern "C" {
#endif
//mesh generation API
__declspec(dllexport) Ng_Mesh* NgMeshGen2D(std::string strFileIn2D);
__declspec(dllexport) Ng_Mesh* NgMeshGen2DByIn2d(std::vector<CBoundWell>& wells, std::vector<CBoundPolygon> limits,
std::vector<Frac>& faults, std::vector<Frac>& fracs, CBoundShape &shape);
__declspec(dllexport) void Ng2VTK(Ng_Mesh *mesh);
//singlePhase solver API
__declspec(dllexport) int singlePhaseSolver(double* pBaseData, std::vector<Point> points, std::vector<Cell> cells,
int m_nWellNum, CBoundWell* pwells, std::vector<std::vector<Point>>& vecBP,
std::vector<std::pair<double, std::vector<double>>>& vecNodePre);
__declspec(dllexport) int singlePhaseSolverNgmesh(const TimeDis time, double* m_pBaseData, nglib::Ng_Mesh* mesh,
std::vector<CBoundWell>& m_wWell, const std::vector<CBoundPolygon>& comPolygon, const std::vector<Frac>& faultPolygon,
std::vector<Frac>& fracs, CBoundShape &outBound, std::vector<std::pair<double, std::vector<double>>>& vecNodePre,
int& timeStep, int& totalSteps);
__declspec(dllexport) int singlePhaseGasSolverNgmesh(double* component, double temperature, const TimeDis time,
double* m_pBaseData, nglib::Ng_Mesh* mesh, std::vector<CBoundWell>& m_wWell,
const std::vector<CBoundPolygon>& comPolygon, const std::vector<Frac>& faultPolygon, std::vector<Frac>& fracs,
CBoundShape &outBound, std::vector<std::pair<double, std::vector<double>>>& vecNodePre, int& timeStep, int& totalSteps);
__declspec(dllexport) bool logLogPre(const std::vector<Point>& pressure, const int& nSection, double* pdTimeQ,
double* pdQ, int nTimeQ, std::vector<Point>& logPre);
#ifdef __cplusplus
}
#endif

@ -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
Loading…
Cancel
Save