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

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

#include "stdafx.h"
#include "singlePhaseSolver.h"
#include "windows.h"
#include <iostream>
#include <fstream>
#include <vector>
#include <string>
#include <sstream>
#include <cstdlib>
#include <algorithm>
#include <stdexcept>
#include "iterator"
using namespace std;
void Ng2VTKWithData(Ng_Mesh *mesh, std::vector<double> vecField, std::string filename)
{
  std::cout << "Convert Ng_Mesh to vtk legacy file version of 4.2" << std::endl;
  int i, nseg, ne, np, matnum;
  int nodes[3];
  double point[2];
  //	int *EleNodeNum;
  //	double *NodeCoors;

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

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

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

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

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

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

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

  // Write cell data
  vtkFile << "CELL_DATA " << ne << std::endl;
  vtkFile << "SCALARS Material int 1" << std::endl;
  vtkFile << "LOOKUP_TABLE default" << std::endl;
  for (int i = 0; i < vecMater.size(); i++)
  {
    vtkFile << vecMater[i] << std::endl;
  }

  // Write point data
  vtkFile << "POINT_DATA " << np << std::endl;
  vtkFile << "SCALARS Pressure double 1" << std::endl;
  vtkFile << "LOOKUP_TABLE default" << std::endl;
  for (int i = 0; i < vecField.size(); i++)
  {
    vtkFile << vecField[i] << std::endl;
  }
  vtkFile.close();
}

