Added riGetCellCorners/riGetActiveCellCorners

Updated how matrix is allocated in Octave for riGetCellCenters
p4#: 21697
This commit is contained in:
Magne Sjaastad 2013-05-23 11:02:10 +02:00
parent e4bc6c9339
commit d5fd09f658
5 changed files with 558 additions and 18 deletions

View File

@ -346,32 +346,36 @@ public:
if (!rimCase || !rimCase->reservoirData() || (argGridIndex >= rimCase->reservoirData()->gridCount()) ) if (!rimCase || !rimCase->reservoirData() || (argGridIndex >= rimCase->reservoirData()->gridCount()) )
{ {
// No data available // No data available
socketStream << (quint64)0 << (quint64)0 ; socketStream << (quint64)0 << (quint64)0 << (quint64)0 << (quint64)0 << (quint64)0;
return true; return true;
} }
RigGridBase* rigGrid = rimCase->reservoirData()->grid(argGridIndex); RigGridBase* rigGrid = rimCase->reservoirData()->grid(argGridIndex);
quint64 cellCount = (quint64)rigGrid->cellCount(); quint64 cellCount = (quint64)rigGrid->cellCount();
std::vector<double> cellCenterValues(cellCount * 3); quint64 cellCountI = (quint64)rigGrid->cellCountI();
quint64 cellCountJ = (quint64)rigGrid->cellCountJ();
quint64 cellCountK = (quint64)rigGrid->cellCountK();
socketStream << cellCount;
socketStream << cellCountI;
socketStream << cellCountJ;
socketStream << cellCountK;
quint64 byteCount = cellCount * 3 * sizeof(double);
socketStream << byteCount;
for (size_t i = 0; i < cellCount; i++) for (size_t i = 0; i < cellCount; i++)
{ {
cvf::Vec3d center = rigGrid->cell(i).center(); cvf::Vec3d center = rigGrid->cell(i).center();
cellCenterValues[i * 3 + 0] = center.x(); socketStream << center.x();
cellCenterValues[i * 3 + 1] = center.y(); socketStream << center.y();
cellCenterValues[i * 3 + 2] = center.z(); socketStream << center.z();
} }
socketStream << cellCount;
quint64 byteCount = cellCount * 3 * sizeof(double);
socketStream << byteCount;
server->currentClient()->write((const char *)cellCenterValues.data(), byteCount);
return true; return true;
} }
}; };
static bool RiaGetCellCenters_init = RiaSocketCommandFactory::instance()->registerCreator<RiaGetCellCenters>(RiaGetCellCenters::commandName()); static bool RiaGetCellCenters_init = RiaSocketCommandFactory::instance()->registerCreator<RiaGetCellCenters>(RiaGetCellCenters::commandName());
@ -468,3 +472,175 @@ public:
}; };
static bool RiaGetActiveCellCenters_init = RiaSocketCommandFactory::instance()->registerCreator<RiaGetActiveCellCenters>(RiaGetActiveCellCenters::commandName()); static bool RiaGetActiveCellCenters_init = RiaSocketCommandFactory::instance()->registerCreator<RiaGetActiveCellCenters>(RiaGetActiveCellCenters::commandName());
class RiaGetCellCorners : public RiaSocketCommand
{
public:
static QString commandName () { return QString("GetCellCorners"); }
virtual bool interpretCommand(RiaSocketServer* server, const QList<QByteArray>& args, QDataStream& socketStream)
{
int argCaseGroupId = -1;
size_t argGridIndex = 0;
if (args.size() == 2)
{
argGridIndex = args[1].toInt();
}
else if (args.size() == 3)
{
argCaseGroupId = args[1].toInt();
argGridIndex = args[2].toUInt();
}
RimCase* rimCase = server->findReservoir(argCaseGroupId);
if (!rimCase || !rimCase->reservoirData() || (argGridIndex >= rimCase->reservoirData()->gridCount()) )
{
// No data available
socketStream << (quint64)0 << (quint64)0 << (quint64)0 << (quint64)0 << (quint64)0;
return true;
}
RigGridBase* rigGrid = rimCase->reservoirData()->grid(argGridIndex);
quint64 cellCount = (quint64)rigGrid->cellCount();
quint64 cellCountI = (quint64)rigGrid->cellCountI();
quint64 cellCountJ = (quint64)rigGrid->cellCountJ();
quint64 cellCountK = (quint64)rigGrid->cellCountK();
size_t doubleValueCount = cellCount * 3 * 8;
quint64 byteCount = doubleValueCount * sizeof(double);
socketStream << cellCount;
socketStream << cellCountI;
socketStream << cellCountJ;
socketStream << cellCountK;
socketStream << byteCount;
std::vector<double> cellCornerValues(doubleValueCount);
cvf::Vec3d cornerVerts[8];
for (size_t localGridCellIdx = 0; localGridCellIdx < rigGrid->cellCount(); localGridCellIdx++)
{
rigGrid->cellCornerVertices(localGridCellIdx, cornerVerts);
for (size_t j = 0; j < 8; j++)
{
/*
cellCornerValues[localGridCellIdx * 3 + 0] = cornerVerts[j].x();
cellCornerValues[localGridCellIdx * 3 + 1] = cornerVerts[j].y();
cellCornerValues[localGridCellIdx * 3 + 2] = cornerVerts[j].z();
*/
socketStream << cornerVerts[j].x();
socketStream << cornerVerts[j].y();
socketStream << cornerVerts[j].z();
}
}
//server->currentClient()->write((const char *)cellCornerValues.data(), byteCount);
return true;
}
};
static bool RiaGetCellCorners_init = RiaSocketCommandFactory::instance()->registerCreator<RiaGetCellCorners>(RiaGetCellCorners::commandName());
class RiaGetActiveCellCorners : public RiaSocketCommand
{
public:
static QString commandName () { return QString("GetActiveCellCorners"); }
virtual bool interpretCommand(RiaSocketServer* server, const QList<QByteArray>& args, QDataStream& socketStream)
{
int argCaseGroupId = -1;
size_t argGridIndex = 0;
QString porosityModelName;
if (args.size() == 2)
{
argGridIndex = args[1].toInt();
}
else if (args.size() == 3)
{
bool numberConversionOk = false;
int tmpValue = args[2].toInt(&numberConversionOk);
if (numberConversionOk)
{
// Two arguments, caseID and gridIndex
argCaseGroupId = args[1].toInt();
argGridIndex = args[2].toUInt();
}
else
{
// Two arguments, gridIndex and porosity model
argGridIndex = args[1].toUInt();
porosityModelName = args[2];
}
}
else if (args.size() > 3)
{
// Two arguments, caseID and gridIndex
argCaseGroupId = args[1].toInt();
argGridIndex = args[2].toUInt();
porosityModelName = args[3];
}
RifReaderInterface::PorosityModelResultType porosityModelEnum = RifReaderInterface::MATRIX_RESULTS;
if (porosityModelName.toUpper() == "FRACTURE")
{
porosityModelEnum = RifReaderInterface::FRACTURE_RESULTS;
}
RimCase* rimCase = server->findReservoir(argCaseGroupId);
if (!rimCase || !rimCase->reservoirData() || (argGridIndex >= rimCase->reservoirData()->gridCount()) )
{
// No data available
socketStream << (quint64)0 << (quint64)0 ;
return true;
}
RigActiveCellInfo* actCellInfo = rimCase->reservoirData()->activeCellInfo(porosityModelEnum);
RigGridBase* rigGrid = rimCase->reservoirData()->grid(argGridIndex);
std::vector<double> cellCornerCoords(rigGrid->cellCount() * 3 * 8);
cvf::Vec3d cornerVerts[8];
quint64 coordCount = 0;
for (size_t localGridCellIdx = 0; localGridCellIdx < rigGrid->cellCount(); localGridCellIdx++)
{
size_t globalCellIdx = rigGrid->globalGridCellIndex(localGridCellIdx);
if (!actCellInfo->isActive(globalCellIdx)) continue;
rigGrid->cellCornerVertices(localGridCellIdx, cornerVerts);
for (size_t j = 0; j < 8; j++)
{
cellCornerCoords[coordCount * 3 + 0] = cornerVerts[j].x();
cellCornerCoords[coordCount * 3 + 1] = cornerVerts[j].y();
cellCornerCoords[coordCount * 3 + 2] = cornerVerts[j].z();
coordCount++;
}
}
cellCornerCoords.resize(coordCount * 3);
socketStream << coordCount;
quint64 byteCount = coordCount * 3 * sizeof(double);
socketStream << byteCount;
server->currentClient()->write((const char *)cellCornerCoords.data(), byteCount);
return true;
}
};
static bool RiaGetActiveCellCorners_init = RiaSocketCommandFactory::instance()->registerCreator<RiaGetActiveCellCorners>(RiaGetActiveCellCorners::commandName());

