Merge pull request #762 from aritorto/globalCellLgr

Level  global cell for CpGrid with LGRs

opm/grid/CpGrid.hpp

@@ -889,10 +889,9 @@ namespace Dune
                                    const std::vector<std::array<int,3>>& cells_per_dim_vec,
                                    const int& preAdaptMaxLevel) const;
 
-        /// @brief Define the cells, cell_to_point_, global_cell_, cell_to_face_, face_to_cell_, for the leaf grid view (or adapted grid).
+        /// @brief Define the cells, cell_to_point_, cell_to_face_, face_to_cell_, for the leaf grid view (or adapted grid).
         void populateLeafGridCells(Dune::cpgrid::EntityVariableBase<cpgrid::Geometry<3,3>>& adapted_cells,
                                    std::vector<std::array<int,8>>& adapted_cell_to_point,
-                                   std::vector<int>& adapted_global_cell,
                                    const int& cell_count,
                                    cpgrid::OrientedEntityTable<0,1>& adapted_cell_to_face,
                                    cpgrid::OrientedEntityTable<1,0>& adapted_face_to_cell,
@@ -922,7 +921,6 @@ namespace Dune
                                            /* Leaf grid View Cells argumemts  */
                                            Dune::cpgrid::EntityVariableBase<cpgrid::Geometry<3,3>>& adapted_cells,
                                            std::vector<std::array<int,8>>& adapted_cell_to_point,
-                                           std::vector<int>& adapted_global_cell,
                                            const int& cell_count,
                                            cpgrid::OrientedEntityTable<0,1>& adapted_cell_to_face,
                                            cpgrid::OrientedEntityTable<1,0>& adapted_face_to_cell,
@@ -950,6 +948,26 @@ namespace Dune
                                        const int& preAdaptMaxLevel,
                                        const int& newLevels);
 
+
+        /// @brief For refined level grids created based on startIJK and endIJK values, compute the "local ijk/Cartesian index" within the LGR.
+        ///
+        /// It's confusing that this "localIJK" is stored in CpGridData member global_cell_. Potential explanation: level zero grid is also called
+        /// "GLOBAL" grid.
+        /// Example: a level zero grid with dimension 4x3x3, an LGR with startIJK = {1,2,2}, endIJK = {3,3,3}, and cells_per_dim = {2,2,2}.
+        /// Then the dimension of the LGR is (3-1)*2 x (3-2)*2 x (3-2)* 2 = 4x2x2 = 16. Therefore the global_cell_lgr minimim value should be 0,
+        /// the maximum should be 15.
+        /// To invoke this method, each refined level grid must have 1. logical_cartesian_size_, 2. cell_to_idxInParentCell_, and 3. cells_per_dim_
+        /// already populated.
+        ///
+        /// @param [in] level    Grid index where LGR is stored
+        /// @param [out] global_cell_lgr
+        void computeGlobalCellLgr(const int& level, const std::array<int,3>& startIJK, std::vector<int>& global_cell_lgr);
+
+        /// @brief For a leaf grid with with LGRs, we assign the global_cell_ values of either the parent cell or the equivalent cell from
+        ///        level zero.
+        ///        For nested refinement, we lookup the oldest ancestor, from level zero.
+        void computeGlobalCellLeafGridViewWithLgrs(std::vector<int>& global_cell_leaf);
+
         /// @brief Get the ijk index of a refined corner, given its corner index of a single-cell-refinement.
         ///
         /// Given a single-cell, we refine it in {nx, ny, nz} refined children cells (per direction). Then, this single-cell-refinement

opm/grid/cpgrid/CartesianIndexMapper.hpp

@@ -52,7 +52,6 @@ namespace Dune
 
         int compressedLevelZeroSize() const
         {
-
             return (*grid_.currentData()[0]).size(0);
         }
 

opm/grid/cpgrid/CpGrid.cpp

@@ -633,6 +633,52 @@ const std::vector<int>& CpGrid::globalCell() const
     return currentData().back() -> global_cell_;
 }
 