int _tmain(int argc, _TCHAR *argv[])
{
  /************************************************************************/
  /* prepare data for mesh generation                                     */
  /************************************************************************/
  // 多边形外边界，通过指定四条边来定义
  CBoundShape shape;
  shape.boundShape = POLYGON;
  shape.nNumSegs = 4;
  Segment seg1;
  seg1.p1 = Point(-1500, -1500);
  seg1.p2 = Point(1500, -1500);
  Segment seg2;
  seg2.p1 = Point(1500, -1500);
  seg2.p2 = Point(1500, 1500);
  Segment seg3;
  seg3.p1 = Point(1500, 1500);
  seg3.p2 = Point(-1500, 1500);
  Segment seg4;
  seg4.p1 = Point(-1500, 1500);
  seg4.p2 = Point(-1500, -1500);

  shape.vecSegments.push_back(seg1);
  shape.vecSegments.push_back(seg2);
  shape.vecSegments.push_back(seg3);
  shape.vecSegments.push_back(seg4);

  // shape.cCenter = Point(0, 0);
  // shape.dRadius = 1000;

  CBoundWell w1, w2;
  // 第一口井为直井（默认为直井）
  w1.cCenter = Point(0, 0); // 井中心坐标(0,0)
  w1.dRadius = 0.09144;     // 井半径0.09144

  // 第二口井为垂直裂缝井
  w2.wellType = VERTICAL_FRACTURED;
  w2.dHf = 20;               // 裂缝半长20米
  w2.dTheata = 30;           // 裂缝角度30°
  w2.cCenter = Point(50, 0); // 井的中心位置（50，0）
  w2.dRadius = 0.1;          // 井的半径0.1

  // 第三口井为多段压裂水平井
  CBoundWell w3;
  w3.wellType = MULTI_FRACED_HORIZONTAL;
  w3.cCenter = Point(142, -226); // 中心点位置
  w3.dRadius = 0.1;              // 半径
  w3.dBeta = 20;                 // 水平井角度
  w3.dLength = 200;              // 水平井长度
  w3.nFracNum = 10;              // 裂缝数量
  w3.dHf = 20;                   // 裂缝半长

  //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  // 在该参数设置条件下，一次只能添加一口井进行计算求解
  //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  std::vector<CBoundWell> wells;
  // wells.push_back(w1);
  // wells.push_back(w2);
  wells.push_back(w3);

  CBoundPolygon limit;
  limit.boundShape = POLYGON;
  limit.nNumSegs = 5;
  // 通过指定的边来定义多边形
  limit.vecSegments.push_back(Segment(Point(-516.2335347838829, -58.37577524138442), Point(-592.72401802758, -180.9965861656542)));
  limit.vecSegments.push_back(Segment(Point(-592.72401802758, -180.9965861656542), Point(-485.12090993338666, -280.2195626173259)));
  limit.vecSegments.push_back(Segment(Point(-485.12090993338666, -280.2195626173259), Point(-328.9162066043042, -229.04373034669038)));
  limit.vecSegments.push_back(Segment(Point(-328.9162066043042, -229.04373034669038), Point(-376.36137535544276, -90.33591813134171)));
  limit.vecSegments.push_back(Segment(Point(-376.36137535544276, -90.33591813134171), Point(-516.2335347838829, -58.37577524138442)));
  limit.bAnchor = true;

  // 通过指定的点来定义多边形，逆时针排列
  std::vector<Point> points;
  points.emplace_back(Point(-9.23023, 31.1986));
  points.emplace_back(Point(-27.0886, 11.26));
  points.emplace_back(Point(-23.5281, -14.9567));
  points.emplace_back(Point(10.0862, -27.2422));
  points.emplace_back(Point(39.6129, -1.9348));
  points.emplace_back(Point(32.122, 15.5807));
  points.emplace_back(Point(22.0438, 31.371));
  CBoundPolygon limit1(points);

  std::vector<Point> points2;
  points2.emplace_back(Point(250.074, -101.492));
  points2.emplace_back(Point(-44, -204.034));
  points2.emplace_back(Point(7, -336.098));
  points2.emplace_back(Point(312.333, -234.724));
  CBoundPolygon limit2(points2);

  std::vector<CBoundPolygon> limits;
  limits.push_back(limit);
  limits.push_back(limit1);
  limits.push_back(limit2);

  std::vector<Frac> faults;
  faults.push_back(Frac(Point(-267.3, 120), Point(316.7, 150)));
  faults.push_back(Frac(Point(316.7, 150), Point(316.7, -50)));

  std::vector<Frac> fracs;
  fracs.push_back(Frac(Point(-326.159, 460.059), Point(165.83631, 459.085)));
  fracs.push_back(Frac(Point(50, -350), Point(358.840, -289.811)));

  /************************************************************************/
  /* load netgen mesh generation library and generate mesh                */
  /************************************************************************/
  HMODULE hMod = LoadLibrary("meshgene.dll");
  if (NULL == hMod)
  {
    std::cout << "meshgene.dll加载失败！\n";
    return 1;
  }

  typedef Ng_Mesh *(*NgMeshGen2DByIn2d)(std::vector<CBoundWell> &wells, const std::vector<CBoundPolygon> limits, const std::vector<Frac> &faults, const std::vector<Frac> &fracs, const CBoundShape &shape);
  typedef Ng_Mesh *(*NgMeshGen2D)(std::string strFileIn2D);
  typedef void *(*Ng2VTK)(Ng_Mesh *mesh);
  // Get the address of the Ng_Mesh function
  NgMeshGen2DByIn2d meshgene1 = (NgMeshGen2DByIn2d)GetProcAddress(hMod, "NgMeshGen2DByIn2d");
  NgMeshGen2D meshgene2 = (NgMeshGen2D)GetProcAddress(hMod, "NgMeshGen2D");
  if (NULL == meshgene1)
  {
    FreeLibrary(hMod);
    std::cout << "meshgene文件加载函数地址获取失败！\n";
    return 1;
  }
  Ng_Mesh *mesh = meshgene1(wells, limits, faults, fracs, shape);

  Ng2VTK ng = (Ng2VTK)GetProcAddress(hMod, "Ng2VTK");
  ng(mesh);

  /************************************************************************/
  /* solve the equation by calling the singlephasesolver                 */
  /************************************************************************/
  // Prepare the reservoir base data
  double pBaseData[8];
  pBaseData[0] = 30.47;
  pBaseData[1] = 15.1454;
  pBaseData[2] = 30.48;
  pBaseData[3] = 1.5;
  pBaseData[4] = 1;
  pBaseData[5] = 0.23;
  pBaseData[6] = 0.00043511;
  pBaseData[7] = 1;

  // prepare the well data
  // well w1
  //	w1.dRc=0.01;
  wells[0].dSkin = 0;
  wells[0].dC = 0.0232056;

  wells[0].bisMultFlow = true;
  wells[0].nTimeNumQ = 5;
  wells[0].pdTimeQ[0] = 12;
  wells[0].pdQ[0] = 158.987;
  wells[0].pdTimeQ[1] = 12;
  wells[0].pdQ[1] = 190.785;
  wells[0].pdTimeQ[2] = 12;
  wells[0].pdQ[2] = 222.582;
  wells[0].pdTimeQ[3] = 48;
  wells[0].pdQ[3] = 238.481;
  wells[0].pdTimeQ[4] = 72;
  wells[0].pdQ[4] = 0;

  // 裂缝井必须指定裂缝的导流能力，通过wFrac变量设置
  wells[0].wFrac[0].dFc = 1000;
  wells[0].wFrac[1].dFc = 1000;

  // 多段压裂水平井必须指定裂缝的导流能力1000
  // 导流能力的设置需要在网格划分后完成
  for (int i = 0; i < wells[0].vecFracs.size(); ++i)
  {
    wells[0].vecFracs[i].dFc = 1000; //
  }

  // prepare the data for limits
  limits[0].dComKr = 1;
  limits[0].dComW = 1;

  HMODULE hMod_solver = LoadLibrary("singlePhaseSolverDll.dll");
  if (NULL == hMod_solver)
  {
    std::cout << "singlePhaseSolverDll.dll加载失败！\n";
    return 1;
  }
  TimeDis timeDis;
  timeDis.dTB = 1.0e-5; //

  typedef int (*Solver)(const TimeDis time, double *m_pBaseData, nglib::Ng_Mesh *mesh, std::vector<CBoundWell> &m_wWell, const std::vector<CBoundPolygon> &comPolygon, const std::vector<Frac> &faultPolygon, std::vector<Frac> &fracs, CBoundShape &outBound, std::vector<std::pair<double, std::vector<double>>> &vecNodePre);
  Solver solverFun = (Solver)GetProcAddress(hMod_solver, "singlePhaseSolverNgmesh");

  if (NULL == solverFun)
  {
    FreeLibrary(hMod_solver);
    std::cout << "singlePhaseSolverNgmesh failed!\n";
    return 1;
  }
  //	std::vector<std::vector<Point>> vecBP;
  std::vector<std::pair<double, std::vector<double>>> vecFieldPre;
  solverFun(timeDis, pBaseData, mesh, wells, limits, faults, fracs, shape, vecFieldPre);

  // export bottom hole pressure history
  std::cout << "Export log-log plot and history plot data to files\n";
  std::ofstream historyFile;
  historyFile.precision(8);
  historyFile.open("output/history.txt");
  if (!historyFile.is_open())
  {
    std::cerr << "Error opening file: " << "history.txt" << std::endl;
    return 1;
  }
  int nWell = 0;
  for (int i = 0; i < wells[nWell].m_PressureData.size(); i++)
  {
    historyFile << wells[nWell].m_PressureData[i].x << "\t" << wells[nWell].m_PressureData[i].y << std::endl;
  }
  historyFile.close();

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

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

  std::vector<Point> logPre;
  preLogFun(wells[0].m_PressureData, 5, wells[0].pdTimeQ, wells[0].pdQ, wells[0].nTimeNumQ, logPre);

  // export the pressure derivative
  std::ofstream logFile;
  logFile.precision(8);
  logFile.open("output/logPre.txt");
  if (!logFile.is_open())
  {
    std::cerr << "Error opening file: " << "logPre.txt" << std::endl;
    return 1;
  }
  for (int i = 0; i < logPre.size(); i++)
  {
    logFile << logPre[i].x << "\t" << logPre[i].y << "\t" << logPre[i].z << "\t" << std::endl;
  }
  logFile.close();

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

  FreeLibrary(hMod);
  FreeLibrary(hMod_solver);
}