diff --git a/Src/nmNum/nmCalculation/nmCalculationDllPebiSolverTask.cpp b/Src/nmNum/nmCalculation/nmCalculationDllPebiSolverTask.cpp index e8db3a5..7951652 100644 --- a/Src/nmNum/nmCalculation/nmCalculationDllPebiSolverTask.cpp +++ b/Src/nmNum/nmCalculation/nmCalculationDllPebiSolverTask.cpp @@ -16,6 +16,7 @@ #include "nmDataTimeStepSetting.h" #include "nmCalculationPebiGrid.h" +#include "nmCalculationUtils.h" #include "nmDataAnalyzeManager.h" #include "nmDataPvtParaForPebi.h" @@ -52,6 +53,111 @@ bool isFiniteSolverNumber(double value) #endif } +// 将启用的数据组插值到全部网格单元中心,并覆盖对应的求解器属性数组. +bool applyPropertyInterpolation( + HX_NWTM_MODEL_INPUT& modelInput, + const QVector& dataSets, + const QString& licensePath, + QString& errorMessage) +{ + bool hasEnabledDataSet = false; + for(int i = 0; i < dataSets.size(); ++i) { + if(dataSets[i].useForCalculation) { + hasEnabledDataSet = true; + break; + } + } + + if(!hasEnabledDataSet) { + return true; + } + + QVector targetPoints; + targetPoints.reserve(static_cast(modelInput.GRID.Trinodexy.size())); + for(size_t cellIndex = 0; + cellIndex < modelInput.GRID.Trinodexy.size(); ++cellIndex) { + const dVec1& cellPosition = modelInput.GRID.Trinodexy[cellIndex]; + if(cellPosition.size() < 2) { + errorMessage = QString("Grid cell %1 has no valid center coordinate.") + .arg(static_cast(cellIndex)); + return false; + } + targetPoints.append(QPointF(cellPosition[0], cellPosition[1])); + } + + bool kApplied = false; + bool phiApplied = false; + bool hApplied = false; + for(int dataSetIndex = 0; dataSetIndex < dataSets.size(); ++dataSetIndex) { + const nmPropertyInterpolationDataSet& dataSet = dataSets[dataSetIndex]; + if(!dataSet.useForCalculation) { + continue; + } + + dVec1* solverValues = NULL; + bool* propertyApplied = NULL; + if(dataSet.property == "k") { + solverValues = &modelInput.Base.k; + propertyApplied = &kApplied; + } else if(dataSet.property == "phi") { + solverValues = &modelInput.Base.phi; + propertyApplied = &phiApplied; + } else if(dataSet.property == "h") { + solverValues = &modelInput.Base.h; + propertyApplied = &hApplied; + } else { + errorMessage = QString("Dataset '%1' has an unknown property.") + .arg(dataSet.name); + return false; + } + + if(*propertyApplied) { + errorMessage = QString( + "More than one dataset is enabled for property %1.") + .arg(dataSet.property); + return false; + } + + QVector measurementPoints; + QVector measurementValues; + measurementPoints.reserve(dataSet.points.size()); + measurementValues.reserve(dataSet.points.size()); + for(int pointIndex = 0; pointIndex < dataSet.points.size(); ++pointIndex) { + const nmPropertyInterpolationPointData& point = dataSet.points[pointIndex]; + measurementPoints.append(QPointF(point.x, point.y)); + measurementValues.append(point.value); + } + + QVector interpolationValues; + QString calculationError; + if(!nmCalculationUtils::calculateKriging( + targetPoints, + measurementPoints, + measurementValues, + dataSet.nugget, + dataSet.sill, + dataSet.range, + dataSet.model, + licensePath, + interpolationValues, + &calculationError)) { + errorMessage = QString("Dataset '%1': %2") + .arg(dataSet.name) + .arg(calculationError); + return false; + } + + solverValues->resize(interpolationValues.size()); + for(int valueIndex = 0; valueIndex < interpolationValues.size(); ++valueIndex) { + // 插值数据已经使用求解器基准单位,结果可直接按网格顺序写入. + (*solverValues)[valueIndex] = interpolationValues[valueIndex]; + } + *propertyApplied = true; + } + + return true; +} + bool isReasonableLogLogValue(double value) { const double maxReasonableAbsValue = 1.0e12; @@ -485,6 +591,21 @@ bool nmCalculationDllPebiSolverTask::execPebiMode() p0.Base.h = dVec1(cellCount, pReservoirData->getThickness().getValue().toDouble()); } + // 未启用的属性保留上面的单值数组;启用后才使用所选数据组覆盖. + QString interpolationError; + if(!applyPropertyInterpolation( + p0, + pDataInstance->getPropertyInterpolationDataSets(), + pDataInstance->getLicensePath(), + interpolationError)) { + QString logMessage = QString("Property interpolation failed: %1") + .arg(interpolationError); + qWarning() << logMessage; + zxLogInstance::getInstance()->writeLogF(logMessage); + FreeLibrary(dll); + return false; + } + // 获取时间步长 nmDataTimeStepSetting* pTimeStepSetting = pDataInstance->getTimeStep(); if(pTimeStepSetting) {