+void CpGrid::computeGlobalCellLgr(const int& level, const std::array<int,3>& startIJK, std::vector<int>& global_cell_lgr)
+{
+    assert(level);
+    for (const auto& element : elements(levelGridView(level))) {
+        // Element belogns to an LGR, therefore has a father. Get IJK of the father in the level grid the father was born.
+        // For CARFIN, parent cells belong to level 0.
+        std::array<int,3> parentIJK = {0,0,0};
+        currentData()[element.father().level()]->getIJK(element.father().index(), parentIJK);
+        // Each parent cell has been refined in cells_per_dim[0]*cells_per_dim[1]*cells_per_dim[2] child cells.
+        // element has certain 'position' inside its parent cell that can be described with 'IJK' indices, let's denote them by ijk,
+        // where 0<= i < cells_per_dim[0], 0<= j < cells_per_dim[1], 0<= k < cells_per_dim[2].
+        const auto& cells_per_dim = currentData()[level]->cells_per_dim_;
+        //
+        // Refined cell (here 'element') has "index in parent cell": k*cells_per_dim[0]*cells_per_dim[1] + j*cells_per_dim[0] + i
+        // and it's stored in  cell_to_idxInParentCell_.
+        auto idx_in_parent_cell =  currentData()[level]-> cell_to_idxInParentCell_[element.index()];
+        // Find ijk.
+        std::array<int,3> childIJK = currentData()[level]-> getIJK(idx_in_parent_cell, cells_per_dim);
+        // The corresponding lgrIJK can be computed as follows:
+        const std::array<int,3>& lgrIJK = { ( (parentIJK[0] - startIJK[0])*cells_per_dim[0] ) + childIJK[0],  // Shift parent index according to the startIJK of the LGR.
+                                            ( (parentIJK[1] - startIJK[1])*cells_per_dim[1] ) + childIJK[1],
+                                            ( (parentIJK[2] - startIJK[2])*cells_per_dim[2] ) + childIJK[2] };
+        // Dimensions of the "patch of cells" formed when providing startIJK and endIJK for an LGR
+        const auto& lgr_logical_cartesian_size = currentData()[level]->logical_cartesian_size_;
+        global_cell_lgr[element.index()] = (lgrIJK[2]*lgr_logical_cartesian_size[0]*lgr_logical_cartesian_size[1]) + (lgrIJK[1]*lgr_logical_cartesian_size[0]) + lgrIJK[0];
+    }
+}
+
+void CpGrid::computeGlobalCellLeafGridViewWithLgrs(std::vector<int>& global_cell_leaf)
+{
+    for (const auto& element: elements(leafGridView()))
+    {
+        // When refine via CpGrid::addLgrsUpdateGridView(/*...*/), level-grid to lookup global_cell_ is equal to level-zero-grid
+        // In the context of allowed nested refinement, we lookup for the oldest ancestor, also belonging to level-zero-grid.
+        auto ancestor = element.getOrigin();
+        int origin_in_level_zero = ancestor.index();
+        while (ancestor.level()>0){
+            ancestor = ancestor.getOrigin();
+            origin_in_level_zero = ancestor.getOrigin().index();
+        }
+        assert(ancestor.level()==0);
+        assert(origin_in_level_zero < currentData().front()->size(0));
+        global_cell_leaf[element.index()] = currentData().front()-> global_cell_[origin_in_level_zero];
+    }
+}
+
 void CpGrid::getIJK(const int c, std::array<int,3>& ijk) const
 {
     current_view_data_->getIJK(c, ijk);
@@ -1809,7 +1855,6 @@ bool CpGrid::adapt(const std::vector<std::array<int,3>>& cells_per_dim_vec,
                                     cornerInMarkedElemWithEquivRefinedCorner,
                                     cells_per_dim_vec);
 
-    std::vector<int> adapted_global_cell(cell_count, 0);
     updateLeafGridViewGeometries( /* Leaf grid View Corners arguments */
                                   adapted_corners,
                                   corner_count,
@@ -1822,7 +1867,6 @@ bool CpGrid::adapt(const std::vector<std::array<int,3>>& cells_per_dim_vec,
                                   /* Leaf grid View Cells argumemts  */
                                   adapted_cells,
                                   adapted_cell_to_point,
-                                  adapted_global_cell,
                                   cell_count,
                                   adapted_cell_to_face,
                                   adapted_face_to_cell,
@@ -1895,7 +1939,6 @@ bool CpGrid::adapt(const std::vector<std::array<int,3>>& cells_per_dim_vec,
         (*data[refinedLevelGridIdx]).child_to_parent_cells_ = refined_child_to_parent_cells_vec[level];
         (*data[refinedLevelGridIdx]).cell_to_idxInParentCell_ = refined_cell_to_idxInParentCell_vec[level];
         (*data[refinedLevelGridIdx]).level_to_leaf_cells_ =  refined_level_to_leaf_cells_vec[level];
-        (*data[refinedLevelGridIdx]).global_cell_.swap(refined_global_cell_vec[level]);
         (*data[refinedLevelGridIdx]).index_set_ = std::make_unique<cpgrid::IndexSet>(data[refinedLevelGridIdx]->size(0),
                                                                                      data[refinedLevelGridIdx]->size(3));
         // Determine the amount of cells per direction, per parent cell, of the corresponding LGR.
@@ -1912,6 +1955,7 @@ bool CpGrid::adapt(const std::vector<std::array<int,3>>& cells_per_dim_vec,
         }
         else {
             (*data[refinedLevelGridIdx]).logical_cartesian_size_ = (*data[0]).logical_cartesian_size_;
+            (*data[refinedLevelGridIdx]).global_cell_.swap(refined_global_cell_vec[level]);
         }
         // One alternative definition for logical_cartesian_size_ in the case where the marked elements for refinement do not form a block of cells,
         // therefore, are not associated with the keyword CARFIN, is to imagine that we put all the marked elements one next to the other, along
@@ -1927,14 +1971,31 @@ bool CpGrid::adapt(const std::vector<std::array<int,3>>& cells_per_dim_vec,
     (*data[levels + preAdaptMaxLevel +1]).child_to_parent_cells_ = adapted_child_to_parent_cells;
     (*data[levels + preAdaptMaxLevel +1]).cell_to_idxInParentCell_ = adapted_cell_to_idxInParentCell;
     (*data[levels + preAdaptMaxLevel +1]).leaf_to_level_cells_ =  leaf_to_level_cells;
-    (*data[levels + preAdaptMaxLevel +1]).global_cell_.swap(adapted_global_cell);
     (*data[levels + preAdaptMaxLevel +1]).index_set_ = std::make_unique<cpgrid::IndexSet>(data[levels + preAdaptMaxLevel +1]->size(0),
                                                                                           data[levels + preAdaptMaxLevel +1]->size(3));
     (*data[levels + preAdaptMaxLevel +1]).logical_cartesian_size_ =  (*data[0]).logical_cartesian_size_;
 
     // Update the leaf grid view
     current_view_data_ = data.back().get();
 
+    // When the refinement is determined by startIJK and endIJK values, the LGR has a (local) Cartesian size.
+    // Therefore, each refined cell belonging to the LGR can be associated with a (local) IJK and its (local) Cartesian index.
+    // If the LGR has NXxNYxNZ dimension, then the Cartesian indices take values
+    // k*NN*NY + j*NX + i, where i<NX, j<Ny, k<NZ.
+    // This index is stored in <refined-level-grid>.global_cell_[ refined cell index (~element.index()) ] =  k*NN*NY + j*NX + i.
+    if (isCARFIN) {
+        for (int level = 0; level < levels; ++level) {
+            const int refinedLevelGridIdx = level + preAdaptMaxLevel +1;
+            std::vector<int> global_cell_lgr(data[refinedLevelGridIdx]->size(0));
+            computeGlobalCellLgr(refinedLevelGridIdx, startIJK_vec[level], global_cell_lgr);
+            (*data[refinedLevelGridIdx]).global_cell_.swap(global_cell_lgr);
+        }
+    }
+
+    std::vector<int> global_cell_leaf( data[levels + preAdaptMaxLevel +1]->size(0));
+    computeGlobalCellLeafGridViewWithLgrs(global_cell_leaf);
+    (*data[levels + preAdaptMaxLevel +1]).global_cell_.swap(global_cell_leaf);
+
     updateCornerHistoryLevels(cornerInMarkedElemWithEquivRefinedCorner,
                               elemLgrAndElemLgrCorner_to_refinedLevelAndRefinedCorner,
                               adaptedCorner_to_elemLgrAndElemLgrCorner,
@@ -2549,7 +2610,7 @@ CpGrid::defineLevelToLeafAndLeafToLevelCells(const std::map<std::array<int,2>,st
                                              const int& cell_count) const
 {
     // If the (level zero) grid has been distributed, then the preAdaptGrid is data_[0]. Otherwise, preApaptGrid is current_view_data_.
-    
+
     // -- Refined to Adapted cells and Adapted-cells to {level where the cell was born, cell index on that level} --
     // Relation between the refined grid and leafview cell indices.
     std::vector<std::vector<int>> refined_level_to_leaf_cells_vec(refined_cell_count_vec.size());
@@ -3237,7 +3298,6 @@ void CpGrid::populateRefinedFaces(std::vector<Dune::cpgrid::EntityVariableBase<c
 
 void CpGrid::populateLeafGridCells(Dune::cpgrid::EntityVariableBase<cpgrid::Geometry<3,3>>& adapted_cells,
                                    std::vector<std::array<int,8>>& adapted_cell_to_point,
-                                   std::vector<int>& adapted_global_cell,
                                    const int& cell_count,
                                    cpgrid::OrientedEntityTable<0,1>& adapted_cell_to_face,
                                    cpgrid::OrientedEntityTable<1,0>& adapted_face_to_cell,
@@ -3265,8 +3325,6 @@ void CpGrid::populateLeafGridCells(Dune::cpgrid::EntityVariableBase<cpgrid::Geom
         const auto& [elemLgr, elemLgrCell] = adaptedCell_to_elemLgrAndElemLgrCell.at(cell);
         const auto& elemLgrCellEntity =  Dune::cpgrid::EntityRep<0>(elemLgrCell, true);
 
-        adapted_global_cell[cell] = current_view_data_->global_cell_[(elemLgr == -1) ? elemLgrCell : elemLgr];
-
         // Auxiliary cell_to_face
         std::vector<cpgrid::EntityRep<1>> aux_cell_to_face;
 
@@ -3574,7 +3632,6 @@ void CpGrid::updateLeafGridViewGeometries( /* Leaf grid View Corners arguments *
                                            /* Leaf grid View Cells argumemts  */
                                            Dune::cpgrid::EntityVariableBase<cpgrid::Geometry<3,3>>& adapted_cells,
                                            std::vector<std::array<int,8>>& adapted_cell_to_point,
-                                           std::vector<int>& adapted_global_cell,
                                            const int& cell_count,
                                            cpgrid::OrientedEntityTable<0,1>& adapted_cell_to_face,
                                            cpgrid::OrientedEntityTable<1,0>& adapted_face_to_cell,
@@ -3617,7 +3674,6 @@ void CpGrid::updateLeafGridViewGeometries( /* Leaf grid View Corners arguments *
     // --- Adapted cells ---
     populateLeafGridCells(adapted_cells,
                           adapted_cell_to_point,
-                          adapted_global_cell,
                           cell_count,
                           adapted_cell_to_face,
                           adapted_face_to_cell,

opm/grid/cpgrid/CpGridData.cpp

@@ -1710,16 +1710,10 @@ void CpGridData::computeCommunicationInterfaces([[maybe_unused]] int noExistingP
 #endif
 }
 
-std::array<Dune::FieldVector<double,3>,8> CpGridData::getReferenceRefinedCorners(int idxInParentCell, const std::array<int,3>& cells_per_dim) const
+std::array<Dune::FieldVector<double,3>,8> CpGridData::getReferenceRefinedCorners(int idx_in_parent_cell, const std::array<int,3>& cells_per_dim) const
 {
     // Refined cells in parent cell: k*cells_per_dim[0]*cells_per_dim[1] + j*cells_per_dim[0] + i
-    std::array<int,3> ijk = {0,0,0};
-    ijk[0] = idxInParentCell % cells_per_dim[0];
-    idxInParentCell -= ijk[0]; // k*cells_per_dim[0]*cells_per_dim[1] + j*cells_per_dim[0]
-    idxInParentCell /= cells_per_dim[0]; // k*cells_per_dim[1] + j
-    ijk[1] = idxInParentCell % cells_per_dim[1];
-    idxInParentCell -= ijk[1]; // k*cells_per_dim[1]
-    ijk[2] = idxInParentCell /cells_per_dim[1];
+    std::array<int,3> ijk = getIJK(idx_in_parent_cell, cells_per_dim);
 
     std::array<Dune::FieldVector<double,3>,8> corners_in_parent_reference_elem = { // corner '0'
         {{ double(ijk[0])/cells_per_dim[0], double(ijk[1])/cells_per_dim[1], double(ijk[2])/cells_per_dim[2] },

opm/grid/cpgrid/CpGridData.hpp

@@ -303,10 +303,28 @@ class CpGridData
     /// @param [out] ijk  Cartesian index triplet
     void getIJK(int c, std::array<int,3>& ijk) const
     {
-        int gc = global_cell_[c];
-        ijk[0] = gc % logical_cartesian_size_[0];  gc /= logical_cartesian_size_[0];
-        ijk[1] = gc % logical_cartesian_size_[1];
-        ijk[2] = gc / logical_cartesian_size_[1];
+        ijk = getIJK(global_cell_[c], logical_cartesian_size_);
+    }
+
+    /// @brief Extract Cartesian index triplet (i,j,k) given an index between 0 and NXxNYxNZ -1
+    ///        where NX, NY, and NZ is the total amoung of cells in each direction x-,y-,and z- respectively.
+    ///
+    /// @param [in] idx      Integer between 0 and cells_per_dim[0]*cells_per_dim[1]*cells_per_dim[2]-1
+    /// @param [in] cells_per_dim
+    /// @return Cartesian index triplet.
+    std::array<int,3> getIJK(int idx_in_parent_cell, const std::array<int,3>& cells_per_dim) const
+    {
+        // idx = k*cells_per_dim_[0]*cells_per_dim_[1] + j*cells_per_dim_[0] + i
+        // with 0<= i < cells_per_dim_[0], 0<= j < cells_per_dim_[1], 0<= k <cells_per_dim_[2].
+        assert(cells_per_dim[0]);
+        assert(cells_per_dim[1]);
+        assert(cells_per_dim[2]);
+
+        std::array<int,3> ijk = {0,0,0};
+        ijk[0] = idx_in_parent_cell % cells_per_dim[0]; idx_in_parent_cell /= cells_per_dim[0];
+        ijk[1] = idx_in_parent_cell % cells_per_dim[1];
+        ijk[2] = idx_in_parent_cell /cells_per_dim[1];
+        return ijk;
     }
 
     /// @brief Determine if a finite amount of patches (of cells) are disjoint, namely, they do not share any corner nor face.
@@ -381,7 +399,7 @@ class CpGridData
     void postAdapt();
 
 private:
-    std::array<Dune::FieldVector<double,3>,8> getReferenceRefinedCorners(int idxInParentCell, const std::array<int,3>& cells_per_dim) const;
+    std::array<Dune::FieldVector<double,3>,8> getReferenceRefinedCorners(int idx_in_parent_cell, const std::array<int,3>& cells_per_dim) const;
 
     /// @brief Compute amount of cells in each direction of a patch of cells. (Cartesian grid required).
     ///

opm/grid/cpgrid/Entity.hpp

@@ -529,11 +529,10 @@ Dune::cpgrid::Geometry<3,3> Dune::cpgrid::Entity<codim>::geometryInFather() cons
     if (!(this->hasFather())){
         OPM_THROW(std::logic_error, "Entity has no father.");
     }
-    if (pgrid_ -> cell_to_idxInParentCell_[this->index()] !=-1) {
-        int idxInParentCell = pgrid_ -> cell_to_idxInParentCell_[this->index()];
-        assert(idxInParentCell>-1);
+    auto idx_in_parent_cell = pgrid_ -> cell_to_idxInParentCell_[this->index()];
+    if (idx_in_parent_cell !=-1) {
         const auto& cells_per_dim =  (*(pgrid_ -> level_data_ptr_))[this->level()] -> cells_per_dim_;
-        const auto& auxArr = pgrid_ -> getReferenceRefinedCorners(idxInParentCell, cells_per_dim);
+        const auto& auxArr = pgrid_ -> getReferenceRefinedCorners(idx_in_parent_cell, cells_per_dim);
         FieldVector<double, 3> corners_in_father_reference_elem_temp[8] =
             { auxArr[0], auxArr[1], auxArr[2], auxArr[3], auxArr[4], auxArr[5], auxArr[6], auxArr[7]};
         auto in_father_reference_elem_corners = std::make_shared<EntityVariable<cpgrid::Geometry<0, 3>, 3>>();
@@ -614,10 +613,9 @@ Dune::cpgrid::Entity<0> Dune::cpgrid::Entity<codim>::getEquivLevelElem() const
 template<int codim>
 int Dune::cpgrid::Entity<codim>::getLevelCartesianIdx() const
 {
-    const auto entityLevel = this -> level();
-    const auto level = (*(pgrid_ -> level_data_ptr_))[entityLevel].get();
-    const auto& elemInLevel = this->getLevelElem(); // throws when the entity does not belong to the leaf grid view.
-    return level -> global_cell_[elemInLevel.index()];
+    const auto& level_data = (*(pgrid_ -> level_data_ptr_))[level()].get();
+    // getLevelElem() throws when the entity does not belong to the leaf grid view.
+    return level_data -> global_cell_[getLevelElem().index()];
 }
 
 } // namespace cpgrid