View File

@ -15,6 +15,8 @@ set(CPP_SOURCES
riGetCoarseningInfo.cpp riGetCoarseningInfo.cpp
riGetCellCenters.cpp riGetCellCenters.cpp
riGetActiveCellCenters.cpp riGetActiveCellCenters.cpp
riGetCellCorners.cpp
riGetActiveCellCorners.cpp
) )
if (${CMAKE_SYSTEM_NAME} MATCHES "Linux") if (${CMAKE_SYSTEM_NAME} MATCHES "Linux")
@ -119,6 +121,8 @@ else()
"${CMAKE_CURRENT_BINARY_DIR}/riGetCoarseningInfo.oct" "${CMAKE_CURRENT_BINARY_DIR}/riGetCoarseningInfo.oct"
"${CMAKE_CURRENT_BINARY_DIR}/riGetCellCenters.oct" "${CMAKE_CURRENT_BINARY_DIR}/riGetCellCenters.oct"
"${CMAKE_CURRENT_BINARY_DIR}/riGetActiveCellCenters.oct" "${CMAKE_CURRENT_BINARY_DIR}/riGetActiveCellCenters.oct"
"${CMAKE_CURRENT_BINARY_DIR}/riGetCellCorners.oct"
"${CMAKE_CURRENT_BINARY_DIR}/riGetActiveCellCorners.oct"
SOURCES ${CPP_SOURCES} SOURCES ${CPP_SOURCES}
) )

