Support using DXV, DYV, DZV in deck to create grid.

opm/core/GridManager.cpp

@@ -22,6 +22,8 @@
 #include <opm/core/grid.h>
 #include <opm/core/grid/cart_grid.h>
 #include <opm/core/grid/cornerpoint_grid.h>
+#include <algorithm>
+#include <numeric>
 
 
 
@@ -33,36 +35,25 @@ namespace Opm
     /// Construct a 3d corner-point grid from a deck.
     GridManager::GridManager(const Opm::EclipseGridParser& deck)
     {
-	// Extract data from deck.
-	const std::vector<double>& zcorn = deck.getFloatingPointValue("ZCORN");
-	const std::vector<double>& coord = deck.getFloatingPointValue("COORD");
-        const int* actnum = 0;
-        if (deck.hasField("ACTNUM")) {
-            actnum = &(deck.getIntegerValue("ACTNUM")[0]);
+        // We accept two different ways to specify the grid.
+        //    1. Corner point format.
+        //       Requires ZCORN, COORDS, DIMENS or SPECGRID, optionally ACTNUM.
+        //       For this format, we will verify that DXV, DYV, DZV,
+        //       DEPTHZ and TOPS are not present.
+        //    2. Tensor grid format.
+        //       Requires DXV, DYV, DZV, optionally DEPTHZ or TOPS.
+        //       For this format, we will verify that ZCORN, COORDS
+        //       and ACTNUM are not present.
+        //       Note that for TOPS, we only allow a uniform vector of values.
+
+        if (deck.hasField("ZCORN") && deck.hasField("COORD")) {
+            initFromDeckCornerpoint(deck);
+        } else if (deck.hasField("DXV") && deck.hasField("DYV") && deck.hasField("DZV")) {
+            initFromDeckTensorgrid(deck);
+        } else {
+            THROW("Could not initialize grid from deck. "
+                  "Need either ZCORN + COORD or DXV + DYV + DZV keywords.");
         }
-	std::vector<int> dims;
-	if (deck.hasField("DIMENS")) {
-	    dims = deck.getIntegerValue("DIMENS");
-	} else if (deck.hasField("SPECGRID")) {
-	    dims = deck.getSPECGRID().dimensions;
-	} else {
-	    THROW("Deck must have either DIMENS or SPECGRID.");
-	}
-
-	// Collect in input struct for preprocessing.
-	struct grdecl grdecl;
-	grdecl.zcorn = &zcorn[0];
-	grdecl.coord = &coord[0];
-	grdecl.actnum = actnum;
-	grdecl.dims[0] = dims[0];
-	grdecl.dims[1] = dims[1];
-	grdecl.dims[2] = dims[2];
-
-	// Process grid.
-	ug_ = create_grid_cornerpoint(&grdecl, 0.0);
-	if (!ug_) {
-	    THROW("Failed to construct grid.");
-	}
     }
 
 
@@ -124,5 +115,94 @@ namespace Opm
 
 
 
+    // Construct corner-point grid from deck.
+    void GridManager::initFromDeckCornerpoint(const Opm::EclipseGridParser& deck)
+    {
+	// Extract data from deck.
+	const std::vector<double>& zcorn = deck.getFloatingPointValue("ZCORN");
+	const std::vector<double>& coord = deck.getFloatingPointValue("COORD");
+        const int* actnum = 0;
+        if (deck.hasField("ACTNUM")) {
+            actnum = &(deck.getIntegerValue("ACTNUM")[0]);
+        }
+	std::vector<int> dims;
+	if (deck.hasField("DIMENS")) {
+	    dims = deck.getIntegerValue("DIMENS");
+	} else if (deck.hasField("SPECGRID")) {
+	    dims = deck.getSPECGRID().dimensions;
+	} else {
+	    THROW("Deck must have either DIMENS or SPECGRID.");
+	}
+
+	// Collect in input struct for preprocessing.
+	struct grdecl grdecl;
+	grdecl.zcorn = &zcorn[0];
+	grdecl.coord = &coord[0];
+	grdecl.actnum = actnum;
+	grdecl.dims[0] = dims[0];
+	grdecl.dims[1] = dims[1];
+	grdecl.dims[2] = dims[2];
+
+	// Process grid.
+	ug_ = create_grid_cornerpoint(&grdecl, 0.0);
+	if (!ug_) {
+	    THROW("Failed to construct grid.");
+	}
+    }
+
+
+    namespace
+    {
+        std::vector<double> coordsFromDeltas(const std::vector<double> deltas)
+        {
+            std::vector<double> coords(deltas.size() + 1);
+            coords[0] = 0.0;
+            std::partial_sum(deltas.begin(), deltas.end(), coords.begin() + 1);
+            return coords;
+        }
+    } // anonymous namespace
+
+
+    // Construct tensor grid from deck.
+    void GridManager::initFromDeckTensorgrid(const Opm::EclipseGridParser& deck)
+    {
+        // Extract coordinates (or offsets from top, in case of z).
+	const std::vector<double>& dxv = deck.getFloatingPointValue("DXV");
+	const std::vector<double>& dyv = deck.getFloatingPointValue("DYV");
+	const std::vector<double>& dzv = deck.getFloatingPointValue("DZV");
+        std::vector<double> x = coordsFromDeltas(dxv);
+        std::vector<double> y = coordsFromDeltas(dyv);
+        std::vector<double> z = coordsFromDeltas(dzv);
+
+        // Extract top corner depths, if available.
+        const double* top_depths = 0;
+        std::vector<double> top_depths_vec;
+        if (deck.hasField("DEPTHZ")) {
+            const std::vector<double>& depthz = deck.getFloatingPointValue("DEPTHZ");
+            if (depthz.size() != x.size()*y.size()) {
+                THROW("Incorrect size of DEPTHZ: " << depthz.size());
+            }
+            top_depths = &depthz[0];
+        } else if (deck.hasField("TOPS")) {
+            // We only support constant values for TOPS.
+            // It is not 100% clear how we best can deal with
+            // varying TOPS (stair-stepping grid, or not).
+            const std::vector<double>& tops = deck.getFloatingPointValue("TOPS");
+            if (std::count(tops.begin(), tops.end(), tops[0]) != int(tops.size())) {
+                THROW("We do not support nonuniform TOPS, please use ZCORN/COORDS instead.");
+            }
+            top_depths_vec.resize(x.size()*y.size(), tops[0]);
+            top_depths = &top_depths_vec[0];
+        }
+
+        // Construct grid.
+        ug_ = create_grid_tensor3d(dxv.size(), dyv.size(), dzv.size(),
+                                   &x[0], &y[0], &z[0], top_depths);
+	if (!ug_) {
+	    THROW("Failed to construct grid.");
+	}
+    }
+
+
 
 } // namespace Opm

opm/core/GridManager.hpp

@@ -33,13 +33,14 @@ namespace Opm
     /// encapsulates creation and destruction of the grid.
     /// The following grid types can be constructed:
     ///   - 3d corner-point grids (from deck input)
+    ///   - 3d tensor grids (from deck input)
     ///   - 2d cartesian grids
     ///   - 3d cartesian grids
     /// The resulting UnstructuredGrid is available through the c_grid() method.
     class GridManager
     {
     public:
-	/// Construct a 3d corner-point grid from a deck.
+	/// Construct a 3d corner-point grid or tensor grid from a deck.
 	GridManager(const Opm::EclipseGridParser& deck);
 
 	/// Construct a 2d cartesian grid with cells of unit size.
@@ -64,6 +65,12 @@ namespace Opm
 	// Disable copying and assignment.
 	GridManager(const GridManager& other);
 	GridManager& operator=(const GridManager& other);
+
+        // Construct corner-point grid from deck.
+        void initFromDeckCornerpoint(const Opm::EclipseGridParser& deck);
+        // Construct tensor grid from deck.
+        void initFromDeckTensorgrid(const Opm::EclipseGridParser& deck);
+
 	// The managed UnstructuredGrid.
 	UnstructuredGrid* ug_;
     };