Add Support for Edge Conformal Corner-Point Grid Processing

CMakeLists_files.cmake

@@ -58,6 +58,7 @@ list (APPEND MAIN_SOURCE_FILES
   opm/grid/cpgpreprocess/facetopology.c
   opm/grid/cpgpreprocess/geometry.c
   opm/grid/cpgpreprocess/preprocess.c
+	opm/grid/cpgpreprocess/make_edge_conformal.cpp
   opm/grid/cpgpreprocess/uniquepoints.c
   opm/grid/UnstructuredGrid.c
   opm/grid/grid_equal.cpp
@@ -250,6 +251,7 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/grid/cpgpreprocess/facetopology.h
   opm/grid/cpgpreprocess/geometry.h
   opm/grid/cpgpreprocess/preprocess.h
+	opm/grid/cpgpreprocess/make_edge_conformal.hpp
   opm/grid/cpgpreprocess/uniquepoints.h
   opm/grid/transmissibility/trans_tpfa.h
   opm/grid/transmissibility/TransTpfa.hpp

opm/grid/CpGrid.hpp

@@ -246,8 +246,11 @@ namespace Dune
         std::vector<std::size_t>
         processEclipseFormat(const Opm::EclipseGrid* ecl_grid,
                              Opm::EclipseState* ecl_state,
-                             bool periodic_extension, bool turn_normals, bool clip_z,
-                             bool pinchActive);
+                             bool periodic_extension,
+			     bool turn_normals,
+			     bool clip_z,
+                             bool pinchActive,
+			     bool edge_conformal);
 
         /// Read the Eclipse grid format ('grdecl').
         ///
@@ -271,14 +274,14 @@ namespace Dune
         std::vector<std::size_t>
         processEclipseFormat(const Opm::EclipseGrid* ecl_grid,
                              Opm::EclipseState* ecl_state,
-                             bool periodic_extension, bool turn_normals = false, bool clip_z = false);
+                             bool periodic_extension, bool turn_normals = false, bool clip_z = false, bool edge_conformal =false);
 
 #endif
 
         /// Read the Eclipse grid format ('grdecl').
         /// \param input_data the data in grdecl format, declared in preprocess.h.
         /// \param remove_ij_boundary if true, will remove (i, j) boundaries. Used internally.
-        void processEclipseFormat(const grdecl& input_data, bool remove_ij_boundary, bool turn_normals = false);
+      void processEclipseFormat(const grdecl& input_data, bool remove_ij_boundary, bool turn_normals = false, bool edge_conformal = false);
 
         //@}
 

opm/grid/GridManager.cpp

@@ -40,18 +40,18 @@ namespace Opm
 
 #if HAVE_ECL_INPUT
     /// Construct a 3d corner-point grid from a deck.
-    GridManager::GridManager(const Opm::EclipseGrid& inputGrid)
+  GridManager::GridManager(const Opm::EclipseGrid& inputGrid,bool edge_conformal)
         : ug_(0)
     {
-        initFromEclipseGrid(inputGrid, std::vector<double>());
+      initFromEclipseGrid(inputGrid, std::vector<double>(), edge_conformal);
     }
 
 
     GridManager::GridManager(const Opm::EclipseGrid& inputGrid,
-                             const std::vector<double>& poreVolumes)
+                             const std::vector<double>& poreVolumes,bool edge_conformal)
         : ug_(0)
     {
-        initFromEclipseGrid(inputGrid, poreVolumes);
+      initFromEclipseGrid(inputGrid, poreVolumes, edge_conformal);
     }
 #endif
 
@@ -133,7 +133,9 @@ namespace Opm
 #if HAVE_ECL_INPUT
     // Construct corner-point grid from EclipseGrid.
     void GridManager::initFromEclipseGrid(const Opm::EclipseGrid& inputGrid,
-                                          const std::vector<double>& poreVolumes)
+                                          const std::vector<double>& poreVolumes,
+					  bool edge_conformal
+					  )
     {
         struct grdecl g;
 
@@ -166,7 +168,7 @@ namespace Opm
                        actnum, opmfil, zcorn.data());
         }
 
-        ug_ = create_grid_cornerpoint(&g, z_tolerance);
+        ug_ = create_grid_cornerpoint(&g, z_tolerance, edge_conformal);
         if (!ug_) {
             OPM_THROW(std::runtime_error, "Failed to construct grid.");
         }

