ResInsight/ApplicationCode/GeoMech/OdbReader/RifOdbReader.cpp
2019-09-06 10:57:22 +02:00

950 lines
34 KiB
C++

/////////////////////////////////////////////////////////////////////////////////
//
// 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.
//
/////////////////////////////////////////////////////////////////////////////////
#ifdef _MSC_VER
// Get rid of warnings from compilation of ODB API
#pragma warning( push )
#pragma warning( disable : 4018 )
#pragma warning( disable : 4482 )
#pragma warning( disable : 4584 )
#pragma warning( disable : 4800 )
#endif
#include <odb_API.h>
#include <odb_Enum.h>
#ifdef _MSC_VER
#pragma warning( pop )
#endif
#include "RifOdbReader.h"
#include "RigFemPart.h"
#include "RigFemPartCollection.h"
#include "cafProgressInfo.h"
#include <QString>
#include <iostream>
#include <limits>
#include <map>
#include <sstream>
//==================================================================================================
//
// Helper class to ensure that ODB bulk data are returned as float. Converting if necessary.
//
//==================================================================================================
class RifOdbBulkDataGetter
{
public:
explicit RifOdbBulkDataGetter( const odb_FieldBulkData& bulkData )
: m_bulkData( bulkData ){};
virtual ~RifOdbBulkDataGetter(){};
float* data()
{
odb_Enum::odb_PrecisionEnum precision = m_bulkData.precision();
if ( precision == odb_Enum::SINGLE_PRECISION )
{
return m_bulkData.data();
}
else if ( precision == odb_Enum::DOUBLE_PRECISION )
{
if ( m_data.size() < 1 )
{
int dataSize = m_bulkData.length() * m_bulkData.width();
m_data.resize( dataSize );
double* doublePtr = m_bulkData.dataDouble();
CVF_ASSERT( doublePtr );
float* dataPtr = m_data.data();
for ( int i = 0; i < dataSize; i++ )
{
dataPtr[i] = (float)doublePtr[i];
}
}
return m_data.data();
}
// Should never end up here
CVF_ASSERT( 0 );
return NULL;
}
private:
const odb_FieldBulkData& m_bulkData;
std::vector<float> m_data;
};
size_t RifOdbReader::sm_instanceCount = 0;
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::map<std::string, RigElementType> initFemTypeMap()
{
std::map<std::string, RigElementType> typeMap;
typeMap["C3D8R"] = HEX8;
typeMap["C3D8"] = HEX8;
typeMap["C3D8P"] = HEX8P;
typeMap["CAX4"] = CAX4;
return typeMap;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RigElementType toRigElementType( const odb_String& odbTypeName )
{
static std::map<std::string, RigElementType> odbElmTypeToRigElmTypeMap = initFemTypeMap();
std::map<std::string, RigElementType>::iterator it = odbElmTypeToRigElmTypeMap.find( odbTypeName.cStr() );
if ( it == odbElmTypeToRigElmTypeMap.end() )
{
#if 0
std::cout << "Unsupported element type :" << odbElm.type().cStr() << std::endl;
#endif
return UNKNOWN_ELM_TYPE;
}
return it->second;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const int* localElmNodeToIntegrationPointMapping( RigElementType elmType )
{
static const int HEX8_Mapping[8] = {0, 1, 3, 2, 4, 5, 7, 6};
switch ( elmType )
{
case HEX8:
case HEX8P:
return HEX8_Mapping;
break;
case CAX4:
return HEX8_Mapping; // First four is identical to HEX8
break;
default:
// assert(false); // Element type not supported
break;
}
return NULL;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RifOdbReader::RifOdbReader()
{
if ( ++sm_instanceCount == 1 )
{
odb_initializeAPI();
}
m_odb = NULL;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RifOdbReader::~RifOdbReader()
{
close();
if ( --sm_instanceCount == 0 )
{
odb_finalizeAPI();
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RifOdbReader::close()
{
if ( m_odb )
{
m_odb->close();
m_odb = NULL;
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RifOdbReader::openFile( const std::string& fileName, std::string* errorMessage )
{
close();
CVF_ASSERT( m_odb == NULL );
odb_String path = fileName.c_str();
try
{
m_odb = &openOdb( path, true );
}
catch ( const nex_Exception& nex )
{
if ( errorMessage )
{
*errorMessage = nex.UserReport().CStr();
}
return false;
}
catch ( ... )
{
if ( errorMessage )
{
std::stringstream errStr;
errStr << "Unable to open file '" << fileName << "'.";
*errorMessage = errStr.str();
}
return false;
}
return true;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RifOdbReader::isOpen() const
{
return m_odb != NULL;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RifOdbReader::assertMetaDataLoaded()
{
CVF_ASSERT( m_odb != NULL );
if ( m_resultsMetaData.empty() )
{
m_resultsMetaData = readResultsMetaData( m_odb );
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::map<RifOdbReader::RifOdbResultKey, std::vector<std::string>> RifOdbReader::readResultsMetaData( odb_Odb* odb )
{
CVF_ASSERT( odb != NULL );
std::map<RifOdbResultKey, std::vector<std::string>> resultsMap;
const odb_StepRepository& stepRepository = odb->steps();
odb_StepRepositoryIT stepIt( stepRepository );
stepIt.first();
const odb_Step& step = stepRepository.constGet( stepIt.currentKey() );
const odb_SequenceFrame& stepFrames = step.frames();
if ( stepFrames.size() > 1 )
{
// Optimization: Get results metadata for the second frame of the first step only
const odb_Frame& frame = stepFrames.constGet( 1 );
const odb_FieldOutputRepository& fieldCon = frame.fieldOutputs();
odb_FieldOutputRepositoryIT fieldConIT( fieldCon );
for ( fieldConIT.first(); !fieldConIT.isDone(); fieldConIT.next() )
{
const odb_FieldOutput& field = fieldCon[fieldConIT.currentKey()];
const odb_SequenceFieldLocation& fieldLocations = field.locations();
for ( int loc = 0; loc < fieldLocations.size(); loc++ )
{
const odb_FieldLocation& fieldLocation = fieldLocations.constGet( loc );
std::string fieldName = field.name().CStr();
odb_SequenceString components = field.componentLabels();
std::vector<std::string> compVec;
int numComp = components.size();
for ( int comp = 0; comp < numComp; comp++ )
{
compVec.push_back( components[comp].CStr() );
}
switch ( fieldLocation.position() )
{
case odb_Enum::NODAL:
resultsMap[RifOdbResultKey( NODAL, fieldName )] = compVec;
break;
case odb_Enum::ELEMENT_NODAL:
resultsMap[RifOdbResultKey( ELEMENT_NODAL, fieldName )] = compVec;
break;
case odb_Enum::INTEGRATION_POINT:
resultsMap[RifOdbResultKey( INTEGRATION_POINT, fieldName )] = compVec;
resultsMap[RifOdbResultKey( ELEMENT_NODAL, fieldName )] = compVec;
break;
default:
break;
}
}
}
}
stepFrames.release();
return resultsMap;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RifOdbReader::readFemParts( RigFemPartCollection* femParts )
{
CVF_ASSERT( femParts );
CVF_ASSERT( m_odb != NULL );
odb_InstanceRepository instanceRepository = m_odb->rootAssembly().instances();
odb_InstanceRepositoryIT iter( instanceRepository );
caf::ProgressInfo modelProgress( instanceRepository.size() * ( 2 + 4 ), "Reading Odb Parts" );
int instanceCount = 0;
for ( iter.first(); !iter.isDone(); iter.next(), instanceCount++ )
{
modelProgress.setProgressDescription( QString( iter.currentKey().cStr() ) + ": Reading Nodes" );
m_nodeIdToIdxMaps.push_back( std::map<int, int>() );
odb_Instance& inst = instanceRepository[iter.currentKey()];
RigFemPart* femPart = new RigFemPart;
// Extract nodes
const odb_SequenceNode& odbNodes = inst.nodes();
std::map<int, int>& nodeIdToIdxMap = m_nodeIdToIdxMaps.back();
int nodeCount = odbNodes.size();
femPart->nodes().nodeIds.resize( nodeCount );
femPart->nodes().coordinates.resize( nodeCount );
for ( int nIdx = 0; nIdx < nodeCount; ++nIdx )
{
const odb_Node odbNode = odbNodes.node( nIdx );
femPart->nodes().nodeIds[nIdx] = odbNode.label();
const float* pos = odbNode.coordinates();
femPart->nodes().coordinates[nIdx].set( pos[0], pos[1], pos[2] );
nodeIdToIdxMap[odbNode.label()] = nIdx;
// Progress reporting
if ( nIdx == nodeCount / 2 )
{
modelProgress.incrementProgress();
}
}
modelProgress.incrementProgress();
modelProgress.setProgressDescription( QString( iter.currentKey().cStr() ) + ": Reading Elements" );
// Extract elements
const odb_SequenceElement& elements = inst.elements();
int elmCount = elements.size();
femPart->preAllocateElementStorage( elmCount );
std::vector<int> indexBasedConnectivities;
m_elementIdToIdxMaps.push_back( std::map<int, int>() );
std::map<int, int>& elementIdToIdxMap = m_elementIdToIdxMaps.back();
for ( int elmIdx = 0; elmIdx < elmCount; ++elmIdx )
{
const odb_Element odbElm = elements.element( elmIdx );
elementIdToIdxMap[odbElm.label()] = elmIdx;
RigElementType elmType = toRigElementType( odbElm.type() );
if ( elmType == UNKNOWN_ELM_TYPE ) continue;
int nodeCount = 0;
const int* idBasedConnectivities = odbElm.connectivity( nodeCount );
CVF_TIGHT_ASSERT( nodeCount == RigFemTypes::elmentNodeCount( elmType ) );
indexBasedConnectivities.resize( nodeCount );
for ( int lnIdx = 0; lnIdx < nodeCount; ++lnIdx )
{
indexBasedConnectivities[lnIdx] = nodeIdToIdxMap[idBasedConnectivities[lnIdx]];
}
femPart->appendElement( elmType, odbElm.label(), indexBasedConnectivities.data() );
// Progress reporting
if ( elmIdx == elmCount / 4 || elmIdx == elmCount / 2 || elmIdx == 3 * elmCount / 4 )
{
modelProgress.incrementProgress();
}
}
femPart->setElementPartId( femParts->partCount() );
femParts->addFemPart( femPart );
modelProgress.incrementProgress();
}
return true;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<std::string> RifOdbReader::allStepNames() const
{
CVF_ASSERT( m_odb != NULL );
std::vector<std::string> stepNames;
odb_StepRepository stepRepository = m_odb->steps();
odb_StepRepositoryIT sIter( stepRepository );
for ( sIter.first(); !sIter.isDone(); sIter.next() )
{
std::string stepName( sIter.currentValue().name().CStr() );
stepNames.push_back( stepName );
}
return stepNames;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<std::string> RifOdbReader::filteredStepNames() const
{
CVF_ASSERT( m_odb != NULL );
std::vector<std::string> stepNames;
odb_StepRepository stepRepository = m_odb->steps();
odb_StepRepositoryIT sIter( stepRepository );
int stepIndex = 0;
for ( sIter.first(); !sIter.isDone(); sIter.next() )
{
std::string stepName( sIter.currentValue().name().CStr() );
if ( this->isTimeStepIncludedByFilter( stepIndex++ ) )
{
stepNames.push_back( stepName );
}
}
return stepNames;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<double> RifOdbReader::frameTimes( int stepIndex ) const
{
CVF_ASSERT( m_odb != NULL );
odb_StepRepository& stepRepository = m_odb->steps();
odb_StepList stepList = stepRepository.stepList();
int stepFileIndex = this->timeStepIndexOnFile( stepIndex );
odb_Step& step = stepList.Get( stepFileIndex );
odb_SequenceFrame& stepFrames = step.frames();
std::vector<double> frameValues;
int numFrames = stepFrames.size();
for ( int f = 0; f < numFrames; f++ )
{
odb_Frame frame = stepFrames.constGet( f );
frameValues.push_back( frame.frameValue() );
}
return frameValues;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<std::string> RifOdbReader::elementSetNames( int partIndex )
{
CVF_ASSERT( m_odb != NULL );
std::map<int, std::vector<std::string>>::const_iterator mapIt = m_partElementSetNames.find( partIndex );
if ( mapIt == m_partElementSetNames.end() )
{
std::vector<std::string> setNames;
const odb_Assembly& rootAssembly = m_odb->constRootAssembly();
const odb_InstanceRepository& instances = rootAssembly.instances();
int currentInstance = 0;
odb_InstanceRepositoryIT instIt( instances );
for ( instIt.first(); !instIt.isDone(); instIt.next(), currentInstance++ )
{
const odb_Instance& instance = instIt.currentValue();
if ( currentInstance == partIndex )
{
const odb_SetRepository& sets = rootAssembly.elementSets();
odb_SetRepositoryIT setIt( sets );
for ( setIt.first(); !setIt.isDone(); setIt.next() )
{
const odb_Set& set = setIt.currentValue();
setNames.push_back( set.name().CStr() );
}
break;
}
}
m_partElementSetNames[partIndex] = setNames;
}
return m_partElementSetNames.at( partIndex );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<size_t> RifOdbReader::elementSet( int partIndex, int setIndex )
{
CVF_ASSERT( m_odb != NULL );
std::vector<std::string> setNames = elementSetNames( partIndex );
const odb_Assembly& rootAssembly = m_odb->constRootAssembly();
const odb_Set& set = rootAssembly.elementSets()[odb_String( setNames[setIndex].c_str() )];
odb_SequenceString instanceNames = set.instanceNames();
const odb_SequenceElement& setElements = set.elements( instanceNames[partIndex] );
int elementCount = setElements.size();
std::vector<size_t> elementIndexes;
elementIndexes.resize( elementCount );
for ( int i = 0; i < elementCount; i++ )
{
elementIndexes[i] = setElements.element( i ).index();
}
return elementIndexes;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::map<std::string, std::vector<std::string>> RifOdbReader::scalarNodeFieldAndComponentNames()
{
return fieldAndComponentNames( NODAL );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::map<std::string, std::vector<std::string>> RifOdbReader::scalarElementNodeFieldAndComponentNames()
{
return fieldAndComponentNames( ELEMENT_NODAL );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::map<std::string, std::vector<std::string>> RifOdbReader::scalarIntegrationPointFieldAndComponentNames()
{
return fieldAndComponentNames( INTEGRATION_POINT );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const odb_Frame& RifOdbReader::stepFrame( int stepIndex, int frameIndex ) const
{
CVF_ASSERT( m_odb );
const odb_StepRepository& stepRepository = m_odb->steps();
const odb_StepList& stepList = stepRepository.stepList();
int stepFileIndex = this->timeStepIndexOnFile( stepIndex );
const odb_Step& step = stepList.ConstGet( stepFileIndex );
const odb_SequenceFrame& stepFrames = step.frames();
return stepFrames.constGet( frameIndex );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
odb_Instance* RifOdbReader::instance( int instanceIndex )
{
CVF_ASSERT( m_odb != NULL );
odb_InstanceRepository& instanceRepository = m_odb->rootAssembly().instances();
odb_InstanceRepositoryIT iter( instanceRepository );
int instanceCount = 0;
for ( iter.first(); !iter.isDone(); iter.next(), instanceCount++ )
{
odb_Instance& inst = instanceRepository[iter.currentKey()];
if ( instanceCount == instanceIndex )
{
return &inst;
}
}
return NULL;
}
//--------------------------------------------------------------------------------------------------
/// Get the number of result items (== #nodes or #elements)
//--------------------------------------------------------------------------------------------------
size_t RifOdbReader::resultItemCount( const std::string& fieldName, int partIndex, int stepIndex, int frameIndex )
{
odb_Instance* partInstance = instance( partIndex );
CVF_ASSERT( partInstance != NULL );
const odb_Frame& frame = stepFrame( stepIndex, frameIndex );
const odb_FieldOutput& instanceFieldOutput = frame.fieldOutputs()[fieldName.c_str()].getSubset( *partInstance );
const odb_SequenceFieldBulkData& seqFieldBulkData = instanceFieldOutput.bulkDataBlocks();
size_t resultItemCount = 0;
int numBlocks = seqFieldBulkData.size();
for ( int block = 0; block < numBlocks; block++ )
{
const odb_FieldBulkData& bulkData = seqFieldBulkData[block];
resultItemCount += bulkData.length();
}
return resultItemCount;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
size_t RifOdbReader::componentsCount( const std::string& fieldName, ResultPosition position )
{
std::vector<std::string> compNames = componentNames( RifOdbResultKey( position, fieldName ) );
return compNames.size() > 0 ? compNames.size() : 1;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<std::string> RifOdbReader::componentNames( const RifOdbResultKey& result )
{
assertMetaDataLoaded();
std::map<RifOdbResultKey, std::vector<std::string>>::const_iterator resMapIt = m_resultsMetaData.find( result );
if ( resMapIt != m_resultsMetaData.end() )
{
std::vector<std::string> compNames;
compNames = resMapIt->second;
return compNames;
}
CVF_ASSERT( false );
return std::vector<std::string>();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::map<std::string, std::vector<std::string>> RifOdbReader::fieldAndComponentNames( ResultPosition position )
{
assertMetaDataLoaded();
std::map<std::string, std::vector<std::string>> fieldsAndComponents;
std::map<RifOdbResultKey, std::vector<std::string>>::const_iterator resMapIt;
for ( resMapIt = m_resultsMetaData.begin(); resMapIt != m_resultsMetaData.end(); ++resMapIt )
{
if ( resMapIt->first.resultPostion == position )
{
fieldsAndComponents[resMapIt->first.fieldName] = resMapIt->second;
}
}
return fieldsAndComponents;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RifOdbReader::readDisplacements( int partIndex, int stepIndex, int frameIndex, std::vector<cvf::Vec3f>* displacements )
{
CVF_ASSERT( displacements );
odb_Instance* partInstance = instance( partIndex );
CVF_ASSERT( partInstance != NULL );
size_t dataSize = resultItemCount( "U", partIndex, stepIndex, frameIndex );
if ( dataSize > 0 )
{
displacements->resize( dataSize );
}
const odb_Frame& frame = stepFrame( stepIndex, frameIndex );
const odb_FieldOutput& instanceFieldOutput = frame.fieldOutputs()["U"].getSubset( *partInstance );
const odb_SequenceFieldBulkData& seqFieldBulkData = instanceFieldOutput.bulkDataBlocks();
size_t dataIndex = 0;
int numBlocks = seqFieldBulkData.size();
for ( int block = 0; block < numBlocks; block++ )
{
const odb_FieldBulkData& bulkData = seqFieldBulkData[block];
RifOdbBulkDataGetter bulkDataGetter( bulkData );
if ( bulkData.numberOfNodes() > 0 )
{
int numNodes = bulkData.length();
int numComp = bulkData.width();
float* data = bulkDataGetter.data();
for ( int i = 0; i < numNodes; i++ )
{
( *displacements )[i + dataIndex].set( data[i * numComp], data[i * numComp + 1], data[i * numComp + 2] );
}
dataIndex += numNodes * numComp;
}
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RifOdbReader::readNodeField( const std::string& fieldName,
int partIndex,
int stepIndex,
int frameIndex,
std::vector<std::vector<float>*>* resultValues )
{
CVF_ASSERT( resultValues );
odb_Instance* partInstance = instance( partIndex );
CVF_ASSERT( partInstance != NULL );
size_t compCount = componentsCount( fieldName, NODAL );
CVF_ASSERT( compCount == resultValues->size() );
std::map<int, int>& nodeIdToIdxMap = m_nodeIdToIdxMaps[partIndex];
size_t dataSize = nodeIdToIdxMap.size();
if ( dataSize > 0 )
{
for ( int comp = 0; comp < compCount; comp++ )
{
CVF_ASSERT( ( *resultValues )[comp] );
( *resultValues )[comp]->resize( dataSize, std::numeric_limits<float>::infinity() );
}
}
const odb_Frame& frame = stepFrame( stepIndex, frameIndex );
const odb_FieldOutput& instanceFieldOutput = frame.fieldOutputs()[fieldName.c_str()].getSubset( *partInstance );
const odb_FieldOutput& fieldOutput = instanceFieldOutput.getSubset( odb_Enum::NODAL );
const odb_SequenceFieldBulkData& seqFieldBulkData = fieldOutput.bulkDataBlocks();
int numBlocks = seqFieldBulkData.size();
for ( int block = 0; block < numBlocks; block++ )
{
const odb_FieldBulkData& bulkData = seqFieldBulkData[block];
RifOdbBulkDataGetter bulkDataGetter( bulkData );
int numNodes = bulkData.length();
int numComp = bulkData.width();
int* nodeLabels = bulkData.nodeLabels();
float* data = bulkDataGetter.data();
for ( int nIdx = 0; nIdx < numNodes; nIdx++ )
{
for ( int comp = 0; comp < numComp; comp++ )
{
std::vector<float>* singleComponentValues = ( *resultValues )[comp];
( *singleComponentValues )[nodeIdToIdxMap[nodeLabels[nIdx]]] = data[nIdx * numComp + comp];
}
}
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RifOdbReader::readElementNodeField( const std::string& fieldName,
int partIndex,
int stepIndex,
int frameIndex,
std::vector<std::vector<float>*>* resultValues )
{
CVF_ASSERT( resultValues );
odb_Instance* partInstance = instance( partIndex );
CVF_ASSERT( partInstance != NULL );
size_t compCount = componentsCount( fieldName, ELEMENT_NODAL );
CVF_ASSERT( compCount == resultValues->size() );
size_t dataSize = resultItemCount( fieldName, partIndex, stepIndex, frameIndex );
if ( dataSize > 0 )
{
for ( int comp = 0; comp < compCount; comp++ )
{
CVF_ASSERT( ( *resultValues )[comp] );
( *resultValues )[comp]->resize( dataSize, std::numeric_limits<float>::infinity() );
}
}
const odb_Frame& frame = stepFrame( stepIndex, frameIndex );
const odb_FieldOutput& instanceFieldOutput = frame.fieldOutputs()[fieldName.c_str()].getSubset( *partInstance );
const odb_FieldOutput& fieldOutput = instanceFieldOutput.getSubset( odb_Enum::ELEMENT_NODAL );
const odb_SequenceFieldBulkData& seqFieldBulkData = fieldOutput.bulkDataBlocks();
std::map<int, int>& elementIdToIdxMap = m_elementIdToIdxMaps[partIndex];
CVF_ASSERT( elementIdToIdxMap.size() > 0 );
int numBlocks = seqFieldBulkData.size();
for ( int block = 0; block < numBlocks; block++ )
{
const odb_FieldBulkData& bulkData = seqFieldBulkData[block];
RifOdbBulkDataGetter bulkDataGetter( bulkData );
int numValues = bulkData.length();
int numComp = bulkData.width();
int elemCount = bulkData.numberOfElements();
int elemNodeCount = numValues / elemCount;
int* elementLabels = bulkData.elementLabels();
float* data = bulkDataGetter.data();
for ( int elem = 0; elem < elemCount; elem++ )
{
int elementIdx = elementIdToIdxMap[elementLabels[elem * elemNodeCount]];
int elementResultStartDestIdx = elementIdx * elemNodeCount; // Ikke generellt riktig !
int elementResultStartSourceIdx = elem * elemNodeCount * numComp;
for ( int elemNode = 0; elemNode < elemNodeCount; elemNode++ )
{
int destIdx = elementResultStartDestIdx + elemNode;
int srcIdx = elementResultStartSourceIdx + elemNode * numComp;
for ( int comp = 0; comp < numComp; comp++ )
{
( *( *resultValues )[comp] )[destIdx] = data[srcIdx + comp];
}
}
}
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RifOdbReader::readIntegrationPointField( const std::string& fieldName,
int partIndex,
int stepIndex,
int frameIndex,
std::vector<std::vector<float>*>* resultValues )
{
CVF_ASSERT( resultValues );
odb_Instance* partInstance = instance( partIndex );
CVF_ASSERT( partInstance != NULL );
size_t compCount = componentsCount( fieldName, INTEGRATION_POINT );
CVF_ASSERT( compCount == resultValues->size() );
size_t dataSize = resultItemCount( fieldName, partIndex, stepIndex, frameIndex );
if ( dataSize > 0 )
{
for ( int comp = 0; comp < compCount; comp++ )
{
CVF_ASSERT( ( *resultValues )[comp] );
( *resultValues )[comp]->resize( dataSize, std::numeric_limits<float>::infinity() );
}
}
const odb_Frame& frame = stepFrame( stepIndex, frameIndex );
const odb_FieldOutput& instanceFieldOutput = frame.fieldOutputs()[fieldName.c_str()].getSubset( *partInstance );
const odb_FieldOutput& fieldOutput = instanceFieldOutput.getSubset( odb_Enum::INTEGRATION_POINT );
const odb_SequenceFieldBulkData& seqFieldBulkData = fieldOutput.bulkDataBlocks();
std::map<int, int>& elementIdToIdxMap = m_elementIdToIdxMaps[partIndex];
CVF_ASSERT( elementIdToIdxMap.size() > 0 );
int numBlocks = seqFieldBulkData.size();
for ( int block = 0; block < numBlocks; block++ )
{
const odb_FieldBulkData& bulkData = seqFieldBulkData[block];
RifOdbBulkDataGetter bulkDataGetter( bulkData );
int numValues = bulkData.length();
int numComp = bulkData.width();
int elemCount = bulkData.numberOfElements();
int ipCount = numValues / elemCount;
int* elementLabels = bulkData.elementLabels();
float* data = bulkDataGetter.data();
RigElementType eType = toRigElementType( bulkData.baseElementType() );
const int* elmNodeToIpResultMapping = localElmNodeToIntegrationPointMapping(
eType ); // Todo: Use the one in RigFemTypes.h, but we need to guard against unknown element types first.
if ( !elmNodeToIpResultMapping ) continue;
for ( int elem = 0; elem < elemCount; elem++ )
{
int elementIdx = elementIdToIdxMap[elementLabels[elem * ipCount]];
int elementResultStartDestIdx = elementIdx * ipCount; // Ikke generellt riktig !
int elementResultStartSourceIdx = elem * ipCount * numComp;
for ( int ipIdx = 0; ipIdx < ipCount; ipIdx++ )
{
int resultIpIdx = elmNodeToIpResultMapping[ipIdx];
int destIdx = elementResultStartDestIdx + ipIdx;
int srcIdx = elementResultStartSourceIdx + resultIpIdx * numComp;
for ( int comp = 0; comp < numComp; comp++ )
{
( *( *resultValues )[comp] )[destIdx] = data[srcIdx + comp];
}
}
}
}
}