View File

@ -0,0 +1,163 @@
#include <QtNetwork>
#include <octave/oct.h>
#include "riSettings.h"
void getActiveCellCorners(NDArray& cellCornerValues, const QString &hostName, quint16 port, const qint32& caseId, const quint32& gridIndex, const QString& porosityModel)
{
QString serverName = hostName;
quint16 serverPort = port;
const int timeout = riOctavePlugin::timeOutMilliSecs;
QTcpSocket socket;
socket.connectToHost(serverName, serverPort);
if (!socket.waitForConnected(timeout))
{
error((("Connection: ") + socket.errorString()).toLatin1().data());
return;
}
// Create command and send it:
QString command = QString("GetActiveCellCorners %1 %2 %3").arg(caseId).arg(gridIndex).arg(porosityModel);
QByteArray cmdBytes = command.toLatin1();
QDataStream socketStream(&socket);
socketStream.setVersion(riOctavePlugin::qtDataStreamVersion);
socketStream << (qint64)(cmdBytes.size());
socket.write(cmdBytes);
// Get response. First wait for the header
while (socket.bytesAvailable() < (int)(2 * sizeof(quint64)))
{
if (!socket.waitForReadyRead(timeout))
{
error((("Wating for header: ") + socket.errorString()).toLatin1().data());
return;
}
}
// Read timestep count and blocksize
quint64 coordCount;
quint64 byteCount;
socketStream >> coordCount;
socketStream >> byteCount;
dim_vector dv (1, 1);
dv(0) = 3;
dv(1) = coordCount;
cellCornerValues.resize(dv);
if (!(byteCount && coordCount))
{
error ("Could not find the requested data in ResInsight");
return;
}
// Wait for available data for each column, then read data for each column
while (socket.bytesAvailable() < (qint64)(byteCount))
{
if (!socket.waitForReadyRead(timeout))
{
error((("Waiting for data: ") + socket.errorString()).toLatin1().data());
return;
}
OCTAVE_QUIT;
}
quint64 bytesRead = 0;
double* internalMatrixData = cellCornerValues.fortran_vec();
bytesRead = socket.read((char*)(internalMatrixData), byteCount);
if (byteCount != bytesRead)
{
error("Could not read binary double data properly from socket");
octave_stdout << "Active cell count: " << coordCount << std::endl;
}
return;
}
DEFUN_DLD (riGetActiveCellCorners, args, nargout,
"Usage:\n"
"\n"
" riGetActiveCellCorners([CaseId], GridIndex, [PorosityModel = “Matrix”|”Fracture”] )\n"
"\n"
"This function returns the UTM coordinates (X, Y, Z) of the 8 corners of each of the active cells.\n"
"If the CaseId is not defined, ResInsights Current Case is used.\n"
)
{
if (nargout < 1)
{
error("riGetActiveCellCorners: Missing output argument.\n");
print_usage();
return octave_value_list ();
}
int nargin = args.length ();
if (nargin < 1)
{
error("riGetActiveCellCorners: Too few arguments. The grid index argument is required.\n");
print_usage();
return octave_value_list ();
}
if (nargin > 3)
{
error("riGetActiveCellCorners: Too many arguments.\n");
print_usage();
return octave_value_list ();
}
qint32 caseId = -1;
quint32 gridIndex = 0;
std::string porosityModel = "Matrix";
if (nargin == 1)
{
gridIndex = args(0).uint_value();
}
else if (nargin == 2)
{
if (args(0).is_numeric_type() && args(1).is_numeric_type())
{
caseId = args(0).uint_value();
gridIndex = args(1).uint_value();
}
else
{
gridIndex = args(0).uint_value();
porosityModel = args(1).string_value();
}
}
else if (nargin == 3)
{
caseId = args(0).uint_value();
gridIndex = args(1).uint_value();
porosityModel = args(2).string_value();
}
if (porosityModel != "Matrix" && porosityModel != "Fracture")
{
error("riGetActiveCellProperty: The value for \"PorosityModel\" is unknown. Please use either \"Matrix\" or \"Fracture\"\n");
print_usage();
return octave_value_list ();
}
NDArray cellCornerValues;
getActiveCellCorners(cellCornerValues, "127.0.0.1", 40001, caseId, gridIndex, porosityModel.c_str());
return octave_value(cellCornerValues);
}