opm/grid/GridManager.hpp

@@ -47,15 +47,15 @@ class EclipseGrid;
 
 #if HAVE_ECL_INPUT
         /// Construct a grid from an EclipseState::EclipseGrid instance.
-        explicit GridManager(const EclipseGrid& inputGrid);
+      explicit GridManager(const EclipseGrid& inputGrid, bool edge_conformal = false);
 
         /// Construct a grid from an EclipseState::EclipseGrid instance,
         /// giving an explicit set of pore volumes to be used for MINPV
         /// considerations.
         /// \input[in] eclipseGrid    encapsulates a corner-point grid given from a deck
         /// \input[in] poreVolumes    one element per logical cartesian grid element
         GridManager(const EclipseGrid& inputGrid,
-                    const std::vector<double>& poreVolumes);
+                    const std::vector<double>& poreVolumes, bool edge_conformal);
 #endif
 
         /// Construct a 2d cartesian grid with cells of unit size.
@@ -92,7 +92,7 @@ class EclipseGrid;
 #if HAVE_ECL_INPUT
         // Construct corner-point grid from EclipseGrid.
         void initFromEclipseGrid(const EclipseGrid& inputGrid,
-                                 const std::vector<double>& poreVolumes);
+                                 const std::vector<double>& poreVolumes,bool edge_conformal);
 #endif
 
         // The managed UnstructuredGrid.

opm/grid/common/GridPartitioning.cpp

@@ -309,7 +309,7 @@ void addOverlapCornerCell(const CpGrid& grid, int owner,
             {
                 // Note: multiple adds for same process are possible
                 exportList.emplace_back(nb_index, owner, AttributeSet::copy);
-                exportList.emplace_back(my_index, cell_part[nb_index],  AttributeSet::copy);
+                //exportList.emplace_back(my_index, cell_part[nb_index],  AttributeSet::copy);
                 return;
             }
         }
