You cannot select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
775 lines
26 KiB
C++
775 lines
26 KiB
C++
#include "nmPebiRegionPressureCalculator.h"
|
|
|
|
#include "nmPebiResultSnapshot.h"
|
|
|
|
#include <QRectF>
|
|
|
|
#include <vtkCell.h>
|
|
#include <vtkCellType.h>
|
|
#include <vtkDoubleArray.h>
|
|
#include <vtkPoints.h>
|
|
#include <vtkUnstructuredGrid.h>
|
|
|
|
#include <algorithm>
|
|
#include <float.h>
|
|
#include <math.h>
|
|
|
|
#if !defined(_MSC_VER)
|
|
#include <cmath>
|
|
#endif
|
|
|
|
namespace {
|
|
|
|
// 所有区域统计几何均投影到 XY 平面，三角形是剪裁和求交的最小单元。
|
|
struct Triangle2D
|
|
{
|
|
QPointF m_a;
|
|
QPointF m_b;
|
|
QPointF m_c;
|
|
};
|
|
|
|
struct PreparedPolygon
|
|
{
|
|
QVector<QPointF> m_vecPoints;
|
|
QVector<Triangle2D> m_vecTriangles;
|
|
QRectF m_oBounds;
|
|
double m_dArea;
|
|
};
|
|
|
|
struct CachedCell
|
|
{
|
|
CachedCell()
|
|
: m_nCellId(-1),
|
|
m_dArea(0.0),
|
|
m_dPoreVolumeWeight(0.0),
|
|
m_bWeightValid(false)
|
|
{
|
|
}
|
|
|
|
vtkIdType m_nCellId;
|
|
double m_dArea;
|
|
double m_dPoreVolumeWeight;
|
|
bool m_bWeightValid;
|
|
};
|
|
|
|
bool isFiniteValue(double dValue)
|
|
{
|
|
#if defined(_MSC_VER)
|
|
return _finite(dValue) != 0;
|
|
#else
|
|
return std::isfinite(dValue);
|
|
#endif
|
|
}
|
|
|
|
double crossProduct(const QPointF& a,
|
|
const QPointF& b,
|
|
const QPointF& c)
|
|
{
|
|
return (b.x() - a.x()) * (c.y() - a.y()) -
|
|
(b.y() - a.y()) * (c.x() - a.x());
|
|
}
|
|
|
|
double squaredDistance(const QPointF& a, const QPointF& b)
|
|
{
|
|
const double dX = a.x() - b.x();
|
|
const double dY = a.y() - b.y();
|
|
return dX * dX + dY * dY;
|
|
}
|
|
|
|
double signedDoubleArea(const QVector<QPointF>& vecPoints)
|
|
{
|
|
if(vecPoints.size() < 3) {
|
|
return 0.0;
|
|
}
|
|
|
|
// 先平移到首点附近再做鞋带求和，避免大地坐标中的大数相减吞掉小区域面积。
|
|
const QPointF oOrigin = vecPoints.first();
|
|
double dArea = 0.0;
|
|
for(int nIndex = 0; nIndex < vecPoints.size(); ++nIndex) {
|
|
const QPointF& a = vecPoints[nIndex];
|
|
const QPointF& b = vecPoints[(nIndex + 1) % vecPoints.size()];
|
|
const double dAx = a.x() - oOrigin.x();
|
|
const double dAy = a.y() - oOrigin.y();
|
|
const double dBx = b.x() - oOrigin.x();
|
|
const double dBy = b.y() - oOrigin.y();
|
|
dArea += dAx * dBy - dBx * dAy;
|
|
}
|
|
return dArea;
|
|
}
|
|
|
|
QRectF polygonBounds(const QVector<QPointF>& vecPoints)
|
|
{
|
|
double dMinimumX = vecPoints[0].x();
|
|
double dMaximumX = dMinimumX;
|
|
double dMinimumY = vecPoints[0].y();
|
|
double dMaximumY = dMinimumY;
|
|
for(int nIndex = 1; nIndex < vecPoints.size(); ++nIndex) {
|
|
dMinimumX = qMin(dMinimumX, vecPoints[nIndex].x());
|
|
dMaximumX = qMax(dMaximumX, vecPoints[nIndex].x());
|
|
dMinimumY = qMin(dMinimumY, vecPoints[nIndex].y());
|
|
dMaximumY = qMax(dMaximumY, vecPoints[nIndex].y());
|
|
}
|
|
return QRectF(QPointF(dMinimumX, dMinimumY),
|
|
QPointF(dMaximumX, dMaximumY));
|
|
}
|
|
|
|
bool pointOnSegment(const QPointF& p,
|
|
const QPointF& a,
|
|
const QPointF& b,
|
|
double dLinearTolerance,
|
|
double dAreaTolerance)
|
|
{
|
|
if(qAbs(crossProduct(a, b, p)) > dAreaTolerance) {
|
|
return false;
|
|
}
|
|
return p.x() >= qMin(a.x(), b.x()) - dLinearTolerance &&
|
|
p.x() <= qMax(a.x(), b.x()) + dLinearTolerance &&
|
|
p.y() >= qMin(a.y(), b.y()) - dLinearTolerance &&
|
|
p.y() <= qMax(a.y(), b.y()) + dLinearTolerance;
|
|
}
|
|
|
|
int orientationSign(const QPointF& a,
|
|
const QPointF& b,
|
|
const QPointF& c,
|
|
double dAreaTolerance)
|
|
{
|
|
const double dCross = crossProduct(a, b, c);
|
|
if(dCross > dAreaTolerance) {
|
|
return 1;
|
|
}
|
|
if(dCross < -dAreaTolerance) {
|
|
return -1;
|
|
}
|
|
return 0;
|
|
}
|
|
|
|
bool segmentsIntersect(const QPointF& a,
|
|
const QPointF& b,
|
|
const QPointF& c,
|
|
const QPointF& d,
|
|
double dLinearTolerance,
|
|
double dAreaTolerance)
|
|
{
|
|
const int nAbC = orientationSign(a, b, c, dAreaTolerance);
|
|
const int nAbD = orientationSign(a, b, d, dAreaTolerance);
|
|
const int nCdA = orientationSign(c, d, a, dAreaTolerance);
|
|
const int nCdB = orientationSign(c, d, b, dAreaTolerance);
|
|
if(nAbC * nAbD < 0 && nCdA * nCdB < 0) {
|
|
return true;
|
|
}
|
|
return (nAbC == 0 && pointOnSegment(c, a, b,
|
|
dLinearTolerance, dAreaTolerance)) ||
|
|
(nAbD == 0 && pointOnSegment(d, a, b,
|
|
dLinearTolerance, dAreaTolerance)) ||
|
|
(nCdA == 0 && pointOnSegment(a, c, d,
|
|
dLinearTolerance, dAreaTolerance)) ||
|
|
(nCdB == 0 && pointOnSegment(b, c, d,
|
|
dLinearTolerance, dAreaTolerance));
|
|
}
|
|
|
|
bool isSimplePolygon(const QVector<QPointF>& vecPoints,
|
|
double dLinearTolerance,
|
|
double dAreaTolerance)
|
|
{
|
|
// 只检查非相邻边；相邻边共享端点是合法多边形的必要条件。
|
|
const int nPointCount = vecPoints.size();
|
|
for(int nFirst = 0; nFirst < nPointCount; ++nFirst) {
|
|
const int nFirstNext = (nFirst + 1) % nPointCount;
|
|
for(int nSecond = nFirst + 1; nSecond < nPointCount; ++nSecond) {
|
|
const int nSecondNext = (nSecond + 1) % nPointCount;
|
|
if(nFirst == nSecond || nFirstNext == nSecond ||
|
|
nSecondNext == nFirst) {
|
|
continue;
|
|
}
|
|
if(segmentsIntersect(vecPoints[nFirst], vecPoints[nFirstNext],
|
|
vecPoints[nSecond], vecPoints[nSecondNext],
|
|
dLinearTolerance, dAreaTolerance)) {
|
|
return false;
|
|
}
|
|
}
|
|
}
|
|
return true;
|
|
}
|
|
|
|
bool pointInsideTriangle(const QPointF& p,
|
|
const QPointF& a,
|
|
const QPointF& b,
|
|
const QPointF& c,
|
|
double dAreaTolerance)
|
|
{
|
|
return crossProduct(a, b, p) >= -dAreaTolerance &&
|
|
crossProduct(b, c, p) >= -dAreaTolerance &&
|
|
crossProduct(c, a, p) >= -dAreaTolerance;
|
|
}
|
|
|
|
bool triangulatePolygon(const QVector<QPointF>& vecPoints,
|
|
double dAreaTolerance,
|
|
QVector<Triangle2D>& vecTriangles)
|
|
{
|
|
// 耳切法同时支持凸、凹简单多边形，无需依赖 VTK 的三维三角化管线。
|
|
vecTriangles.clear();
|
|
QVector<int> vecIndices;
|
|
for(int nIndex = 0; nIndex < vecPoints.size(); ++nIndex) {
|
|
vecIndices.append(nIndex);
|
|
}
|
|
|
|
while(vecIndices.size() > 3) {
|
|
bool bEarFound = false;
|
|
for(int nPosition = 0; nPosition < vecIndices.size(); ++nPosition) {
|
|
const int nPrevious = vecIndices[
|
|
(nPosition + vecIndices.size() - 1) % vecIndices.size()];
|
|
const int nCurrent = vecIndices[nPosition];
|
|
const int nNext = vecIndices[(nPosition + 1) % vecIndices.size()];
|
|
const QPointF& a = vecPoints[nPrevious];
|
|
const QPointF& b = vecPoints[nCurrent];
|
|
const QPointF& c = vecPoints[nNext];
|
|
if(crossProduct(a, b, c) <= dAreaTolerance) {
|
|
continue;
|
|
}
|
|
|
|
bool bContainsOtherPoint = false;
|
|
for(int nOtherPosition = 0;
|
|
nOtherPosition < vecIndices.size(); ++nOtherPosition) {
|
|
const int nOther = vecIndices[nOtherPosition];
|
|
if(nOther == nPrevious || nOther == nCurrent ||
|
|
nOther == nNext) {
|
|
continue;
|
|
}
|
|
if(pointInsideTriangle(vecPoints[nOther], a, b, c,
|
|
dAreaTolerance)) {
|
|
bContainsOtherPoint = true;
|
|
break;
|
|
}
|
|
}
|
|
if(bContainsOtherPoint) {
|
|
continue;
|
|
}
|
|
|
|
Triangle2D oTriangle;
|
|
oTriangle.m_a = a;
|
|
oTriangle.m_b = b;
|
|
oTriangle.m_c = c;
|
|
vecTriangles.append(oTriangle);
|
|
vecIndices.remove(nPosition);
|
|
bEarFound = true;
|
|
break;
|
|
}
|
|
if(!bEarFound) {
|
|
vecTriangles.clear();
|
|
return false;
|
|
}
|
|
}
|
|
|
|
if(vecIndices.size() != 3) {
|
|
vecTriangles.clear();
|
|
return false;
|
|
}
|
|
Triangle2D oTriangle;
|
|
oTriangle.m_a = vecPoints[vecIndices[0]];
|
|
oTriangle.m_b = vecPoints[vecIndices[1]];
|
|
oTriangle.m_c = vecPoints[vecIndices[2]];
|
|
if(crossProduct(oTriangle.m_a, oTriangle.m_b, oTriangle.m_c) <=
|
|
dAreaTolerance) {
|
|
vecTriangles.clear();
|
|
return false;
|
|
}
|
|
vecTriangles.append(oTriangle);
|
|
return true;
|
|
}
|
|
|
|
bool preparePolygon(const QVector<QPointF>& vecInput,
|
|
PreparedPolygon& oPolygon)
|
|
{
|
|
oPolygon = PreparedPolygon();
|
|
if(vecInput.size() < 3) {
|
|
return false;
|
|
}
|
|
for(int nIndex = 0; nIndex < vecInput.size(); ++nIndex) {
|
|
if(!isFiniteValue(vecInput[nIndex].x()) ||
|
|
!isFiniteValue(vecInput[nIndex].y())) {
|
|
return false;
|
|
}
|
|
}
|
|
|
|
// 容差跟随当前多边形尺度，兼容小尺寸网格和大地坐标。
|
|
const QRectF oInputBounds = polygonBounds(vecInput);
|
|
const double dScale = qMax(1.0, qMax(oInputBounds.width(),
|
|
oInputBounds.height()));
|
|
const double dLinearTolerance = dScale * 1.0e-10;
|
|
const double dAreaTolerance = dScale * dScale * 1.0e-12;
|
|
const double dSquaredTolerance =
|
|
dLinearTolerance * dLinearTolerance;
|
|
|
|
// 先去掉相邻重复点及显式闭合点，避免产生零长边。
|
|
QVector<QPointF> vecPoints;
|
|
for(int nIndex = 0; nIndex < vecInput.size(); ++nIndex) {
|
|
if(vecPoints.isEmpty() ||
|
|
squaredDistance(vecPoints.last(), vecInput[nIndex]) >
|
|
dSquaredTolerance) {
|
|
vecPoints.append(vecInput[nIndex]);
|
|
}
|
|
}
|
|
if(vecPoints.size() > 1 &&
|
|
squaredDistance(vecPoints.first(), vecPoints.last()) <=
|
|
dSquaredTolerance) {
|
|
vecPoints.remove(vecPoints.size() - 1);
|
|
}
|
|
|
|
// 反复清理共线中间点，直到点集稳定，降低耳切法退化风险。
|
|
bool bRemovedPoint = true;
|
|
while(bRemovedPoint && vecPoints.size() >= 3) {
|
|
bRemovedPoint = false;
|
|
for(int nIndex = 0; nIndex < vecPoints.size(); ++nIndex) {
|
|
const QPointF& a = vecPoints[
|
|
(nIndex + vecPoints.size() - 1) % vecPoints.size()];
|
|
const QPointF& p = vecPoints[nIndex];
|
|
const QPointF& b = vecPoints[(nIndex + 1) % vecPoints.size()];
|
|
if(pointOnSegment(p, a, b,
|
|
dLinearTolerance, dAreaTolerance)) {
|
|
vecPoints.remove(nIndex);
|
|
bRemovedPoint = true;
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
if(vecPoints.size() < 3 ||
|
|
!isSimplePolygon(vecPoints,
|
|
dLinearTolerance, dAreaTolerance)) {
|
|
return false;
|
|
}
|
|
|
|
double dDoubleArea = signedDoubleArea(vecPoints);
|
|
if(!isFiniteValue(dDoubleArea) ||
|
|
qAbs(dDoubleArea) <= dAreaTolerance) {
|
|
return false;
|
|
}
|
|
// 后续半平面剪裁约定顶点为逆时针顺序。
|
|
if(dDoubleArea < 0.0) {
|
|
std::reverse(vecPoints.begin(), vecPoints.end());
|
|
dDoubleArea = -dDoubleArea;
|
|
}
|
|
|
|
QVector<Triangle2D> vecTriangles;
|
|
if(!triangulatePolygon(vecPoints,
|
|
dAreaTolerance, vecTriangles)) {
|
|
return false;
|
|
}
|
|
oPolygon.m_vecPoints = vecPoints;
|
|
oPolygon.m_vecTriangles = vecTriangles;
|
|
oPolygon.m_oBounds = polygonBounds(vecPoints);
|
|
oPolygon.m_dArea = dDoubleArea * 0.5;
|
|
return isFiniteValue(oPolygon.m_dArea) && oPolygon.m_dArea > 0.0;
|
|
}
|
|
|
|
QPointF segmentLineIntersection(const QPointF& s,
|
|
const QPointF& e,
|
|
double dStartDistance,
|
|
double dEndDistance)
|
|
{
|
|
const double dDenominator = dStartDistance - dEndDistance;
|
|
const double dT = dDenominator == 0.0
|
|
? 0.0 : dStartDistance / dDenominator;
|
|
return QPointF(s.x() + dT * (e.x() - s.x()),
|
|
s.y() + dT * (e.y() - s.y()));
|
|
}
|
|
|
|
QVector<QPointF> clipByEdge(const QVector<QPointF>& vecSubject,
|
|
const QPointF& a,
|
|
const QPointF& b,
|
|
double dAreaTolerance)
|
|
{
|
|
// Sutherland-Hodgman 半平面剪裁，裁剪边左侧视为三角形内部。
|
|
QVector<QPointF> vecOutput;
|
|
if(vecSubject.isEmpty()) {
|
|
return vecOutput;
|
|
}
|
|
QPointF s = vecSubject.last();
|
|
double dStartDistance = crossProduct(a, b, s);
|
|
bool bStartInside = dStartDistance >= -dAreaTolerance;
|
|
for(int nIndex = 0; nIndex < vecSubject.size(); ++nIndex) {
|
|
const QPointF e = vecSubject[nIndex];
|
|
const double dEndDistance = crossProduct(a, b, e);
|
|
const bool bEndInside = dEndDistance >= -dAreaTolerance;
|
|
if(bEndInside) {
|
|
if(!bStartInside) {
|
|
vecOutput.append(segmentLineIntersection(
|
|
s, e, dStartDistance, dEndDistance));
|
|
}
|
|
vecOutput.append(e);
|
|
} else if(bStartInside) {
|
|
vecOutput.append(segmentLineIntersection(
|
|
s, e, dStartDistance, dEndDistance));
|
|
}
|
|
s = e;
|
|
dStartDistance = dEndDistance;
|
|
bStartInside = bEndInside;
|
|
}
|
|
return vecOutput;
|
|
}
|
|
|
|
double triangleIntersectionArea(const Triangle2D& oSubject,
|
|
const Triangle2D& oClip,
|
|
double dAreaTolerance)
|
|
{
|
|
QVector<QPointF> vecClipped;
|
|
vecClipped.append(oSubject.m_a);
|
|
vecClipped.append(oSubject.m_b);
|
|
vecClipped.append(oSubject.m_c);
|
|
vecClipped = clipByEdge(vecClipped,
|
|
oClip.m_a, oClip.m_b, dAreaTolerance);
|
|
vecClipped = clipByEdge(vecClipped,
|
|
oClip.m_b, oClip.m_c, dAreaTolerance);
|
|
vecClipped = clipByEdge(vecClipped,
|
|
oClip.m_c, oClip.m_a, dAreaTolerance);
|
|
return vecClipped.size() >= 3
|
|
? qAbs(signedDoubleArea(vecClipped)) * 0.5 : 0.0;
|
|
}
|
|
|
|
bool rectanglesOverlap(const QRectF& a,
|
|
const QRectF& b,
|
|
double dTolerance)
|
|
{
|
|
return a.right() >= b.left() - dTolerance &&
|
|
b.right() >= a.left() - dTolerance &&
|
|
a.bottom() >= b.top() - dTolerance &&
|
|
b.bottom() >= a.top() - dTolerance;
|
|
}
|
|
|
|
bool polygonsHavePositiveAreaIntersection(
|
|
const PreparedPolygon& a,
|
|
const PreparedPolygon& b)
|
|
{
|
|
// 业务语义是“面积相交”：仅边界或顶点接触不把单元纳入统计。
|
|
const double dFirstScale =
|
|
qMax(a.m_oBounds.width(), a.m_oBounds.height());
|
|
const double dSecondScale =
|
|
qMax(b.m_oBounds.width(), b.m_oBounds.height());
|
|
// 交叠容差跟随较小多边形，避免大选区吞掉小单元的真实窄交叠。
|
|
const double dScale = qMax(1.0, qMin(dFirstScale, dSecondScale));
|
|
const double dLinearTolerance = dScale * 1.0e-10;
|
|
if(!rectanglesOverlap(a.m_oBounds, b.m_oBounds,
|
|
dLinearTolerance)) {
|
|
return false;
|
|
}
|
|
const double dAreaTolerance = dScale * dScale * 1.0e-13;
|
|
const double dPositiveAreaThreshold = qMax(
|
|
dAreaTolerance,
|
|
qMin(a.m_dArea, b.m_dArea) * 1.0e-12);
|
|
for(int nA = 0; nA < a.m_vecTriangles.size(); ++nA) {
|
|
for(int nB = 0; nB < b.m_vecTriangles.size(); ++nB) {
|
|
if(triangleIntersectionArea(
|
|
a.m_vecTriangles[nA], b.m_vecTriangles[nB],
|
|
dAreaTolerance) > dPositiveAreaThreshold) {
|
|
return true;
|
|
}
|
|
}
|
|
}
|
|
return false;
|
|
}
|
|
|
|
bool readCellPolygon(vtkCell* pCell,
|
|
QVector<QPointF>& vecPoints)
|
|
{
|
|
// 当前 PEBI 结果只接受可直接解释为平面多边形的单元类型。
|
|
vecPoints.clear();
|
|
if(pCell == NULL ||
|
|
(pCell->GetCellType() != VTK_TRIANGLE &&
|
|
pCell->GetCellType() != VTK_QUAD &&
|
|
pCell->GetCellType() != VTK_POLYGON) ||
|
|
pCell->GetPoints() == NULL ||
|
|
pCell->GetNumberOfPoints() < 3) {
|
|
return false;
|
|
}
|
|
double aPoint[3] = { 0.0, 0.0, 0.0 };
|
|
for(vtkIdType nPointIndex = 0;
|
|
nPointIndex < pCell->GetNumberOfPoints(); ++nPointIndex) {
|
|
pCell->GetPoints()->GetPoint(nPointIndex, aPoint);
|
|
if(!isFiniteValue(aPoint[0]) || !isFiniteValue(aPoint[1])) {
|
|
vecPoints.clear();
|
|
return false;
|
|
}
|
|
vecPoints.append(QPointF(aPoint[0], aPoint[1]));
|
|
}
|
|
return true;
|
|
}
|
|
|
|
void addCompensated(double dValue,
|
|
double& dSum,
|
|
double& dCorrection)
|
|
{
|
|
// Kahan 补偿求和用于缓解单元数量较多或权重差异较大时的累计误差。
|
|
const double dAdjusted = dValue - dCorrection;
|
|
const double dCandidate = dSum + dAdjusted;
|
|
dCorrection = (dCandidate - dSum) - dAdjusted;
|
|
dSum = dCandidate;
|
|
}
|
|
|
|
}
|
|
|
|
class nmPebiRegionPressureCalculatorPrivate
|
|
{
|
|
public:
|
|
nmPebiRegionPressureCalculatorPrivate()
|
|
: m_bPrepared(false),
|
|
m_dSelectedCellArea(0.0)
|
|
{
|
|
}
|
|
|
|
bool m_bPrepared;
|
|
QSharedPointer<const nmPebiResultSnapshot> m_pSnapshot;
|
|
QVector<QPointF> m_vecPolygonPoints;
|
|
QVector<vtkIdType> m_vecIntersectingCellIds;
|
|
QVector<CachedCell> m_vecCells;
|
|
double m_dSelectedCellArea;
|
|
};
|
|
|
|
nmPebiRegionPressureStatistics::nmPebiRegionPressureStatistics()
|
|
{
|
|
clear();
|
|
}
|
|
|
|
void nmPebiRegionPressureStatistics::clear()
|
|
{
|
|
m_bAvailable = false;
|
|
m_nPressureFrameIndex = -1;
|
|
m_dPressureTimeHours = 0.0;
|
|
m_nIntersectingCellCount = 0;
|
|
m_nValidCellCount = 0;
|
|
m_nIgnoredCellCount = 0;
|
|
m_dSelectedCellAreaSquareMeters = 0.0;
|
|
m_dEffectivePoreVolumeCubicMeters = 0.0;
|
|
m_dMinimumPressureMpa = 0.0;
|
|
m_dMaximumPressureMpa = 0.0;
|
|
m_dPoreVolumeWeightedPressureMpa = 0.0;
|
|
}
|
|
|
|
nmPebiRegionPressureCalculator::nmPebiRegionPressureCalculator()
|
|
: m_pPrivate(new nmPebiRegionPressureCalculatorPrivate)
|
|
{
|
|
}
|
|
|
|
nmPebiRegionPressureCalculator::~nmPebiRegionPressureCalculator()
|
|
{
|
|
delete m_pPrivate;
|
|
m_pPrivate = NULL;
|
|
}
|
|
|
|
bool nmPebiRegionPressureCalculator::prepare(
|
|
const QSharedPointer<const nmPebiResultSnapshot>& pSnapshot,
|
|
const QVector<QPointF>& vecPolygonPoints,
|
|
QString* pError)
|
|
{
|
|
if(pError != NULL) {
|
|
pError->clear();
|
|
}
|
|
const vtkIdType nCellCount = pSnapshot.isNull()
|
|
? 0 : pSnapshot->getResultCellCount();
|
|
if(m_pPrivate == NULL || pSnapshot.isNull() ||
|
|
!pSnapshot->m_bComplete || pSnapshot->m_pResultGrid == NULL ||
|
|
nCellCount <= 0 || pSnapshot->m_pPorosity == NULL ||
|
|
pSnapshot->m_pThickness == NULL ||
|
|
pSnapshot->m_pPorosity->GetNumberOfComponents() != 1 ||
|
|
pSnapshot->m_pThickness->GetNumberOfComponents() != 1 ||
|
|
pSnapshot->m_pPorosity->GetNumberOfTuples() != nCellCount ||
|
|
pSnapshot->m_pThickness->GetNumberOfTuples() != nCellCount ||
|
|
pSnapshot->m_vecPressureFrames.isEmpty()) {
|
|
if(pError != NULL) {
|
|
*pError = QString::fromLatin1(
|
|
"The PEBI result snapshot is invalid.");
|
|
}
|
|
return false;
|
|
}
|
|
|
|
PreparedPolygon oSelection;
|
|
if(!preparePolygon(vecPolygonPoints, oSelection)) {
|
|
if(pError != NULL) {
|
|
*pError = QString::fromLatin1(
|
|
"The region polygon is invalid or self-intersecting.");
|
|
}
|
|
return false;
|
|
}
|
|
|
|
// 先在候选状态中完成全部校验和缓存，失败时不破坏上一个可用区域。
|
|
nmPebiRegionPressureCalculatorPrivate oCandidate;
|
|
oCandidate.m_pSnapshot = pSnapshot;
|
|
oCandidate.m_vecPolygonPoints = oSelection.m_vecPoints;
|
|
|
|
vtkUnstructuredGrid* pGrid = pSnapshot->m_pResultGrid;
|
|
QVector<QPointF> vecCellPoints;
|
|
double dAreaCorrection = 0.0;
|
|
for(vtkIdType nCellId = 0;
|
|
nCellId < pGrid->GetNumberOfCells(); ++nCellId) {
|
|
vtkCell* pCell = pGrid->GetCell(nCellId);
|
|
if(!readCellPolygon(pCell, vecCellPoints)) {
|
|
continue;
|
|
}
|
|
PreparedPolygon oCellPolygon;
|
|
if(!preparePolygon(vecCellPoints, oCellPolygon) ||
|
|
!polygonsHavePositiveAreaIntersection(
|
|
oSelection, oCellPolygon)) {
|
|
continue;
|
|
}
|
|
|
|
CachedCell oCell;
|
|
oCell.m_nCellId = nCellId;
|
|
// 选中规则是“有正面积相交即纳入”，纳入后使用整个单元面积计算孔隙体积权重。
|
|
oCell.m_dArea = oCellPolygon.m_dArea;
|
|
const double dPorosity =
|
|
pSnapshot->m_pPorosity->GetValue(nCellId);
|
|
const double dThickness =
|
|
pSnapshot->m_pThickness->GetValue(nCellId);
|
|
const double dWeight = oCell.m_dArea * dThickness * dPorosity;
|
|
// 几何相交单元保留在计数中；非法物性只在压力统计时被忽略。
|
|
oCell.m_bWeightValid = isFiniteValue(dPorosity) &&
|
|
isFiniteValue(dThickness) && isFiniteValue(dWeight) &&
|
|
dPorosity > 0.0 && dThickness > 0.0 && dWeight > 0.0;
|
|
oCell.m_dPoreVolumeWeight = oCell.m_bWeightValid ? dWeight : 0.0;
|
|
oCandidate.m_vecCells.append(oCell);
|
|
oCandidate.m_vecIntersectingCellIds.append(nCellId);
|
|
addCompensated(oCell.m_dArea,
|
|
oCandidate.m_dSelectedCellArea,
|
|
dAreaCorrection);
|
|
}
|
|
if(!isFiniteValue(oCandidate.m_dSelectedCellArea)) {
|
|
if(pError != NULL) {
|
|
*pError = QString::fromLatin1(
|
|
"The selected cell area is invalid.");
|
|
}
|
|
return false;
|
|
}
|
|
|
|
oCandidate.m_bPrepared = true;
|
|
*m_pPrivate = oCandidate;
|
|
return true;
|
|
}
|
|
|
|
void nmPebiRegionPressureCalculator::clear()
|
|
{
|
|
if(m_pPrivate != NULL) {
|
|
*m_pPrivate = nmPebiRegionPressureCalculatorPrivate();
|
|
}
|
|
}
|
|
|
|
bool nmPebiRegionPressureCalculator::isPrepared() const
|
|
{
|
|
return m_pPrivate != NULL && m_pPrivate->m_bPrepared &&
|
|
!m_pPrivate->m_pSnapshot.isNull();
|
|
}
|
|
|
|
QString nmPebiRegionPressureCalculator::getSnapshotId() const
|
|
{
|
|
return isPrepared()
|
|
? m_pPrivate->m_pSnapshot->getSnapshotId() : QString();
|
|
}
|
|
|
|
const QVector<QPointF>&
|
|
nmPebiRegionPressureCalculator::getPolygonPoints() const
|
|
{
|
|
return m_pPrivate->m_vecPolygonPoints;
|
|
}
|
|
|
|
const QVector<vtkIdType>&
|
|
nmPebiRegionPressureCalculator::getIntersectingCellIds() const
|
|
{
|
|
return m_pPrivate->m_vecIntersectingCellIds;
|
|
}
|
|
|
|
bool nmPebiRegionPressureCalculator::calculateFrame(
|
|
int nPressureFrameIndex,
|
|
nmPebiRegionPressureStatistics& oStatistics,
|
|
QString* pError) const
|
|
{
|
|
oStatistics.clear();
|
|
if(pError != NULL) {
|
|
pError->clear();
|
|
}
|
|
if(!isPrepared() || nPressureFrameIndex < 0 ||
|
|
nPressureFrameIndex >=
|
|
m_pPrivate->m_pSnapshot->m_vecPressureFrames.size()) {
|
|
if(pError != NULL) {
|
|
*pError = QString::fromLatin1(
|
|
"The pressure frame is invalid.");
|
|
}
|
|
return false;
|
|
}
|
|
|
|
const nmPebiResultSnapshot::PressureFrame& oFrame =
|
|
m_pPrivate->m_pSnapshot->m_vecPressureFrames[
|
|
nPressureFrameIndex];
|
|
vtkDoubleArray* pPressure = oFrame.m_pPressure;
|
|
if(pPressure == NULL || pPressure->GetNumberOfComponents() != 1 ||
|
|
pPressure->GetNumberOfTuples() !=
|
|
m_pPrivate->m_pSnapshot->getResultCellCount() ||
|
|
!isFiniteValue(oFrame.m_dTime)) {
|
|
if(pError != NULL) {
|
|
*pError = QString::fromLatin1(
|
|
"The pressure frame data is invalid.");
|
|
}
|
|
return false;
|
|
}
|
|
|
|
oStatistics.m_nPressureFrameIndex = nPressureFrameIndex;
|
|
oStatistics.m_dPressureTimeHours = oFrame.m_dTime;
|
|
oStatistics.m_nIntersectingCellCount =
|
|
m_pPrivate->m_vecCells.size();
|
|
oStatistics.m_dSelectedCellAreaSquareMeters =
|
|
m_pPrivate->m_dSelectedCellArea;
|
|
|
|
// 交集及静态权重已在 prepare() 缓存，逐帧只读压力并做线性汇总。
|
|
double dPoreVolumeSum = 0.0;
|
|
double dPoreVolumeCorrection = 0.0;
|
|
double dWeightedPressureSum = 0.0;
|
|
double dWeightedPressureCorrection = 0.0;
|
|
double dMinimumPressure = DBL_MAX;
|
|
double dMaximumPressure = -DBL_MAX;
|
|
for(int nIndex = 0;
|
|
nIndex < m_pPrivate->m_vecCells.size(); ++nIndex) {
|
|
const CachedCell& oCell = m_pPrivate->m_vecCells[nIndex];
|
|
const double dPressure = pPressure->GetValue(oCell.m_nCellId);
|
|
if(!oCell.m_bWeightValid || !isFiniteValue(dPressure)) {
|
|
continue;
|
|
}
|
|
addCompensated(oCell.m_dPoreVolumeWeight,
|
|
dPoreVolumeSum, dPoreVolumeCorrection);
|
|
addCompensated(oCell.m_dPoreVolumeWeight * dPressure,
|
|
dWeightedPressureSum,
|
|
dWeightedPressureCorrection);
|
|
dMinimumPressure = qMin(dMinimumPressure, dPressure);
|
|
dMaximumPressure = qMax(dMaximumPressure, dPressure);
|
|
++oStatistics.m_nValidCellCount;
|
|
}
|
|
oStatistics.m_nIgnoredCellCount =
|
|
oStatistics.m_nIntersectingCellCount -
|
|
oStatistics.m_nValidCellCount;
|
|
// 空选区或全部权重无效都是可展示结果，不把它们当成计算错误。
|
|
if(oStatistics.m_nValidCellCount == 0) {
|
|
return true;
|
|
}
|
|
if(!isFiniteValue(dPoreVolumeSum) || dPoreVolumeSum <= 0.0 ||
|
|
!isFiniteValue(dWeightedPressureSum)) {
|
|
if(pError != NULL) {
|
|
*pError = QString::fromLatin1(
|
|
"The pore-volume weighted pressure is invalid.");
|
|
}
|
|
return false;
|
|
}
|
|
|
|
const double dAveragePressure =
|
|
dWeightedPressureSum / dPoreVolumeSum;
|
|
if(!isFiniteValue(dMinimumPressure) ||
|
|
!isFiniteValue(dMaximumPressure) ||
|
|
!isFiniteValue(dAveragePressure)) {
|
|
if(pError != NULL) {
|
|
*pError = QString::fromLatin1(
|
|
"The regional pressure statistics are invalid.");
|
|
}
|
|
return false;
|
|
}
|
|
oStatistics.m_bAvailable = true;
|
|
oStatistics.m_dEffectivePoreVolumeCubicMeters = dPoreVolumeSum;
|
|
oStatistics.m_dMinimumPressureMpa = dMinimumPressure;
|
|
oStatistics.m_dMaximumPressureMpa = dMaximumPressure;
|
|
oStatistics.m_dPoreVolumeWeightedPressureMpa = dAveragePressure;
|
|
return true;
|
|
}
|