接入网格属性插值求解流程

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<nmPropertyInterpolationDataSet>& 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<QPointF> targetPoints;
+    targetPoints.reserve(static_cast<int>(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<qulonglong>(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<QPointF> measurementPoints;
+        QVector<double> 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<double> 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) {