@@ -440,7 +440,7 @@ void addOverlapLayer(const CpGrid& grid, int index, const CpGrid::Codim<0>::Enti
                     {
                         // Note: multiple adds for same process are possible
                         exportList.emplace_back(nb_index, owner, AttributeSet::copy);
-                        exportList.emplace_back(index, cell_part[nb_index],  AttributeSet::copy);
+                        //exportList.emplace_back(index, cell_part[nb_index],  AttributeSet::copy);
                         if ( recursion_deps>0 )
                         {
                             // Add another layer

opm/grid/cornerpoint_grid.c

@@ -18,6 +18,7 @@
 #include <opm/grid/cornerpoint_grid.h>
 #include <opm/grid/cpgpreprocess/geometry.h>
 #include <opm/grid/cpgpreprocess/preprocess.h>
+#include <opm/grid/cpgpreprocess/make_edge_conformal.hpp>
 #include <opm/grid/UnstructuredGrid.h>
 
 
@@ -161,7 +162,7 @@ void compute_geometry(struct UnstructuredGrid *g)
 
 
 struct UnstructuredGrid *
-create_grid_cornerpoint(const struct grdecl *in, double tol)
+create_grid_cornerpoint(const struct grdecl *in, double tol, bool edge_conformal)
 {
     struct UnstructuredGrid *g;
    int                      ok;
@@ -173,10 +174,16 @@ create_grid_cornerpoint(const struct grdecl *in, double tol)
        return NULL;
    }
 
-   ok = process_grdecl(in, tol, NULL, &pg, false);
+   ok = process_grdecl(in, tol, NULL, &pg, false, edge_conformal);
+   if(edge_conformal){
+     // if add cells is done one could skip som of the code later
+     add_cells(&pg);
+     make_edge_conformal(&pg);
+   }
+
    if (!ok)
    {
-       free_processed_grid(&pg);
+     free_processed_grid(&pg, edge_conformal);
        destroy_grid(g);
        return NULL;
    }
@@ -235,7 +242,7 @@ create_grid_cornerpoint(const struct grdecl *in, double tol)
        pg.local_cell_index = NULL;
    }
 
-   free_processed_grid(&pg);
+   free_processed_grid(&pg, edge_conformal);
 
    return g;
 }

opm/grid/cornerpoint_grid.h

@@ -68,7 +68,7 @@ extern "C" {
      * destroy_grid().
      */
     struct UnstructuredGrid *
-    create_grid_cornerpoint(const struct grdecl *in, double tol);
+    create_grid_cornerpoint(const struct grdecl *in, double tol, bool edge_conformal);
 
 
     /**

opm/grid/cpgpreprocess/facetopology.c

@@ -64,11 +64,14 @@
 static int *
 computeFaceTopology(const int *a1, const int *a2,
                     const int *b1, const int *b2,
-                    int intersect[4], int *faces)
+                    int intersect[4], int *faces,
+		    bool edge_conformal)
 {
     int mask[8];
     int k;
     int *f;
+    int zn;
+    int i;
 
     for (k = 0; k < 8; k++) { mask[k] = -1; }
 
@@ -139,12 +142,74 @@ computeFaceTopology(const int *a1, const int *a2,
         }
     }
     f = faces;
+    if(!edge_conformal){
     for (k=7; k>=0; --k){
         if(mask[k] != -1){
             *f++ = mask[k];
         }
     }
+    }else{
+        k=7;
+        if(mask[k] != -1){
+            *f++ = mask[k];
+        }
+        /* intersection top top*/
+	k=6;
+        if(mask[k] != -1){
+            *f++ = mask[k];
+        }
+        /*upper right*/
+	k=5;
+        if(mask[k] != -1){
+            *f++ = mask[k];
+        }
+        /* intersection bottom-top or top bottum */
+	k=4;
+        if(mask[k] != -1){
+            *f++ = mask[k];
+        }
+        /* add possible hanging nodes */
+
+        if((mask[2] != -1) && (mask[4] != -1)){
+            assert(mask[3] == -1);
+            if((mask[2]-mask[4])> 1){
+                for(zn = mask[4]+1; zn<mask[2]; ++zn){
+                    *f++ = zn;
+                }
+            }
+        }
+
+        /* lower right */
+	k=3;
+        if(mask[k] != -1){
+            *f++ = mask[k];
+        }
+	k=2;
+        if(mask[k] != -1){
+            *f++ = mask[k];
+        }
+	k=1;
+        if(mask[k] != -1){
+            *f++ = mask[k];
+        }	
+	k=0;
+        if(mask[k] != -1){
+            *f++ = mask[k];
+        }	
+        if((mask[6] != -1) && (mask[0] != -1)){
+            assert(mask[7] == -1);
+            if((mask[0]-mask[6])> 1){
+		fprintf(stdout, "mask 6 %d \n", mask[6]);
+                for(zn = mask[0]-1; zn>mask[6]; --zn){
+                    *f++ = zn;
+		    fprintf(stderr, "adde mpde %d \n", zn);
+                }
+		fprintf(stderr, "mask 0 %d \n", mask[0]);
+            }
+        }
+    }
 
+
 #if DEBUG>1
     /* Check for repeated nodes:*/
     int i;
@@ -206,7 +271,8 @@ faceintersection(const int *a1, const int *a2,
 void findconnections(int n, int *pts[4],
                      int *intersectionlist,
                      int *work,
-                     struct processed_grid *out)
+                     struct processed_grid *out,
+		     bool edge_conformal)
 {
     /* vectors of point numbers for faces a(b) on pillar 1(2) */
     int *a1 = pts[0];
@@ -226,6 +292,7 @@ void findconnections(int n, int *pts[4],
     int i,j=0;
     int intersect[4];
     int *tmp;
+    int zn;
     /* for (i=0; i<2*n; work[i++]=-1); */
 
     for (i = 0; i < 4; i++) { intersect[i] = -1; }
@@ -332,7 +399,7 @@ void findconnections(int n, int *pts[4],
                             *c++ = cell_a;
                             *c++ = cell_b;
 
-                            f = computeFaceTopology(a1+i, a2+i, b1+j, b2+j, intersect, f);
+                            f = computeFaceTopology(a1+i, a2+i, b1+j, b2+j, intersect, f, edge_conformal);
 
                             out->face_ptr[++out->number_of_faces] = f - out->face_nodes;
                         }

opm/grid/cpgpreprocess/facetopology.h

@@ -39,7 +39,8 @@
 void findconnections(int n, int *pts[4],
                      int *intersectionlist,
                      int *work,
-                     struct processed_grid *out);
+                     struct processed_grid *out,
+		     bool edge_conformal);
 
 #endif /* OPM_FACETOPOLOGY_HEADER */