nmWTAI-Platform/Src/nmNum/nmCalculation/nmPebiStaticPropertyCalculator.cpp

#include "nmPebiStaticPropertyCalculator.h"

#include "nmCalculationUtils.h"

#include <QtGlobal>

#include <float.h>
#include <math.h>

#if !defined(_MSC_VER)
#include <cmath>
#endif

namespace {

enum StaticPropertyKind
{
    StaticProperty_Permeability = 0,
    StaticProperty_Porosity,
    StaticProperty_Thickness,
    StaticProperty_Invalid
};

bool isCancellationRequested(const QAtomicInt* pCancelRequested)
{
    return pCancelRequested != NULL &&
            static_cast<int>(*pCancelRequested) != 0;
}

bool isFiniteValue(double dValue)
{
#if defined(_MSC_VER)
    return _finite(dValue) != 0;
#else
    return std::isfinite(dValue);
#endif
}

StaticPropertyKind propertyKindForLegacyKey(const QString& sPropertyKey)
{
    // 旧工程持久化键必须保持封闭；正式属性目录中的 pressure 不能进入此路径。
    if(sPropertyKey == QLatin1String("k")) {
        return StaticProperty_Permeability;
    }
    if(sPropertyKey == QLatin1String("phi")) {
        return StaticProperty_Porosity;
    }
    if(sPropertyKey == QLatin1String("h")) {
        return StaticProperty_Thickness;
    }
    return StaticProperty_Invalid;
}

nmPebiStaticPropertyCalculationStatus cancelled(QString* pErrorMessage)
{
    // 主动停止不是插值错误，交由任务状态单独报告。
    if(pErrorMessage != NULL) {
        pErrorMessage->clear();
    }
    return NM_PebiStaticPropertyCalculation_Cancelled;
}

nmPebiStaticPropertyCalculationStatus failed(
        QString* pErrorMessage,
        const QString& sErrorMessage)
{
    if(pErrorMessage != NULL) {
        *pErrorMessage = sErrorMessage;
    }
    return NM_PebiStaticPropertyCalculation_Failed;
}

QVector<double>* valuesForKind(
        nmPebiStaticPropertyValues& oValues,
        StaticPropertyKind eKind)
{
    switch(eKind) {
    case StaticProperty_Permeability:
        return &oValues.m_vecPermeabilityMilliDarcy;
    case StaticProperty_Porosity:
        return &oValues.m_vecPorosity;
    case StaticProperty_Thickness:
        return &oValues.m_vecThicknessMeters;
    default:
        return NULL;
    }
}

}

nmPebiStaticPropertyCalculationInput::
nmPebiStaticPropertyCalculationInput()
    : m_bCaptured(false),
      m_dReferencePermeabilityMilliDarcy(0.0),
      m_dReferencePorosity(0.0),
      m_dReferenceThicknessMeters(0.0)
{
}