View File

@ -33,7 +33,7 @@ void getCellCenters(NDArray& cellCenterValues, const QString &hostName, quint16
// Get response. First wait for the header // Get response. First wait for the header
while (socket.bytesAvailable() < (int)(2 * sizeof(quint64))) while (socket.bytesAvailable() < (int)(5 * sizeof(quint64)))
{ {
if (!socket.waitForReadyRead(timeout)) if (!socket.waitForReadyRead(timeout))
{ {
@ -44,16 +44,28 @@ void getCellCenters(NDArray& cellCenterValues, const QString &hostName, quint16
// Read timestep count and blocksize // Read timestep count and blocksize
quint64 cellCountI;
quint64 cellCountJ;
quint64 cellCountK;
quint64 cellCount; quint64 cellCount;
quint64 byteCount; quint64 byteCount;
socketStream >> cellCount; socketStream >> cellCount;
socketStream >> cellCountI;
socketStream >> cellCountJ;
socketStream >> cellCountK;
socketStream >> byteCount; socketStream >> byteCount;
dim_vector dv (1, 1); // Create a 4D matrix, with the a column with the tree double value coords running as fastest index, then I, J, K
// Octave script to access coords
// coords = riGetCellCenters
// coords(:,i, j, k) # Will return the coords for given ijk location
dim_vector dv;
dv.resize(4);
dv(0) = 3; dv(0) = 3;
dv(1) = cellCount; dv(1) = cellCountI;
dv(2) = cellCountJ;
dv(3) = cellCountK;
cellCenterValues.resize(dv); cellCenterValues.resize(dv);
if (!(byteCount && cellCount)) if (!(byteCount && cellCount))
@ -73,8 +85,21 @@ void getCellCenters(NDArray& cellCenterValues, const QString &hostName, quint16
OCTAVE_QUIT; OCTAVE_QUIT;
} }
quint64 bytesRead = 0; octave_idx_type valueCount = cellCenterValues.length();
octave_stdout << " riGetCellCenters : I = " << cellCountI <<" J = " << cellCountJ << " K = " << cellCountK << std::endl;
octave_stdout << " riGetCellCenters : numDoubles = " << valueCount << std::endl;
double* internalMatrixData = cellCenterValues.fortran_vec(); double* internalMatrixData = cellCenterValues.fortran_vec();
#if 1
double val;
for (octave_idx_type i = 0; i < valueCount; i++)
{
socketStream >> internalMatrixData[i];
}
#else
quint64 bytesRead = 0;
bytesRead = socket.read((char*)(internalMatrixData), byteCount); bytesRead = socket.read((char*)(internalMatrixData), byteCount);
if (byteCount != bytesRead) if (byteCount != bytesRead)
@ -83,6 +108,8 @@ void getCellCenters(NDArray& cellCenterValues, const QString &hostName, quint16
octave_stdout << "Cell count: " << cellCount << std::endl; octave_stdout << "Cell count: " << cellCount << std::endl;
} }
#endif
return; return;
} }

View File

@ -0,0 +1,170 @@
#include <QtNetwork>
#include <octave/oct.h>
#include "riSettings.h"
void getCellCorners(NDArray& cellCornerValues, const QString &hostName, quint16 port, const qint32& caseId, const quint32& gridIndex)
{
QString serverName = hostName;
quint16 serverPort = port;
const int timeout = riOctavePlugin::timeOutMilliSecs;
QTcpSocket socket;
socket.connectToHost(serverName, serverPort);
if (!socket.waitForConnected(timeout))
{
error((("Connection: ") + socket.errorString()).toLatin1().data());
return;
}
// Create command and send it:
QString command = QString("GetCellCorners %1 %2").arg(caseId).arg(gridIndex);
QByteArray cmdBytes = command.toLatin1();
QDataStream socketStream(&socket);
socketStream.setVersion(riOctavePlugin::qtDataStreamVersion);
socketStream << (qint64)(cmdBytes.size());
socket.write(cmdBytes);
// Get response. First wait for the header
while (socket.bytesAvailable() < (int)(5 * sizeof(quint64)))
{
if (!socket.waitForReadyRead(timeout))
{
error((("Wating for header: ") + socket.errorString()).toLatin1().data());
return;
}
}
// Read timestep count and blocksize
quint64 cellCountI;
quint64 cellCountJ;
quint64 cellCountK;
quint64 cellCount;
quint64 byteCount;
socketStream >> cellCount;
socketStream >> cellCountI;
socketStream >> cellCountJ;
socketStream >> cellCountK;
socketStream >> byteCount;
// Create a 4D matrix, with the a column with the tree double value coords running as fastest index, then I, J, K
// Octave script to access coords
// coords = riGetCellCenters
// coords(:,i, j, k) # Will return the coords for given ijk location
dim_vector dv;
dv.resize(5);
dv(0) = 3;
dv(1) = 8;
dv(2) = cellCountI;
dv(3) = cellCountJ;
dv(4) = cellCountK;
cellCornerValues.resize(dv);
// octave_stdout << "GetCellCorners - coord count: " << coordCount << ", byteCount: " << byteCount << std::endl;
if (!(byteCount && cellCount))
{
error ("Could not find the requested data in ResInsight");
return;
}
octave_stdout << "GetCellCorners - before wait for data" << std::endl;
// Wait for available data for each column, then read data for each column
while (socket.bytesAvailable() < (qint64)(byteCount))
{
if (!socket.waitForReadyRead(timeout))
{
error((("Waiting for data: ") + socket.errorString()).toLatin1().data());
return;
}
OCTAVE_QUIT;
}
octave_idx_type valueCount = cellCornerValues.length();
octave_stdout << "GetCellCorners - after wait for data" << std::endl;
double* internalMatrixData = cellCornerValues.fortran_vec();
#if 1
double val;
for (octave_idx_type i = 0; i < valueCount; i++)
{
socketStream >> internalMatrixData[i];
}
#else
quint64 bytesRead = 0;
bytesRead = socket.read((char*)(internalMatrixData), byteCount);
if (byteCount != bytesRead)
{
error("Could not read binary double data properly from socket");
octave_stdout << "Cell count: " << cellCount << std::endl;
}
#endif
return;
}
DEFUN_DLD (riGetCellCorners, args, nargout,
"Usage:\n"
"\n"
" riGetCellCorners([CaseId], GridIndex )\n"
"\n"
"This function returns the UTM coordinates (X, Y, Z) of the 8 corners of all the cells in the grid.\n"
"If the CaseId is not defined, ResInsights Current Case is used.\n"
)
{
int nargin = args.length ();
if (nargin > 2)
{
error("riGetCellCorners: Too many arguments. CaseId is optional input argument.\n");
print_usage();
}
else if (nargout < 1)
{
error("riGetCellCorners: Missing output argument.\n");
print_usage();
}
else
{
NDArray cellCornerValues;
qint32 caseId = -1;
quint32 gridIndex = 0;
if (nargin == 1)
{
gridIndex = args(0).uint_value();
}
else if (nargin == 2)
{
unsigned int argCaseId = args(0).uint_value();
caseId = argCaseId;
gridIndex = args(1).uint_value();
}
getCellCorners(cellCornerValues, "127.0.0.1", 40001, caseId, gridIndex);
return octave_value(cellCornerValues);
}
return octave_value();
}