ResInsight/ApplicationLibCode/GeoMech/GeoMechVisualization/RivFemElmVisibilityCalculator.cpp
jonjenssen d09ae4e1cb
Improved ODB support (#8046)
* Experiments for supporting visualization of new ODB files from WIA workflow

* Some more experiments to get odb working for wia results

* More work in progress, experimenting to get wellIA result files to load properly

* Make sure all part geometries use the same global bounding box

* Clean up code

* Add some safeguards for data calculations
Move parts below grid in project tree

* Fix warnings

* Add support for C3D8RT elements
Add some more safeguards for missing data
Remove strange part handling

* Support elements with reduced number of integration points by pretending to have 8.

* Change integration point mapping to correct order (ref. Stein and Abaqus 2019 doc)

* Do not allocate too much memory for element nodal results for 20 element node types

* Code cleanup.
Revert back to old integration point numbering scheme (ref. Stein)

* And, another integration point order update...

* Update comments
2021-09-27 12:44:29 +02:00

282 lines
12 KiB
C++

/////////////////////////////////////////////////////////////////////////////////
//
// Copyright (C) 2015- Statoil ASA
// Copyright (C) 2015- 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 "RivFemElmVisibilityCalculator.h"
#include "RigCaseToCaseCellMapper.h"
#include "RigFemPart.h"
#include "RigFemPartGrid.h"
#include "RigFemPartResultsCollection.h"
#include "RigGeoMechCaseData.h"
#include "RimEclipsePropertyFilter.h"
#include "RimGeoMechCase.h"
#include "RimGeoMechPropertyFilter.h"
#include "RimGeoMechPropertyFilterCollection.h"
#include "RimGeoMechResultDefinition.h"
#include "RimGeoMechView.h"
#include "RimViewController.h"
#include "RimViewLinker.h"
#include "cvfStructGrid.h"
#include "cvfStructGridGeometryGenerator.h"
#include <cmath>
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RivFemElmVisibilityCalculator::computeAllVisible( cvf::UByteArray* elmVisibilities, const RigFemPart* femPart )
{
elmVisibilities->resize( femPart->elementCount() );
elmVisibilities->setAll( true );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RivFemElmVisibilityCalculator::computeRangeVisibility( cvf::UByteArray* elmVisibilities,
const RigFemPart* femPart,
const cvf::CellRangeFilter& rangeFilter )
{
elmVisibilities->resize( femPart->elementCount() );
const RigFemPartGrid* grid = femPart->getOrCreateStructGrid();
if ( rangeFilter.hasIncludeRanges() )
{
for ( int elmIdx = 0; elmIdx < femPart->elementCount(); ++elmIdx )
{
size_t mainGridI;
size_t mainGridJ;
size_t mainGridK;
grid->ijkFromCellIndex( elmIdx, &mainGridI, &mainGridJ, &mainGridK );
( *elmVisibilities )[elmIdx] = rangeFilter.isCellVisible( mainGridI, mainGridJ, mainGridK, false );
}
}
else
{
for ( int elmIdx = 0; elmIdx < femPart->elementCount(); ++elmIdx )
{
size_t mainGridI;
size_t mainGridJ;
size_t mainGridK;
grid->ijkFromCellIndex( elmIdx, &mainGridI, &mainGridJ, &mainGridK );
( *elmVisibilities )[elmIdx] = !rangeFilter.isCellExcluded( mainGridI, mainGridJ, mainGridK, false );
}
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RivFemElmVisibilityCalculator::computePropertyVisibility( cvf::UByteArray* cellVisibility,
const RigFemPart* part,
int timeStepIndex,
const cvf::UByteArray* rangeFilterVisibility,
RimGeoMechPropertyFilterCollection* propFilterColl )
{
CVF_ASSERT( cellVisibility != nullptr );
CVF_ASSERT( rangeFilterVisibility != nullptr );
CVF_ASSERT( propFilterColl != nullptr );
CVF_ASSERT( part->elementCount() > 0 );
CVF_ASSERT( rangeFilterVisibility->size() == static_cast<size_t>( part->elementCount() ) );
// Copy if not equal
if ( cellVisibility != rangeFilterVisibility ) ( *cellVisibility ) = *rangeFilterVisibility;
const int elementCount = part->elementCount();
if ( !propFilterColl->hasActiveFilters() ) return;
for ( size_t i = 0; i < propFilterColl->propertyFilters().size(); i++ )
{
RimGeoMechPropertyFilter* propertyFilter = propFilterColl->propertyFilters()[i];
if ( !propertyFilter->isActiveAndHasResult() ) continue;
const RimCellFilter::FilterModeType filterType = propertyFilter->filterMode();
RigGeoMechCaseData* caseData = propFilterColl->reservoirView()->geoMechCase()->geoMechData();
RigFemResultAddress resVarAddress = propertyFilter->resultDefinition->resultAddress();
// Do a "Hack" to use elm nodal and not nodal POR results
if ( resVarAddress.resultPosType == RIG_NODAL && resVarAddress.fieldName == "POR-Bar" )
resVarAddress.resultPosType = RIG_ELEMENT_NODAL;
const std::vector<float>& resVals =
caseData->femPartResults()->resultValues( resVarAddress, part->elementPartId(), timeStepIndex );
if ( !propertyFilter->isActive() ) continue;
if ( !propertyFilter->resultDefinition->hasResult() ) continue;
if ( resVals.size() == 0 ) continue;
const double lowerBound = propertyFilter->lowerBound();
const double upperBound = propertyFilter->upperBound();
if ( propertyFilter->resultDefinition->resultAddress().resultPosType == RIG_FORMATION_NAMES )
{
std::vector<int> integerVector = propertyFilter->selectedCategoryValues();
std::set<int> integerSet;
for ( auto val : integerVector )
{
integerSet.insert( val );
}
for ( int cellIndex = 0; cellIndex < elementCount; cellIndex++ )
{
if ( !( *cellVisibility )[cellIndex] ) continue;
size_t resultValueIndex = part->elementNodeResultIdx( cellIndex, 0 );
double scalarValue = resVals[resultValueIndex];
int intValue = nearbyint( scalarValue );
if ( integerSet.find( intValue ) != integerSet.end() )
{
if ( filterType == RimCellFilter::EXCLUDE )
{
( *cellVisibility )[cellIndex] = false;
}
}
else
{
if ( filterType == RimCellFilter::INCLUDE )
{
( *cellVisibility )[cellIndex] = false;
}
}
}
}
else if ( resVarAddress.resultPosType == RIG_ELEMENT )
{
#pragma omp parallel for schedule( dynamic )
for ( int cellIndex = 0; cellIndex < elementCount; cellIndex++ )
{
if ( !( *cellVisibility )[cellIndex] ) continue;
double scalarValue = resVals[cellIndex];
evaluateAndSetCellVisibiliy( cellIndex, scalarValue, lowerBound, upperBound, filterType, cellVisibility );
}
}
else if ( resVarAddress.resultPosType == RIG_ELEMENT_NODAL_FACE )
{
#pragma omp parallel for schedule( dynamic )
for ( int cellIndex = 0; cellIndex < elementCount; cellIndex++ )
{
if ( !( *cellVisibility )[cellIndex] ) continue;
for ( int fpIdx = 0; fpIdx < 24; ++fpIdx )
{
double scalarValue = resVals[cellIndex * 24 + fpIdx];
evaluateAndSetCellVisibiliy( cellIndex, scalarValue, lowerBound, upperBound, filterType, cellVisibility );
}
}
}
else
{
#pragma omp parallel for schedule( dynamic )
for ( int cellIndex = 0; cellIndex < elementCount; cellIndex++ )
{
if ( !( *cellVisibility )[cellIndex] ) continue;
RigElementType eType = part->elementType( cellIndex );
int elmNodeCount = RigFemTypes::elementNodeCount( eType );
const int* elmNodeIndices = part->connectivities( cellIndex );
for ( int enIdx = 0; enIdx < elmNodeCount; ++enIdx )
{
size_t resultValueIndex = cvf::UNDEFINED_SIZE_T;
if ( resVarAddress.resultPosType == RIG_NODAL )
{
resultValueIndex = elmNodeIndices[enIdx];
}
else
{
resultValueIndex = part->elementNodeResultIdx( cellIndex, enIdx );
}
double scalarValue = resVals[resultValueIndex];
evaluateAndSetCellVisibiliy( cellIndex, scalarValue, lowerBound, upperBound, filterType, cellVisibility );
}
}
}
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RivFemElmVisibilityCalculator::evaluateAndSetCellVisibiliy( int cellIndex,
double scalarValue,
double lowerBound,
double upperBound,
const RimCellFilter::FilterModeType filterType,
cvf::UByteArray* cellVisibility )
{
if ( lowerBound <= scalarValue && scalarValue <= upperBound )
{
if ( filterType == RimCellFilter::EXCLUDE )
{
( *cellVisibility )[cellIndex] = false;
}
}
else
{
if ( filterType == RimCellFilter::INCLUDE )
{
( *cellVisibility )[cellIndex] = false;
}
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RivFemElmVisibilityCalculator::computeOverriddenCellVisibility( cvf::UByteArray* elmVisibilities,
const RigFemPart* femPart,
RimViewController* masterViewLink )
{
CVF_ASSERT( elmVisibilities != nullptr );
CVF_ASSERT( femPart != nullptr );
RimGridView* masterView = masterViewLink->ownerViewLinker()->masterView();
cvf::ref<cvf::UByteArray> totCellVisibility = masterView->currentTotalCellVisibility();
int elmCount = femPart->elementCount();
elmVisibilities->resize( elmCount );
elmVisibilities->setAll( false );
const RigCaseToCaseCellMapper* cellMapper = masterViewLink->cellMapper();
for ( int elmIdx = 0; elmIdx < elmCount; ++elmIdx )
{
// We are assuming that there is only one part.
int cellCount = 0;
const int* cellIndicesInMasterCase = cellMapper->masterCaseCellIndices( elmIdx, &cellCount );
for ( int mcIdx = 0; mcIdx < cellCount; ++mcIdx )
{
( *elmVisibilities )[elmIdx] |= ( *totCellVisibility )[cellIndicesInMasterCase[mcIdx]]; // If any is
// visible, show
}
}
}