diff --git a/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.cpp b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.cpp index d46f9f9b60..0f6bdead28 100644 --- a/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.cpp +++ b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.cpp @@ -53,8 +53,31 @@ #include #include +#include "RigHexIntersectionTools.h" +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +const std::string RigFemPartResultsCollection::FIELD_NAME_COMPACTION = "COMPACTION"; + +//-------------------------------------------------------------------------------------------------- +/// Internal definitions +//-------------------------------------------------------------------------------------------------- +class RefElement +{ +public: + size_t elementIdx; + float intersectionPointToCurrentNodeDistance; + cvf::Vec3f intersectionPoint; + std::vector elementFaceNodeIdxs; +}; + +static std::vector coordsFromNodeIndices(const RigFemPart& part, const std::vector& nodeIdxs); +static std::vector nodesForElement(const RigFemPart& part, size_t elementIdx); +static float horizontalDistance(const cvf::Vec3f& p1, const cvf::Vec3f& p2); +static void findReferenceElementForNode(const RigFemPart& part, size_t nodeIdx, size_t kRefLayer, RefElement *refElement); + //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- @@ -295,6 +318,7 @@ std::map > RigFemPartResultsCollection::sc { fieldCompNames = m_readerInterface->scalarNodeFieldAndComponentNames(); fieldCompNames["POR-Bar"]; + fieldCompNames[FIELD_NAME_COMPACTION]; } else if (resPos == RIG_ELEMENT_NODAL) { @@ -559,7 +583,7 @@ RigFemScalarResultFrames* RigFemPartResultsCollection::calculateTimeLapseResult( frameCountProgress.setProgressDescription("Calculating " + QString::fromStdString(resVarAddr.fieldName + ": " + resVarAddr.componentName)); frameCountProgress.setNextProgressIncrement(this->frameCount()); - RigFemResultAddress resVarNative(resVarAddr.resultPosType, resVarAddr.fieldName, resVarAddr.componentName); + RigFemResultAddress resVarNative(resVarAddr.resultPosType, resVarAddr.fieldName, resVarAddr.componentName, RigFemResultAddress::NO_TIME_LAPSE, resVarAddr.refKLayerIndex); RigFemScalarResultFrames * srcDataFrames = this->findOrLoadScalarResult(partIndex, resVarNative); RigFemScalarResultFrames * dstDataFrames = m_femPartResults[partIndex]->createScalarResult(resVarAddr); @@ -1427,6 +1451,63 @@ RigFemScalarResultFrames* RigFemPartResultsCollection::calculatePrincipalStrainV return requestedPrincipal; } +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RigFemScalarResultFrames* RigFemPartResultsCollection::calculateCompactionValues(int partIndex, const RigFemResultAddress &resVarAddr) +{ + CVF_ASSERT(resVarAddr.fieldName == FIELD_NAME_COMPACTION); + + caf::ProgressInfo frameCountProgress(this->frameCount() + 1, ""); + frameCountProgress.setProgressDescription("Calculating " + QString::fromStdString(resVarAddr.fieldName)); + + RigFemScalarResultFrames * u3Frames = this->findOrLoadScalarResult(partIndex, RigFemResultAddress(resVarAddr.resultPosType, "U", "U3")); + frameCountProgress.incrementProgress(); + + RigFemScalarResultFrames* compactionFrames = m_femPartResults[partIndex]->createScalarResult(resVarAddr); + + const RigFemPart* part = m_femParts->part(partIndex); + for (int t = 0; t < u3Frames->frameCount(); t++) + { + std::vector& compactionFrame = compactionFrames->frameData(t); + size_t nodeCount = part->nodes().nodeIds.size(); + + frameCountProgress.incrementProgress(); + + compactionFrame.resize(nodeCount); + for (int n = 0; n < nodeCount; n++) + { + RefElement refElement; + findReferenceElementForNode(*part, n, resVarAddr.refKLayerIndex, &refElement); + + if (refElement.elementIdx != cvf::UNDEFINED_SIZE_T) + { + float shortestDist = HUGE_VALF; + size_t closestNodeIdx = cvf::UNDEFINED_SIZE_T; + + for (size_t nodeIdx : refElement.elementFaceNodeIdxs) + { + float dist = horizontalDistance(refElement.intersectionPoint, part->nodes().coordinates[nodeIdx]); + if (dist < shortestDist) + { + shortestDist = dist; + closestNodeIdx = nodeIdx; + } + } + + compactionFrame[n] = u3Frames->frameData(t)[n] - u3Frames->frameData(t)[closestNodeIdx]; + } + else + { + compactionFrame[n] = HUGE_VAL; + } + } + } + + RigFemScalarResultFrames* requestedPrincipal = this->findOrLoadScalarResult(partIndex, resVarAddr); + return requestedPrincipal; +} + //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- @@ -1453,6 +1534,11 @@ RigFemScalarResultFrames* RigFemPartResultsCollection::calculateDerivedResult(in } } + if (resVarAddr.fieldName == FIELD_NAME_COMPACTION) + { + return calculateCompactionValues(partIndex, resVarAddr); + } + if (resVarAddr.fieldName == "SE" && resVarAddr.componentName == "SFI") { return calculateSFI(partIndex, resVarAddr); @@ -2372,4 +2458,95 @@ RigFemClosestResultIndexCalculator::RigFemClosestResultIndexCalculator(RigFemPar m_resultIndexToClosestResult = elmNodFaceResIdx; m_closestNodeId = femPart->nodes().nodeIds[closestNodeIdx]; } -} \ No newline at end of file +} + +//-------------------------------------------------------------------------------------------------- +/// Internal functions +//-------------------------------------------------------------------------------------------------- + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::vector coordsFromNodeIndices(const RigFemPart& part, const std::vector& nodeIdxs) +{ + std::vector out; + for (const auto& nodeIdx : nodeIdxs) out.push_back(cvf::Vec3d(part.nodes().coordinates[nodeIdx])); + return out; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::vector nodesForElement(const RigFemPart& part, size_t elementIdx) +{ + std::vector nodeIdxs; + const int* nodeConn = part.connectivities(elementIdx); + for (int n = 0; n < 8; n++) nodeIdxs.push_back(nodeConn[n]); + return nodeIdxs; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +float horizontalDistance(const cvf::Vec3f& p1, const cvf::Vec3f& p2) +{ + cvf::Vec3f p1_ = p1; + cvf::Vec3f p2_ = p2; + p1_.z() = p2_.z() = 0; + return p1_.pointDistance(p2_); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void findReferenceElementForNode(const RigFemPart& part, size_t nodeIdx, size_t kRefLayer, RefElement *refElement) +{ + static const double zMin = -1e6, zMax = 1e6; + + cvf::BoundingBox bb; + cvf::Vec3f currentNodeCoord = part.nodes().coordinates[nodeIdx]; + cvf::Vec3f p1 = cvf::Vec3f(currentNodeCoord.x(), currentNodeCoord.y(), zMin); + cvf::Vec3f p2 = cvf::Vec3f(currentNodeCoord.x(), currentNodeCoord.y(), zMax); + bb.add(p1); + bb.add(p2); + + std::vector refElementCandidates; + part.findIntersectingCells(bb, &refElementCandidates); + + const RigFemPartGrid* grid = part.structGrid(); + const std::vector& nodeCoords = part.nodes().coordinates; + + refElement->elementIdx = cvf::UNDEFINED_SIZE_T; + refElement->intersectionPointToCurrentNodeDistance = HUGE_VALF; + size_t i, j, k; + for (const size_t elemIdx : refElementCandidates) + { + grid->ijkFromCellIndex(elemIdx, &i, &j, &k); + if (k == kRefLayer) + { + const std::vector nodeIndices = nodesForElement(part, elemIdx); + CVF_ASSERT(nodeIndices.size() == 8); + + std::vector intersections; + RigHexIntersectionTools::lineHexCellIntersection(cvf::Vec3d(p1), cvf::Vec3d(p2), coordsFromNodeIndices(part, nodeIndices).data(), elemIdx, &intersections); + + for (const auto& intersection : intersections) + { + cvf::Vec3f intersectionPoint = cvf::Vec3f(intersection.m_intersectionPoint); + + float nodeToIntersectionDistance = currentNodeCoord.pointDistance(intersectionPoint); + if (nodeToIntersectionDistance < refElement->intersectionPointToCurrentNodeDistance) + { + cvf::ubyte faceNodes[4]; + grid->cellFaceVertexIndices(intersection.m_face, faceNodes); + std::vector topFaceCoords({ nodeIndices[faceNodes[0]], nodeIndices[faceNodes[1]], nodeIndices[faceNodes[2]], nodeIndices[faceNodes[3]] }); + + refElement->elementIdx = elemIdx; + refElement->intersectionPointToCurrentNodeDistance = nodeToIntersectionDistance; + refElement->intersectionPoint = intersectionPoint; + refElement->elementFaceNodeIdxs = topFaceCoords; + } + } + } + } +} diff --git a/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.h b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.h index 2cf9f147b8..b494d177be 100644 --- a/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.h +++ b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.h @@ -48,6 +48,8 @@ namespace caf class RigFemPartResultsCollection: public cvf::Object { public: + static const std::string FIELD_NAME_COMPACTION; + RigFemPartResultsCollection(RifGeoMechReaderInterface* readerInterface, RifElementPropertyReader* elementPropertyReader, const RigFemPartCollection * femPartCollection); ~RigFemPartResultsCollection(); @@ -115,6 +117,7 @@ private: RigFemScalarResultFrames* calculateSurfaceAngles(int partIndex, const RigFemResultAddress& resVarAddr); RigFemScalarResultFrames* calculatePrincipalStressValues(int partIndex, const RigFemResultAddress &resVarAddr); RigFemScalarResultFrames* calculatePrincipalStrainValues(int partIndex, const RigFemResultAddress &resVarAddr); + RigFemScalarResultFrames* calculateCompactionValues(int partIndex, const RigFemResultAddress &resVarAddr); cvf::Collection m_femPartResults; cvf::ref m_readerInterface; diff --git a/ApplicationCode/GeoMech/GeoMechDataModel/RigFemResultAddress.h b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemResultAddress.h index 942b6e1fbb..2e9f4d92f9 100644 --- a/ApplicationCode/GeoMech/GeoMechDataModel/RigFemResultAddress.h +++ b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemResultAddress.h @@ -36,24 +36,16 @@ public: componentName = ""; } - RigFemResultAddress(RigFemResultPosEnum resPosType, - const std::string& aFieldName, - const std::string& aComponentName) - : resultPosType(resPosType), - fieldName(aFieldName), - componentName(aComponentName), - timeLapseBaseFrameIdx(-1) - { - } - RigFemResultAddress(RigFemResultPosEnum resPosType, const std::string& aFieldName, const std::string& aComponentName, - int aTimeLapseBaseFrame) + int timeLapseBaseFrameIdx = NO_TIME_LAPSE, + int refKLayerIndex = NO_COMPACTION) : resultPosType(resPosType), fieldName(aFieldName), componentName(aComponentName), - timeLapseBaseFrameIdx(aTimeLapseBaseFrame) + timeLapseBaseFrameIdx(timeLapseBaseFrameIdx), + refKLayerIndex(refKLayerIndex) { } @@ -61,12 +53,17 @@ public: std::string fieldName; std::string componentName; int timeLapseBaseFrameIdx; + int refKLayerIndex; static const int ALL_TIME_LAPSES = -2; + static const int NO_TIME_LAPSE = -1; + static const int NO_COMPACTION = -1; - bool isTimeLapse() const { return timeLapseBaseFrameIdx >= 0;} + bool isTimeLapse() const { return timeLapseBaseFrameIdx > NO_TIME_LAPSE;} bool representsAllTimeLapses() const { return timeLapseBaseFrameIdx == ALL_TIME_LAPSES;} + bool isCompaction() const { return refKLayerIndex > NO_COMPACTION; } + bool isValid() const { bool isTypeValid = resultPosType == RIG_NODAL @@ -97,6 +94,11 @@ public: return (fieldName < other.fieldName); } + if (refKLayerIndex != other.refKLayerIndex) + { + return refKLayerIndex < other.refKLayerIndex; + } + return (componentName < other.componentName); } diff --git a/ApplicationCode/ProjectDataModel/RimGeoMechResultDefinition.cpp b/ApplicationCode/ProjectDataModel/RimGeoMechResultDefinition.cpp index d799a8846a..7e306b7341 100644 --- a/ApplicationCode/ProjectDataModel/RimGeoMechResultDefinition.cpp +++ b/ApplicationCode/ProjectDataModel/RimGeoMechResultDefinition.cpp @@ -24,6 +24,8 @@ #include "RigFemPartResultsCollection.h" #include "RigFemResultAddress.h" #include "RigGeoMechCaseData.h" +#include "RigFemPartCollection.h" +#include "RigFemPartGrid.h" #include "RiaDefines.h" #include "RimGeoMechCase.h" @@ -77,6 +79,12 @@ RimGeoMechResultDefinition::RimGeoMechResultDefinition(void) CAF_PDM_InitField(&m_timeLapseBaseTimestep, "TimeLapseBaseTimeStep", 0, "Base Time Step", "", "", ""); m_timeLapseBaseTimestep.uiCapability()->setUiHidden(true); + CAF_PDM_InitField(&m_isCompaction, "IsCompaction", false, "Compaction", "", "", ""); + m_isCompaction.uiCapability()->setUiHidden(true); + + CAF_PDM_InitField(&m_compactionRefLayer, "CompactionRefLayer", 0, "Compaction Ref Layer", "", "", ""); + m_compactionRefLayer.uiCapability()->setUiHidden(true); + CAF_PDM_InitFieldNoDefault(&m_resultPositionTypeUiField, "ResultPositionTypeUi", "Result Position", "", "", ""); m_resultPositionTypeUiField.xmlCapability()->setIOWritable(false); m_resultPositionTypeUiField.xmlCapability()->setIOReadable(false); @@ -96,6 +104,14 @@ RimGeoMechResultDefinition::RimGeoMechResultDefinition(void) m_resultVariableUiField.uiCapability()->setUiEditorTypeName(caf::PdmUiListEditor::uiEditorTypeName()); m_resultVariableUiField.uiCapability()->setUiLabelPosition(caf::PdmUiItemInfo::TOP); + CAF_PDM_InitField(&m_isCompactionUiField, "IsCompactionUi", false, "Enable Compaction", "", "", ""); + m_isCompactionUiField.xmlCapability()->setIOWritable(false); + m_isCompactionUiField.xmlCapability()->setIOReadable(false); + + CAF_PDM_InitField(&m_compactionRefLayerUiField, "CompactionRefLayerUi", 0, "Compaction Ref Layer", "", "", ""); + m_compactionRefLayerUiField.xmlCapability()->setIOWritable(false); + m_compactionRefLayerUiField.xmlCapability()->setIOReadable(false); + m_isChangedByField = false; } @@ -123,6 +139,14 @@ void RimGeoMechResultDefinition::defineUiOrdering(QString uiConfigName, caf::Pdm timeLapseGr->add(&m_timeLapseBaseTimestepUiField); } + if (m_resultPositionTypeUiField() == RIG_NODAL) + { + caf::PdmUiGroup * compactionGroup = uiOrdering.addNewGroup("Compaction Options"); + compactionGroup->add(&m_isCompactionUiField); + if (m_isCompactionUiField()) + compactionGroup->add(&m_compactionRefLayerUiField); + } + if (!m_isChangedByField) { m_resultPositionTypeUiField = m_resultPositionType; @@ -156,6 +180,8 @@ QList RimGeoMechResultDefinition::calculateValueOptions( for (int oIdx = 0; oIdx < uiVarNames.size(); ++oIdx) { + if (uiVarNames[oIdx].startsWith(RigFemPartResultsCollection::FIELD_NAME_COMPACTION.data()) && !m_isCompactionUiField()) continue; + options.push_back(caf::PdmOptionItemInfo(uiVarNames[oIdx], varNames[oIdx])); } } @@ -177,6 +203,14 @@ QList RimGeoMechResultDefinition::calculateValueOptions( options.push_back(caf::PdmOptionItemInfo(QString::fromStdString(stepNames[stepIdx]), static_cast(stepIdx))); } } + else if (&m_compactionRefLayerUiField == fieldNeedingOptions) + { + size_t kCount = m_geomCase->geoMechData()->femParts()->part(0)->structGrid()->gridPointCountK() - 1; + for (size_t layerIdx = 0; layerIdx < kCount; ++layerIdx) + { + options.push_back(caf::PdmOptionItemInfo(QString::number(layerIdx + 1), (int)layerIdx)); + } + } } return options; @@ -229,7 +263,7 @@ void RimGeoMechResultDefinition::fieldChangedByUi(const caf::PdmFieldHandle* cha RimPlotCurve* curve = nullptr; this->firstAncestorOrThisOfType(curve); - if (&m_resultVariableUiField == changedField) + if (&m_resultVariableUiField == changedField || &m_compactionRefLayerUiField == changedField) { QStringList fieldComponentNames = m_resultVariableUiField().split(QRegExp("\\s+")); if (fieldComponentNames.size() > 0) @@ -255,8 +289,10 @@ void RimGeoMechResultDefinition::fieldChangedByUi(const caf::PdmFieldHandle* cha m_resultComponentName = ""; } - m_isTimeLapseResult = m_isTimeLapseResultUiField(); + m_isTimeLapseResult = m_isTimeLapseResultUiField(); m_timeLapseBaseTimestep = m_timeLapseBaseTimestepUiField(); + m_isCompaction = m_isCompactionUiField(); + m_compactionRefLayer = m_compactionRefLayerUiField(); } if (m_geomCase->geoMechData()->femPartResults()->assertResultsLoaded(this->resultAddress())) @@ -378,6 +414,8 @@ void RimGeoMechResultDefinition::initAfterRead() m_resultVariableUiField = composeFieldCompString(m_resultFieldName(), m_resultComponentName()); m_isTimeLapseResultUiField = m_isTimeLapseResult; m_timeLapseBaseTimestepUiField = m_timeLapseBaseTimestep; + m_isCompactionUiField = m_isCompaction; + m_compactionRefLayerUiField = m_compactionRefLayer; } @@ -392,16 +430,16 @@ void RimGeoMechResultDefinition::loadResult() } } - //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- RigFemResultAddress RimGeoMechResultDefinition::resultAddress() { - if(m_isTimeLapseResult) - return RigFemResultAddress(resultPositionType(), resultFieldName().toStdString(), resultComponentName().toStdString(), m_timeLapseBaseTimestep); - else - return RigFemResultAddress(resultPositionType(), resultFieldName().toStdString(), resultComponentName().toStdString()); + return RigFemResultAddress(resultPositionType(), + resultFieldName().toStdString(), + resultComponentName().toStdString(), + m_isTimeLapseResult() ? m_timeLapseBaseTimestep() : RigFemResultAddress::NO_TIME_LAPSE, + m_isCompaction() && resultFieldName().toStdString() == RigFemPartResultsCollection::FIELD_NAME_COMPACTION ? m_compactionRefLayer() : RigFemResultAddress::NO_COMPACTION); } //-------------------------------------------------------------------------------------------------- @@ -479,11 +517,15 @@ void RimGeoMechResultDefinition::setResultAddress( const RigFemResultAddress& re m_resultFieldName = QString::fromStdString(resultAddress.fieldName); m_resultComponentName = QString::fromStdString(resultAddress.componentName); m_isTimeLapseResult = resultAddress.isTimeLapse(); - - m_timeLapseBaseTimestep = m_isTimeLapseResult ? resultAddress.timeLapseBaseFrameIdx: 0; + m_isCompaction = resultAddress.isCompaction(); + + m_timeLapseBaseTimestep = m_isTimeLapseResult ? resultAddress.timeLapseBaseFrameIdx: -1; + m_compactionRefLayer = m_isCompaction ? resultAddress.refKLayerIndex : -1; m_resultPositionTypeUiField = m_resultPositionType; m_resultVariableUiField = composeFieldCompString(m_resultFieldName(), m_resultComponentName()); m_isTimeLapseResultUiField = m_isTimeLapseResult; m_timeLapseBaseTimestepUiField = m_timeLapseBaseTimestep; + m_isCompactionUiField = m_isCompaction; + m_compactionRefLayerUiField = m_compactionRefLayer; } diff --git a/ApplicationCode/ProjectDataModel/RimGeoMechResultDefinition.h b/ApplicationCode/ProjectDataModel/RimGeoMechResultDefinition.h index e8e2b3da5b..dd10008d92 100644 --- a/ApplicationCode/ProjectDataModel/RimGeoMechResultDefinition.h +++ b/ApplicationCode/ProjectDataModel/RimGeoMechResultDefinition.h @@ -101,6 +101,8 @@ private: caf::PdmField m_resultComponentName; caf::PdmField m_isTimeLapseResult; caf::PdmField m_timeLapseBaseTimestep; + caf::PdmField m_isCompaction; + caf::PdmField m_compactionRefLayer; // UI Fields only @@ -112,7 +114,8 @@ private: caf::PdmField m_resultVariableUiField; caf::PdmField m_isTimeLapseResultUiField; caf::PdmField m_timeLapseBaseTimestepUiField; - + caf::PdmField m_isCompactionUiField; + caf::PdmField m_compactionRefLayerUiField; caf::PdmPointer m_geomCase; bool m_isChangedByField;