nmPebiStaticPropertyCalculationStatus
nmPebiStaticPropertyCalculator::calculate(
        const QVector<QPointF>& vecOriginalCellPropertyLocations,
        const nmPebiStaticPropertyCalculationInput& oInput,
        nmPebiStaticPropertyValues& oOutput,
        QString* pErrorMessage,
        const QAtomicInt* pCancelRequested)
{
    oOutput.clear();
    if(pErrorMessage != NULL) {
        pErrorMessage->clear();
    }

    if(isCancellationRequested(pCancelRequested)) {
        return cancelled(pErrorMessage);
    }
    if(!oInput.m_bCaptured) {
        return failed(pErrorMessage,
                QLatin1String("Static property input was not captured."));
    }
    if(vecOriginalCellPropertyLocations.isEmpty()) {
        return failed(pErrorMessage,
                QLatin1String("No PEBI property locations are available."));
    }
    if(!isFiniteValue(oInput.m_dReferencePermeabilityMilliDarcy) ||
       !isFiniteValue(oInput.m_dReferencePorosity) ||
       !isFiniteValue(oInput.m_dReferenceThicknessMeters)) {
        return failed(pErrorMessage,
                QLatin1String("A reference static property value is invalid."));
    }

    for(int nIndex = 0;
        nIndex < vecOriginalCellPropertyLocations.size(); ++nIndex) {
        if((nIndex % 256) == 0 &&
           isCancellationRequested(pCancelRequested)) {
            return cancelled(pErrorMessage);
        }
        const QPointF& oLocation =
                vecOriginalCellPropertyLocations[nIndex];
        if(!isFiniteValue(oLocation.x()) || !isFiniteValue(oLocation.y())) {
            return failed(pErrorMessage,
                    QString::fromLatin1(
                        "PEBI property location %1 contains an invalid coordinate.")
                    .arg(nIndex));
        }
    }

    // 先完整检查启用键及唯一性，避免完成一次昂贵 Kriging 后才发现重复数据组。
    bool aPropertyEnabled[3] = { false, false, false };
    for(int nDataSetIndex = 0;
        nDataSetIndex < oInput.m_vecDataSets.size(); ++nDataSetIndex) {
        if(isCancellationRequested(pCancelRequested)) {
            return cancelled(pErrorMessage);
        }
        const nmPropertyInterpolationDataSet& oDataSet =
                oInput.m_vecDataSets[nDataSetIndex];
        if(!oDataSet.useForCalculation) {
            continue;
        }

        const StaticPropertyKind eKind =
                propertyKindForLegacyKey(oDataSet.property);
        if(eKind == StaticProperty_Invalid) {
            return failed(pErrorMessage,
                    QString::fromLatin1(
                        "Dataset '%1' has an unknown static property key.")
                    .arg(oDataSet.name));
        }
        if(aPropertyEnabled[static_cast<int>(eKind)]) {
            return failed(pErrorMessage,
                    QString::fromLatin1(
                        "More than one dataset is enabled for property %1.")
                    .arg(oDataSet.property));
        }
        aPropertyEnabled[static_cast<int>(eKind)] = true;
    }

    const int nCellCount = vecOriginalCellPropertyLocations.size();
    nmPebiStaticPropertyValues oCandidate;
    oCandidate.m_vecPermeabilityMilliDarcy.fill(
                oInput.m_dReferencePermeabilityMilliDarcy, nCellCount);
    oCandidate.m_vecPorosity.fill(
                oInput.m_dReferencePorosity, nCellCount);
    oCandidate.m_vecThicknessMeters.fill(
                oInput.m_dReferenceThicknessMeters, nCellCount);

    for(int nDataSetIndex = 0;
        nDataSetIndex < oInput.m_vecDataSets.size(); ++nDataSetIndex) {
        if(isCancellationRequested(pCancelRequested)) {
            return cancelled(pErrorMessage);
        }
        const nmPropertyInterpolationDataSet& oDataSet =
                oInput.m_vecDataSets[nDataSetIndex];
        if(!oDataSet.useForCalculation) {
            continue;
        }

        QVector<QPointF> vecMeasurementPoints;
        QVector<double> vecMeasurementValues;
        vecMeasurementPoints.reserve(oDataSet.points.size());
        vecMeasurementValues.reserve(oDataSet.points.size());
        for(int nPointIndex = 0;
            nPointIndex < oDataSet.points.size(); ++nPointIndex) {
            const nmPropertyInterpolationPointData& oPoint =
                    oDataSet.points[nPointIndex];
            vecMeasurementPoints.append(QPointF(oPoint.x, oPoint.y));
            vecMeasurementValues.append(oPoint.value);
        }

        QVector<double> vecInterpolationValues;
        QString sCalculationError;
        if(!nmCalculationUtils::calculateKriging(
               vecOriginalCellPropertyLocations,
               vecMeasurementPoints,
               vecMeasurementValues,
               oDataSet.nugget,
               oDataSet.sill,
               oDataSet.range,
               oDataSet.model,
               oInput.m_sLicensePath,
               vecInterpolationValues,
               &sCalculationError,
               pCancelRequested)) {
            // calculateKriging() 只有 bool 返回值，取消必须读取原子状态判定。
            if(isCancellationRequested(pCancelRequested)) {
                return cancelled(pErrorMessage);
            }
            return failed(pErrorMessage,
                    QString::fromLatin1("Dataset '%1': %2")
                    .arg(oDataSet.name)
                    .arg(sCalculationError));
        }

        const StaticPropertyKind eKind =
                propertyKindForLegacyKey(oDataSet.property);
        QVector<double>* pTargetValues = valuesForKind(oCandidate, eKind);
        if(pTargetValues == NULL ||
           vecInterpolationValues.size() != nCellCount) {
            return failed(pErrorMessage,
                    QString::fromLatin1(
                        "Dataset '%1' returned an invalid value count.")
                    .arg(oDataSet.name));
        }
        *pTargetValues = vecInterpolationValues;
    }

    if(!oCandidate.isValidForCellCount(nCellCount)) {
        return failed(pErrorMessage,
                QLatin1String("The calculated static properties are invalid."));
    }
    if(isCancellationRequested(pCancelRequested)) {
        return cancelled(pErrorMessage);
    }

    // 三个 qSwap 均为容器所有权交换；到达此处后不再执行可能失败的操作。
    qSwap(oOutput.m_vecPermeabilityMilliDarcy,
          oCandidate.m_vecPermeabilityMilliDarcy);
    qSwap(oOutput.m_vecPorosity, oCandidate.m_vecPorosity);
    qSwap(oOutput.m_vecThicknessMeters, oCandidate.m_vecThicknessMeters);
    return NM_PebiStaticPropertyCalculation_Succeeded;
}