#4259 Implement mapaxes keyword with possibility to export in local coordinates

ApplicationCode/Commands/ExportCommands/RicExportEclipseSectorModelFeature.cpp

@@ -100,6 +100,7 @@ void RicExportEclipseSectorModelFeature::executeCommand(RimEclipseView* view,
 
         bool worked = RifEclipseInputFileTools::exportGrid(exportSettings.exportGridFilename(),
                                                            view->eclipseCase()->eclipseCaseData(),
+                                                           exportSettings.exportInLocalCoordinates(),
                                                            cellVisibilityForActnum,
                                                            min,
                                                            max,

ApplicationCode/Commands/ExportCommands/RicExportEclipseSectorModelUi.cpp

@@ -64,7 +64,7 @@ RicExportEclipseSectorModelUi::RicExportEclipseSectorModelUi(RigEclipseCaseData*
     CAF_PDM_InitField(&exportGrid, "ExportGrid", true, "Export Grid", "", "Includes COORD, ZCORN and ACTNUM", "");
     CAF_PDM_InitField(&exportGridFilename, "ExportGridFilename", QString(), "Grid File Name", "", "", "");
     exportGridFilename.uiCapability()->setUiEditorTypeName(caf::PdmUiFilePathEditor::uiEditorTypeName());
-
+    CAF_PDM_InitField(&exportInLocalCoordinates, "ExportInLocalCoords", false, "Export in Local Coordinates", "", "Remove UTM location on export", "");
     CAF_PDM_InitField(&makeInvisibleCellsInactive, "InvisibleCellActnum", false, "Make Invisible Cells Inactive", "", "", "");
 
     CAF_PDM_InitFieldNoDefault(&exportFaults, "ExportFaults", "Export Faults", "", "", "");
@@ -156,6 +156,7 @@ void RicExportEclipseSectorModelUi::defineUiOrdering(QString uiConfigName, caf::
     caf::PdmUiGroup* gridGroup = uiOrdering.addNewGroup("Grid Export");
     gridGroup->add(&exportGrid);
     gridGroup->add(&exportGridFilename);
+    gridGroup->add(&exportInLocalCoordinates);
     exportGridFilename.uiCapability()->setUiReadOnly(!exportGrid());        
     gridGroup->add(&makeInvisibleCellsInactive);
 

ApplicationCode/Commands/ExportCommands/RicExportEclipseSectorModelUi.h

@@ -55,7 +55,7 @@ public:
 
     caf::PdmField<bool>                    exportGrid;
     caf::PdmField<QString>                 exportGridFilename;
-
+    caf::PdmField<bool>                    exportInLocalCoordinates;
     caf::PdmField<bool>                    makeInvisibleCellsInactive;
 
     caf::PdmField<ResultExportOptionsEnum> exportFaults;

ApplicationCode/FileInterface/RifEclipseInputFileTools.cpp

@@ -47,6 +47,7 @@
 
 #include "ert/ecl/ecl_box.hpp"
 #include "ert/ecl/ecl_kw.h"
+#include "ert/ecl/ecl_grid.hpp"
 
 QString includeKeyword("INCLUDE");
 QString faultsKeyword("FAULTS");
@@ -212,6 +213,16 @@ bool RifEclipseInputFileTools::openGridFile(const QString&      fileName,
         mainGrid->setFaults(faults);
     }
 
+    bool useMapAxes = ecl_grid_use_mapaxes(inputGrid);
+    eclipseCase->mainGrid()->setUseMapAxes(useMapAxes);
+
+    if (useMapAxes)
+    {
+        std::array<double, 6> mapAxesValues;
+        ecl_grid_init_mapaxes_data_double(inputGrid, mapAxesValues.data());
+        eclipseCase->mainGrid()->setMapAxes(mapAxesValues);
+    }    
+
     progress.setProgress(8);
     progress.setProgressDescription("Cleaning up ...");
 
@@ -233,6 +244,7 @@ bool RifEclipseInputFileTools::openGridFile(const QString&      fileName,
 //--------------------------------------------------------------------------------------------------
 bool RifEclipseInputFileTools::exportGrid(const QString&         fileName,
                                           RigEclipseCaseData*    eclipseCase,
+                                          bool                   exportInLocalCoordinates,
                                           const cvf::UByteArray* cellVisibilityOverrideForActnum,
                                           const cvf::Vec3st&     min,
                                           const cvf::Vec3st&     maxIn,
@@ -271,6 +283,20 @@ bool RifEclipseInputFileTools::exportGrid(const QString&         fileName,
     std::vector<int*> ecl_coords;
     ecl_coords.reserve(mainGrid->cellCount() * cellsPerOriginal);
 
+    std::array<float, 6> mapAxes = mainGrid->mapAxesF();
+    cvf::Mat4d mapAxisTrans = mainGrid->mapAxisTransform();
+    if (exportInLocalCoordinates)
+    {
+        cvf::Vec3d minPoint3d(mainGrid->boundingBox().min());
+        cvf::Vec2f minPoint2f(minPoint3d.x(), minPoint3d.y());
+        cvf::Vec2f origin(mapAxes[2] - minPoint2f.x(), mapAxes[3] - minPoint2f.y());
+        cvf::Vec2f xPoint = cvf::Vec2f(mapAxes[4], mapAxes[5]) - minPoint2f;
+        cvf::Vec2f yPoint = cvf::Vec2f(mapAxes[0], mapAxes[1]) - minPoint2f;
+        mapAxes = { yPoint.x(), yPoint.y(), origin.x(), origin.y(), xPoint.x(), xPoint.y() };
+
+        mapAxisTrans.setTranslation(mapAxisTrans.translation() - minPoint3d);
+    }
+
     const size_t* cellMappingECLRi = RifReaderEclipseOutput::eclipseCellIndexMapping();
 
     int incrementalIndex = 0;
@@ -306,6 +332,14 @@ bool RifEclipseInputFileTools::exportGrid(const QString&         fileName,
                 std::array<cvf::Vec3d, 8> cellCorners;
                 mainGrid->cellCornerVertices(mainIndex, cellCorners.data());
 
+                if (mainGrid->useMapAxes())
+                {
+                    for (cvf::Vec3d& corner : cellCorners)
+                    {
+                        corner.transformPoint(mapAxisTrans);
+                    }
+                }
+
                 auto refinedCoords =
                     RiaCellDividingTools::createHexCornerCoords(cellCorners, refinement.x(), refinement.y(), refinement.z());
 
@@ -332,7 +366,7 @@ bool RifEclipseInputFileTools::exportGrid(const QString&         fileName,
     }
 
     ecl_grid_type* mainEclGrid =
-        ecl_grid_alloc_GRID_data((int)ecl_coords.size(), ecl_nx, ecl_ny, ecl_nz, 5, &ecl_coords[0], &ecl_corners[0], false, NULL);
+        ecl_grid_alloc_GRID_data((int)ecl_coords.size(), ecl_nx, ecl_ny, ecl_nz, 5, &ecl_coords[0], &ecl_corners[0], false, mapAxes.data());
     progress.setProgress(mainGrid->cellCount());
 
     for (float* floatArray : ecl_corners)
@@ -399,7 +433,12 @@ bool RifEclipseInputFileTools::exportKeywords(const QString&              result
         std::vector<double> resultValues;
 
         RigEclipseResultAddress resAddr(RiaDefines::STATIC_NATIVE, keyword);
+        CVF_ASSERT(cellResultsData->hasResultEntry(resAddr));
+
         cellResultsData->ensureKnownResultLoaded(resAddr);
+        
+        CVF_ASSERT(!cellResultsData->cellScalarResults(resAddr).empty());
+
         resultValues = cellResultsData->cellScalarResults(resAddr)[0];
         CVF_ASSERT(!resultValues.empty());
 

ApplicationCode/FileInterface/RifEclipseInputFileTools.h

@@ -63,6 +63,7 @@ public:
 
     static bool exportGrid(const QString&         gridFileName,
                            RigEclipseCaseData*    eclipseCase,
+                           bool                   exportInLocalCoordinates,
                            const cvf::UByteArray* cellVisibilityOverrideForActnum = nullptr,
                            const cvf::Vec3st&     min                             = cvf::Vec3st::ZERO,
                            const cvf::Vec3st&     max                             = cvf::Vec3st::UNDEFINED,

ApplicationCode/ReservoirDataModel/RigMainGrid.cpp

@@ -39,6 +39,9 @@ RigMainGrid::RigMainGrid()
 
     m_flipXAxis = false;
     m_flipYAxis = false;
+
+    m_useMapAxes = false;
+    m_mapAxes    = defaultMapAxes();
 }
 
 RigMainGrid::~RigMainGrid() {}
@@ -787,3 +790,84 @@ const std::string& RigMainGrid::associatedWellPathName() const
     static const std::string EMPTY_STRING;
     return EMPTY_STRING;
 }
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+void RigMainGrid::setUseMapAxes(bool useMapAxes)
+{
+    m_useMapAxes = useMapAxes;
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+bool RigMainGrid::useMapAxes() const
+{
+    return m_useMapAxes;
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+void RigMainGrid::setMapAxes(const std::array<double, 6>& mapAxes)
+{
+    m_mapAxes = mapAxes;
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+const std::array<double, 6>& RigMainGrid::mapAxes() const
+{
+    return m_mapAxes;
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+std::array<float, 6> RigMainGrid::mapAxesF() const
+{
+    std::array<float, 6> floatAxes;
+    for (size_t i = 0; i < 6; ++i)
+    {
+        floatAxes[i] = (float)m_mapAxes[i];
+    }
+    return floatAxes;
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+cvf::Mat4d RigMainGrid::mapAxisTransform() const
+{
+    cvf::Mat4d mapAxisTrans;
+    if (m_useMapAxes)
+    {
+        cvf::Vec3d origin(m_mapAxes[2], m_mapAxes[3], 0.0);
+        cvf::Vec3d xAxis = cvf::Vec3d(m_mapAxes[4] - origin[0], m_mapAxes[5] - origin[1], 0.0).getNormalized();
+        cvf::Vec3d yAxis = cvf::Vec3d(m_mapAxes[0] - origin[0], m_mapAxes[1] - origin[1], 0.0).getNormalized();
+        cvf::Vec3d zAxis(0.0, 0.0, 1.0);
+        mapAxisTrans = cvf::Mat4d::fromCoordSystemAxes(&xAxis, &yAxis, &zAxis);
+        mapAxisTrans.setTranslation(origin);
+        mapAxisTrans.invert();
+    }
+    else
+    {
+        mapAxisTrans.setIdentity();
+    }
+    return mapAxisTrans;
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+std::array<double, 6> RigMainGrid::defaultMapAxes()
+{
+    const double origin[2] = {0.0, 0.0};
+    const double xPoint[2] = {1.0, 0.0};
+    const double yPoint[2] = {0.0, 1.0};
+
+    // Order (see Elipse Reference Manual for keyword MAPAXES): Y_x, Y_y, O_x, O_y, X_x, X_y
+    return { yPoint[0], yPoint[1], origin[0], origin[1], xPoint[0], xPoint[1] };
+}

ApplicationCode/ReservoirDataModel/RigMainGrid.h

@@ -94,11 +94,19 @@ public:
     bool                                    isTempGrid() const override;
     const std::string&                      associatedWellPathName() const override;
 
+    void                                    setUseMapAxes(bool useMapAxes);
+    bool                                    useMapAxes() const;
+    void                                    setMapAxes(const std::array<double, 6>& mapAxes);
+    const std::array<double, 6>&            mapAxes() const;
+    std::array<float, 6>                    mapAxesF() const;
+
+    cvf::Mat4d                              mapAxisTransform() const;
 private:
     void                                    initAllSubCellsMainGridCellIndex();
     void                                    buildCellSearchTree();
     bool                                    hasFaultWithName(const QString& name) const;
 
+    static std::array<double, 6>            defaultMapAxes();
 private:
     std::vector<cvf::Vec3d>                 m_nodes;        ///< Global vertex table
     std::vector<RigCell>                    m_cells;        ///< Global array of all cells in the reservoir (including the ones in LGR's)
@@ -116,5 +124,8 @@ private:
 
     bool                                    m_flipXAxis;
     bool                                    m_flipYAxis;
+
+    bool                                    m_useMapAxes;
+    std::array<double, 6>                   m_mapAxes;
 };