ResInsight/ApplicationLibCode/UserInterface/RiuFemTimeHistoryResultAccessor.cpp

/////////////////////////////////////////////////////////////////////////////////
//
//  Copyright (C) Statoil ASA
//  Copyright (C) Ceetron Solutions AS
//
//  ResInsight is free software: you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation, either version 3 of the License, or
//  (at your option) any later version.
//
//  ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY
//  WARRANTY; without even the implied warranty of MERCHANTABILITY or
//  FITNESS FOR A PARTICULAR PURPOSE.
//
//  See the GNU General Public License at <http://www.gnu.org/licenses/gpl.html>
//  for more details.
//
/////////////////////////////////////////////////////////////////////////////////

#include "RiuFemTimeHistoryResultAccessor.h"

#include "RiaNumberFormat.h"

#include "RigFemClosestResultIndexCalculator.h"
#include "RigFemPart.h"
#include "RigFemPartCollection.h"
#include "RigFemPartGrid.h"
#include "RigFemPartResultsCollection.h"
#include "RigFemTypes.h"
#include "RigGeoMechCaseData.h"

#include "RiuGeoMechXfTensorResultAccessor.h"

#include <cmath> // Needed for HUGE_VAL on Linux

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RiuFemTimeHistoryResultAccessor::RiuFemTimeHistoryResultAccessor( RigGeoMechCaseData* geomData,
                                                                  RigFemResultAddress femResultAddress,
                                                                  size_t              gridIndex,
                                                                  int                 elementIndex,
                                                                  int                 face,
                                                                  const cvf::Vec3d&   intersectionPointInDomain )
    : m_geoMechCaseData( geomData )
    , m_femResultAddress( new RigFemResultAddress( femResultAddress ) )
    , m_gridIndex( gridIndex )
    , m_elementIndex( elementIndex )
    , m_face( face )
    , m_intersectionPointInDomain( intersectionPointInDomain )
    , m_hasIntersectionTriangle( false )
{
    computeTimeHistoryData();
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RiuFemTimeHistoryResultAccessor::RiuFemTimeHistoryResultAccessor( RigGeoMechCaseData*              geomData,
                                                                  RigFemResultAddress              femResultAddress,
                                                                  size_t                           gridIndex,
                                                                  int                              elementIndex,
                                                                  int                              face,
                                                                  const cvf::Vec3d&                intersectionPointInDomain,
                                                                  const std::array<cvf::Vec3f, 3>& intersectionTriangle )
    : m_geoMechCaseData( geomData )
    , m_femResultAddress( new RigFemResultAddress( femResultAddress ) )
    , m_gridIndex( gridIndex )
    , m_elementIndex( elementIndex )
    , m_face( face )
    , m_intersectionPointInDomain( intersectionPointInDomain )
    , m_hasIntersectionTriangle( true )
    , m_intersectionTriangle( intersectionTriangle )
{
    computeTimeHistoryData();
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
QString RiuFemTimeHistoryResultAccessor::geometrySelectionText() const
{
    QString text;

    if ( m_geoMechCaseData )
    {
        RigFemPart* femPart   = m_geoMechCaseData->femParts()->part( m_gridIndex );
        int         elementId = femPart->elmId( m_elementIndex );
        text += QString( "Element : Id[%1]" ).arg( elementId );

        size_t i = 0;
        size_t j = 0;
        size_t k = 0;
        if ( m_geoMechCaseData->femParts()->part( m_gridIndex )->getOrCreateStructGrid()->ijkFromCellIndex( m_elementIndex, &i, &j, &k ) )
        {
            // Adjust to 1-based Eclipse indexing
            i++;
            j++;
            k++;

            cvf::Vec3d domainCoord = m_intersectionPointInDomain;
            text += QString( ", ijk[%1, %2, %3] " ).arg( i ).arg( j ).arg( k );

            auto xTxt = RiaNumberFormat::valueToText( domainCoord.x(), RiaNumberFormat::NumberFormatType::FIXED, 2 );
            auto yTxt = RiaNumberFormat::valueToText( domainCoord.y(), RiaNumberFormat::NumberFormatType::FIXED, 2 );
            auto zTxt = RiaNumberFormat::valueToText( -domainCoord.z(), RiaNumberFormat::NumberFormatType::FIXED, 2 );

            QString formattedText = QString( "Intersection point : [E: %1, N: %2, Depth: %3]" ).arg( xTxt ).arg( yTxt ).arg( zTxt );

            text += formattedText;
        }
    }

    return text;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<double> RiuFemTimeHistoryResultAccessor::timeHistoryValues() const
{
    return m_timeHistoryValues;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
int RiuFemTimeHistoryResultAccessor::closestNodeId() const
{
    return m_closestNodeId;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RiuFemTimeHistoryResultAccessor::computeTimeHistoryData()
{
    m_timeHistoryValues.clear();

    RigFemClosestResultIndexCalculator closestCalc( m_geoMechCaseData->femParts()->part( m_gridIndex ),
                                                    m_femResultAddress->resultPosType,
                                                    m_elementIndex,
                                                    m_face,
                                                    m_intersectionPointInDomain );

    int scalarResultIndex = closestCalc.resultIndexToClosestResult();
    m_closestNodeId       = closestCalc.closestNodeId();

    RigFemPartResultsCollection* femPartResultsColl = m_geoMechCaseData->femPartResults();

    if ( m_femResultAddress->resultPosType == RIG_ELEMENT_NODAL_FACE && m_hasIntersectionTriangle )
    {
        int closestElmNodeResIndex = closestCalc.closestElementNodeResIdx();

        const int timeSteps = femPartResultsColl->timeStepCount();
        for ( int stepIdx = 0; stepIdx < timeSteps; stepIdx++ )
        {
            const int frameCount = femPartResultsColl->frameCount( stepIdx );
            for ( int frameIdx = 0; frameIdx < frameCount; frameIdx++ )
            {
                RiuGeoMechXfTensorResultAccessor stressXfAccessor( femPartResultsColl,
                                                                   *m_femResultAddress,
                                                                   static_cast<int>( m_gridIndex ),
                                                                   stepIdx,
                                                                   frameIdx );
                float scalarValue = stressXfAccessor.calculateElmNodeValue( m_intersectionTriangle, closestElmNodeResIndex );
                m_timeHistoryValues.push_back( scalarValue );
            }
        }
    }
    else
    {
        if ( scalarResultIndex < 0 ) return;

        const int timeSteps = femPartResultsColl->timeStepCount();
        for ( int stepIdx = 0; stepIdx < timeSteps; stepIdx++ )
        {
            const int frameCount = femPartResultsColl->frameCount( stepIdx );
            for ( int frameIdx = 0; frameIdx < frameCount; frameIdx++ )
            {
                const std::vector<float>& scalarResults =
                    m_geoMechCaseData->femPartResults()->resultValues( *m_femResultAddress, static_cast<int>( m_gridIndex ), stepIdx, frameIdx );
                if ( !scalarResults.empty() )
                {
                    float scalarValue = scalarResults[scalarResultIndex];

                    m_timeHistoryValues.push_back( scalarValue );
                }
            }
        }
    }
}