diff --git a/opm/grid/CpGrid.hpp b/opm/grid/CpGrid.hpp index a188dca8d..c602d63fa 100644 --- a/opm/grid/CpGrid.hpp +++ b/opm/grid/CpGrid.hpp @@ -889,10 +889,9 @@ namespace Dune const std::vector>& 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>& adapted_cells, std::vector>& adapted_cell_to_point, - std::vector& 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>& adapted_cells, std::vector>& adapted_cell_to_point, - std::vector& 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& startIJK, std::vector& 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& 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 diff --git a/opm/grid/cpgrid/CartesianIndexMapper.hpp b/opm/grid/cpgrid/CartesianIndexMapper.hpp index 05f28affc..2c50cadb8 100644 --- a/opm/grid/cpgrid/CartesianIndexMapper.hpp +++ b/opm/grid/cpgrid/CartesianIndexMapper.hpp @@ -52,7 +52,6 @@ namespace Dune int compressedLevelZeroSize() const { - return (*grid_.currentData()[0]).size(0); } diff --git a/opm/grid/cpgrid/CpGrid.cpp b/opm/grid/cpgrid/CpGrid.cpp index 5188020d3..d48b85105 100644 --- a/opm/grid/cpgrid/CpGrid.cpp +++ b/opm/grid/cpgrid/CpGrid.cpp @@ -633,6 +633,52 @@ const std::vector& CpGrid::globalCell() const return currentData().back() -> global_cell_; } +void CpGrid::computeGlobalCellLgr(const int& level, const std::array& startIJK, std::vector& 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 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 childIJK = currentData()[level]-> getIJK(idx_in_parent_cell, cells_per_dim); + // The corresponding lgrIJK can be computed as follows: + const std::array& 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& 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& ijk) const { current_view_data_->getIJK(c, ijk); @@ -1809,7 +1855,6 @@ bool CpGrid::adapt(const std::vector>& cells_per_dim_vec, cornerInMarkedElemWithEquivRefinedCorner, cells_per_dim_vec); - std::vector 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>& 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>& 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(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>& 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,7 +1971,6 @@ bool CpGrid::adapt(const std::vector>& 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(data[levels + preAdaptMaxLevel +1]->size(0), data[levels + preAdaptMaxLevel +1]->size(3)); (*data[levels + preAdaptMaxLevel +1]).logical_cartesian_size_ = (*data[0]).logical_cartesian_size_; @@ -1935,6 +1978,24 @@ bool CpGrid::adapt(const std::vector>& cells_per_dim_vec, // 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.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 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 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,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> refined_level_to_leaf_cells_vec(refined_cell_count_vec.size()); @@ -3237,7 +3298,6 @@ void CpGrid::populateRefinedFaces(std::vector>& adapted_cells, std::vector>& adapted_cell_to_point, - std::vector& 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(elemLgrCell, true); - adapted_global_cell[cell] = current_view_data_->global_cell_[(elemLgr == -1) ? elemLgrCell : elemLgr]; - // Auxiliary cell_to_face std::vector> aux_cell_to_face; @@ -3574,7 +3632,6 @@ void CpGrid::updateLeafGridViewGeometries( /* Leaf grid View Corners arguments * /* Leaf grid View Cells argumemts */ Dune::cpgrid::EntityVariableBase>& adapted_cells, std::vector>& adapted_cell_to_point, - std::vector& 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, diff --git a/opm/grid/cpgrid/CpGridData.cpp b/opm/grid/cpgrid/CpGridData.cpp index f9afd1c17..d63d80d3f 100644 --- a/opm/grid/cpgrid/CpGridData.cpp +++ b/opm/grid/cpgrid/CpGridData.cpp @@ -1710,16 +1710,10 @@ void CpGridData::computeCommunicationInterfaces([[maybe_unused]] int noExistingP #endif } -std::array,8> CpGridData::getReferenceRefinedCorners(int idxInParentCell, const std::array& cells_per_dim) const +std::array,8> CpGridData::getReferenceRefinedCorners(int idx_in_parent_cell, const std::array& 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 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 ijk = getIJK(idx_in_parent_cell, cells_per_dim); std::array,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] }, diff --git a/opm/grid/cpgrid/CpGridData.hpp b/opm/grid/cpgrid/CpGridData.hpp index 5bd148dfc..136da9a38 100644 --- a/opm/grid/cpgrid/CpGridData.hpp +++ b/opm/grid/cpgrid/CpGridData.hpp @@ -303,10 +303,28 @@ class CpGridData /// @param [out] ijk Cartesian index triplet void getIJK(int c, std::array& 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 getIJK(int idx_in_parent_cell, const std::array& 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 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,8> getReferenceRefinedCorners(int idxInParentCell, const std::array& cells_per_dim) const; + std::array,8> getReferenceRefinedCorners(int idx_in_parent_cell, const std::array& cells_per_dim) const; /// @brief Compute amount of cells in each direction of a patch of cells. (Cartesian grid required). /// diff --git a/opm/grid/cpgrid/Entity.hpp b/opm/grid/cpgrid/Entity.hpp index 737764fec..ac95f51c0 100644 --- a/opm/grid/cpgrid/Entity.hpp +++ b/opm/grid/cpgrid/Entity.hpp @@ -529,11 +529,10 @@ Dune::cpgrid::Geometry<3,3> Dune::cpgrid::Entity::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 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, 3>>(); @@ -614,10 +613,9 @@ Dune::cpgrid::Entity<0> Dune::cpgrid::Entity::getEquivLevelElem() const template int Dune::cpgrid::Entity::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 diff --git a/tests/cpgrid/adapt_cpgrid_test.cpp b/tests/cpgrid/adapt_cpgrid_test.cpp index e62bfa3a5..6d08ae94b 100644 --- a/tests/cpgrid/adapt_cpgrid_test.cpp +++ b/tests/cpgrid/adapt_cpgrid_test.cpp @@ -105,6 +105,7 @@ void markAndAdapt_check(Dune::CpGrid& coarse_grid, BOOST_CHECK(static_cast(data.size()) == coarse_grid.maxLevel() +2); const auto& leafGridIdx = coarse_grid.maxLevel() +1; const auto& adapted_leaf = *data[leafGridIdx]; + if(isBlockShape) { // For a mixed grid that gets refined a second time, isBlockShape == false, even though the marked elements form a block. const auto& blockRefinement_data = other_grid.currentData(); const auto& blockRefinement_leaf = *blockRefinement_data.back(); @@ -207,6 +208,12 @@ void markAndAdapt_check(Dune::CpGrid& coarse_grid, const auto& grid_view = coarse_grid.leafGridView(); Dune::MultipleCodimMultipleGeomTypeMapper adaptMapper(grid_view, Dune::mcmgElementLayout()); + auto itMin = std::min_element((data.back() -> global_cell_).begin(), (data.back()-> global_cell_).end()); + auto itMax = std::max_element((data.back() -> global_cell_).begin(), (data.back() -> global_cell_).end()); + BOOST_CHECK_EQUAL( *itMin, 0); + const auto& maxCartesianIdx = coarse_grid.logicalCartesianSize()[0]*coarse_grid.logicalCartesianSize()[1]*coarse_grid.logicalCartesianSize()[2] -1; + BOOST_CHECK_EQUAL( *itMax, maxCartesianIdx); + for(const auto& element: elements(grid_view)) { // postAdapt() has been called, therefore every element gets marked with 0 BOOST_CHECK( coarse_grid.getMark(element) == 0); @@ -241,10 +248,6 @@ void markAndAdapt_check(Dune::CpGrid& coarse_grid, BOOST_CHECK_EQUAL( child_to_parent[1], element.father().index()); BOOST_CHECK( element.father() == element.getOrigin()); - const auto& auxLevel = (element.father().level() == 0) ? element.getOrigin().level() : element.getLevelElem().level(); - const auto& auxLevelIdx = (element.father().level() == 0) ? element.getOrigin().index() : element.getLevelElem().index(); - BOOST_CHECK( ( adapted_leaf.global_cell_[element.index()]) == (data[auxLevel]->global_cell_[auxLevelIdx]) ); - BOOST_CHECK( std::get<0>(data[element.father().level()]->parent_to_children_cells_[element.father().index()]) == element.level()); // Check amount of children cells of the parent cell BOOST_CHECK_EQUAL(std::get<1>(data[element.father().level()]->parent_to_children_cells_[child_to_parent[1]]).size(), @@ -370,6 +373,15 @@ void markAndAdapt_check(Dune::CpGrid& coarse_grid, } // end-preAdaptElements-for-loop } // end-startingGridIdx==0 + for (int level = 1; level < coarse_grid.maxLevel(); ++level) + { + auto itMinLevel = std::min_element((data[level] -> global_cell_).begin(), (data[level] -> global_cell_).end()); + auto itMaxLevel = std::max_element((data[level] -> global_cell_).begin(), (data[level] -> global_cell_).end()); + BOOST_CHECK_EQUAL( *itMinLevel, 0); + const auto& maxCartesianIdxLevel = data[level]->logical_cartesian_size_[0]*data[level]->logical_cartesian_size_[1]* data[level]->logical_cartesian_size_[2] -1; + BOOST_CHECK_EQUAL( *itMaxLevel, maxCartesianIdxLevel); + } + std::set allIds_set; std::vector allIds_vec; allIds_vec.reserve(data.back()->size(0) + data.back()->size(3)); diff --git a/tests/cpgrid/addLgrsOnDistributedGrid_test.cpp b/tests/cpgrid/addLgrsOnDistributedGrid_test.cpp index 5215886ce..8b0696bf3 100644 --- a/tests/cpgrid/addLgrsOnDistributedGrid_test.cpp +++ b/tests/cpgrid/addLgrsOnDistributedGrid_test.cpp @@ -175,14 +175,23 @@ void refinePatch_and_check(Dune::CpGrid& coarse_grid, } BOOST_CHECK( entity.level() == 0); } - + + if (!(data[level] -> global_cell_.empty())) + { + auto itMin = std::min_element((data[level] -> global_cell_).begin(), (data[level] -> global_cell_).end()); + auto itMax = std::max_element((data[level] -> global_cell_).begin(), (data[level] -> global_cell_).end()); + BOOST_CHECK_EQUAL( *itMin, 0); + const auto& maxCartesianIdxLevel = data[level]->logical_cartesian_size_[0]*data[level]->logical_cartesian_size_[1]* data[level]->logical_cartesian_size_[2]; + BOOST_CHECK( *itMax < maxCartesianIdxLevel); + } + + // LGRs for (int cell = 0; cell < data[level]-> size(0); ++cell) { Dune::cpgrid::Entity<0> entity = Dune::cpgrid::Entity<0>(*data[level], cell, true); BOOST_CHECK( entity.hasFather() == true); BOOST_CHECK( entity.getOrigin() == entity.father()); - BOOST_CHECK( entity.index() == (data[level] -> global_cell_[entity.index()])); // global_cell_ = {0,1,..., total cells -1} BOOST_CHECK( entity.getOrigin().level() == 0); BOOST_CHECK_CLOSE(entity.geometryInFather().volume(), 1./(cells_per_dim_vec[level-1][0]*cells_per_dim_vec[level-1][1]*cells_per_dim_vec[level-1][2]), 1e-6); @@ -220,6 +229,12 @@ void refinePatch_and_check(Dune::CpGrid& coarse_grid, BOOST_CHECK((*data[startIJK_vec.size() +1]).face_to_cell_[faceEntity].size() < 3); } + auto itMin = std::min_element((data.back() -> global_cell_).begin(), (data.back()-> global_cell_).end()); + auto itMax = std::max_element((data.back() -> global_cell_).begin(), (data.back() -> global_cell_).end()); + BOOST_CHECK( *itMin >= 0); + const auto& maxCartesianIdx = coarse_grid.logicalCartesianSize()[0]*coarse_grid.logicalCartesianSize()[1]*coarse_grid.logicalCartesianSize()[2]; + BOOST_CHECK( *itMax < maxCartesianIdx); + // LeafView for (int cell = 0; cell < data[startIJK_vec.size()+1]-> size(0); ++cell) { @@ -249,8 +264,6 @@ void refinePatch_and_check(Dune::CpGrid& coarse_grid, BOOST_CHECK_EQUAL( child_to_parent[0] == 0, true); BOOST_CHECK_EQUAL( child_to_parent[1], entity.father().index()); BOOST_CHECK( entity.father() == entity.getOrigin()); - BOOST_CHECK( (data[startIJK_vec.size() +1] -> global_cell_[entity.index()]) == - (data[0] -> global_cell_[entity.getOrigin().index()]) ); BOOST_CHECK( entity.getOrigin().level() == 0); BOOST_CHECK( std::get<0>((*data[0]).parent_to_children_cells_[child_to_parent[1]]) == entity.level()); BOOST_CHECK_EQUAL((std::find(std::get<1>((*data[0]).parent_to_children_cells_[child_to_parent[1]]).begin(), @@ -285,6 +298,13 @@ void refinePatch_and_check(Dune::CpGrid& coarse_grid, BOOST_CHECK( static_cast(startIJK_vec.size()) == coarse_grid.maxLevel()); BOOST_CHECK( (*data[data.size()-1]).parent_to_children_cells_.empty()); + auto it_min = std::min_element((data.back() -> global_cell_).begin(), (data.back()-> global_cell_).end()); + auto it_max = std::max_element((data.back() -> global_cell_).begin(), (data.back() -> global_cell_).end()); + auto it_min_level_zero = std::min_element((data.front() -> global_cell_).begin(), (data.front() -> global_cell_).end()); + auto it_max_level_zero = std::max_element((data.front() -> global_cell_).begin(), (data.front() -> global_cell_).end()); + BOOST_CHECK_EQUAL( *it_min, *it_min_level_zero); + BOOST_CHECK_EQUAL( *it_max, *it_max_level_zero); + for (long unsigned int l = 0; l < startIJK_vec.size() +1; ++l) // level 0,1,2,... , last patch { const auto& view = coarse_grid.levelGridView(l); diff --git a/tests/cpgrid/grid_lgr_test.cpp b/tests/cpgrid/grid_lgr_test.cpp index f6ffe89f0..93357f8c8 100644 --- a/tests/cpgrid/grid_lgr_test.cpp +++ b/tests/cpgrid/grid_lgr_test.cpp @@ -124,6 +124,84 @@ void refinePatch_and_check(Dune::CpGrid& coarse_grid, BOOST_CHECK(coarse_grid.getLgrNameToLevel().at("GLOBAL") == 0); const auto& all_parent_cell_indices = (*data[0]).getPatchesCells(startIJK_vec, endIJK_vec); + // GLOBAL grid + for (int cell = 0; cell < data[0]-> size(0); ++cell) + { + Dune::cpgrid::Entity<0> entity = Dune::cpgrid::Entity<0>(*data[0], cell, true); + BOOST_CHECK( entity.hasFather() == false); + BOOST_CHECK_THROW(entity.father(), std::logic_error); + BOOST_CHECK_THROW(entity.geometryInFather(), std::logic_error); + BOOST_CHECK( entity.getOrigin() == entity); + BOOST_CHECK( entity.getOrigin().level() == 0); + + auto it = entity.hbegin(coarse_grid.maxLevel()); // With element.level(), fails + auto endIt = entity.hend(coarse_grid.maxLevel()); + const auto& [lgr, childrenList] = (*data[0]).parent_to_children_cells_[entity.index()]; + // If it == endIt, then entity.isLeaf() true (when dristibuted_data_ is empty) + if (it == endIt){ + BOOST_CHECK_EQUAL(lgr, -1); + BOOST_CHECK(childrenList.empty()); + BOOST_CHECK( entity.isLeaf() == true); + BOOST_CHECK( entity.mightVanish() == false); + } + else{ + BOOST_CHECK(lgr != -1); + // If it != endIt, then entity.isLeaf() false (when dristibuted_data_ is empty) + BOOST_CHECK_EQUAL( it == endIt, false); + BOOST_CHECK( entity.mightVanish() == true); + BOOST_CHECK( entity.isNew() == false); + BOOST_CHECK_EQUAL( entity.isLeaf(), false); // parent cells do not appear in the LeafView + + double referenceElemOneParent_volume = 0.; + std::array referenceElem_entity_center = {0.,0.,0.}; // Expected {.5,.5,.5} + for (const auto& child : childrenList) { + BOOST_CHECK( child != -1); + BOOST_CHECK( data[lgr]-> child_to_parent_cells_[child][0] == 0); + BOOST_CHECK( data[lgr]-> child_to_parent_cells_[child][1] == entity.index()); + + const auto& childElem = Dune::cpgrid::Entity<0>(*data[lgr], child, true); + const auto& childGeoInFather = childElem.geometryInFather(); + + BOOST_CHECK(childElem.hasFather() == true); + BOOST_CHECK(childElem.level() == lgr); + referenceElemOneParent_volume += childGeoInFather.volume(); + for (int c = 0; c < 3; ++c) { + referenceElem_entity_center[c] += childGeoInFather.center()[c]; + } + } + + double referenceElemOneParent_volume_it = 0.; + std::array referenceElem_entity_center_it = {0.,0.,0.}; // Expected {.5,.5,.5} + for (; it != endIt; ++it) + { + BOOST_CHECK(it ->hasFather() == true); + BOOST_CHECK(it ->level() == lgr); + const auto& geoInFather = it-> geometryInFather(); + + referenceElemOneParent_volume_it += geoInFather.volume(); + for (int c = 0; c < 3; ++c) + { + referenceElem_entity_center_it[c] += geoInFather.center()[c]; + } + } + for (int c = 0; c < 3; ++c) { + referenceElem_entity_center[c] /= cells_per_dim_vec[lgr-1][0]*cells_per_dim_vec[lgr-1][1]*cells_per_dim_vec[lgr-1][2]; + referenceElem_entity_center_it[c] /= cells_per_dim_vec[lgr-1][0]*cells_per_dim_vec[lgr-1][1]*cells_per_dim_vec[lgr-1][2]; + } + BOOST_CHECK_CLOSE(referenceElemOneParent_volume, 1, 1e-13); + BOOST_CHECK_CLOSE(referenceElem_entity_center[0], .5, 1e-13); + BOOST_CHECK_CLOSE(referenceElem_entity_center[1], .5, 1e-13); + BOOST_CHECK_CLOSE(referenceElem_entity_center[2], .5, 1e-13); + + BOOST_CHECK_CLOSE(referenceElemOneParent_volume_it, 1, 1e-13); + BOOST_CHECK_CLOSE(referenceElem_entity_center_it[0], .5, 1e-13); + BOOST_CHECK_CLOSE(referenceElem_entity_center_it[1], .5, 1e-13); + BOOST_CHECK_CLOSE(referenceElem_entity_center_it[2], .5, 1e-13); + } + BOOST_CHECK( entity.level() == 0); + } // end-GLOBAL/LevelZero-grid + + for (long unsigned int level = 1; level < startIJK_vec.size() +1; ++level) // only 1 when there is only 1 patch { check_refinedPatch_grid(cells_per_dim_vec[level-1], startIJK_vec[level-1], endIJK_vec[level-1], @@ -133,81 +211,11 @@ void refinePatch_and_check(Dune::CpGrid& coarse_grid, BOOST_CHECK( (*data[level]).parent_to_children_cells_.empty()); BOOST_CHECK(coarse_grid.getLgrNameToLevel().at(lgr_name_vec[level-1]) == static_cast(level)); - const auto& patch_cells = (*data[0]).getPatchCells(startIJK_vec[level-1], endIJK_vec[level-1]); - - // GLOBAL grid - for (int cell = 0; cell < data[0]-> size(0); ++cell) - { - Dune::cpgrid::Entity<0> entity = Dune::cpgrid::Entity<0>(*data[0], cell, true); - BOOST_CHECK( entity.hasFather() == false); - BOOST_CHECK_THROW(entity.father(), std::logic_error); - BOOST_CHECK_THROW(entity.geometryInFather(), std::logic_error); - BOOST_CHECK( entity.getOrigin() == entity); - BOOST_CHECK( entity.getOrigin().level() == 0); - auto it = entity.hbegin(coarse_grid.maxLevel()); - auto endIt = entity.hend(coarse_grid.maxLevel()); - const auto& [lgr, childrenList] = (*data[0]).parent_to_children_cells_[cell]; - if (std::find(all_parent_cell_indices.begin(), all_parent_cell_indices.end(), cell) == all_parent_cell_indices.end()){ - BOOST_CHECK_EQUAL(lgr, -1); - BOOST_CHECK(childrenList.empty()); - BOOST_CHECK( entity.isLeaf() == true); - // If it == endIt, then entity.isLeaf() true (when dristibuted_data_ is empty) - BOOST_CHECK( it == endIt); - } - else{ - BOOST_CHECK(lgr != -1); - BOOST_CHECK(childrenList.size() > 1); - // Auxiliary int to check amount of children - double referenceElemOneParent_volume = 0.; - std::array referenceElem_entity_center = {0.,0.,0.}; // Expected {.5,.5,.5} - for (const auto& child : childrenList) { - BOOST_CHECK( child != -1); - BOOST_CHECK( data[lgr]-> child_to_parent_cells_[child][0] == 0); - BOOST_CHECK( data[lgr]-> child_to_parent_cells_[child][1] == cell); - - const auto& childElem = Dune::cpgrid::Entity<0>(*data[lgr], child, true); - BOOST_CHECK(childElem.hasFather() == true); - BOOST_CHECK(childElem.level() == lgr); - referenceElemOneParent_volume += childElem.geometryInFather().volume(); - for (int c = 0; c < 3; ++c) { - referenceElem_entity_center[c] += (childElem.geometryInFather().center())[c]; - } - } - BOOST_CHECK_EQUAL( entity.isLeaf(), false); // parent cells do not appear in the LeafView - // Auxiliary int to check hierarchic iterator functionality - double referenceElemOneParent_volume_it = 0.; - std::array referenceElem_entity_center_it = {0.,0.,0.}; // Expected {.5,.5,.5} - // If it != endIt, then entity.isLeaf() false (when dristibuted_data_ is empty) - BOOST_CHECK( it != endIt ); - for (; it != endIt; ++it) - { - // Do something with the son available through it-> - BOOST_CHECK(it ->hasFather() == true); - BOOST_CHECK(it ->level() == lgr); - referenceElemOneParent_volume_it += it-> geometryInFather().volume(); - for (int c = 0; c < 3; ++c) - { - referenceElem_entity_center_it[c] += (it-> geometryInFather().center())[c]; - } - } - for (int c = 0; c < 3; ++c) - { - referenceElem_entity_center[c] - /= cells_per_dim_vec[lgr-1][0]*cells_per_dim_vec[lgr-1][1]*cells_per_dim_vec[lgr-1][2]; - referenceElem_entity_center_it[c] - /= cells_per_dim_vec[lgr-1][0]*cells_per_dim_vec[lgr-1][1]*cells_per_dim_vec[lgr-1][2]; - } - BOOST_CHECK_CLOSE(referenceElemOneParent_volume, 1, 1e-13); - BOOST_CHECK_CLOSE(referenceElem_entity_center[0], .5, 1e-13); - BOOST_CHECK_CLOSE(referenceElem_entity_center[1], .5, 1e-13); - BOOST_CHECK_CLOSE(referenceElem_entity_center[2], .5, 1e-13); - BOOST_CHECK_CLOSE(referenceElemOneParent_volume_it, 1, 1e-13); - BOOST_CHECK_CLOSE(referenceElem_entity_center_it[0], .5, 1e-13); - BOOST_CHECK_CLOSE(referenceElem_entity_center_it[1], .5, 1e-13); - BOOST_CHECK_CLOSE(referenceElem_entity_center_it[2], .5, 1e-13); - } - BOOST_CHECK( entity.level() == 0); - } + auto itMin = std::min_element((data[level] -> global_cell_).begin(), (data[level] -> global_cell_).end()); + auto itMax = std::max_element((data[level] -> global_cell_).begin(), (data[level] -> global_cell_).end()); + BOOST_CHECK_EQUAL( *itMin, 0); + const auto& maxCartesianIdxLevel = data[level]->logical_cartesian_size_[0]*data[level]->logical_cartesian_size_[1]* data[level]->logical_cartesian_size_[2] -1; + BOOST_CHECK_EQUAL( *itMax, maxCartesianIdxLevel); // LGRs for (int cell = 0; cell < data[level]-> size(0); ++cell) @@ -215,14 +223,10 @@ void refinePatch_and_check(Dune::CpGrid& coarse_grid, Dune::cpgrid::Entity<0> entity = Dune::cpgrid::Entity<0>(*data[level], cell, true); BOOST_CHECK( entity.hasFather() == true); BOOST_CHECK( entity.getOrigin() == entity.father()); - BOOST_CHECK( entity.index() == (data[level] -> global_cell_[entity.index()])); // global_cell_ = {0,1,..., total cells -1} BOOST_CHECK( entity.getOrigin().level() == 0); BOOST_CHECK_CLOSE(entity.geometryInFather().volume(), 1./(cells_per_dim_vec[level-1][0]*cells_per_dim_vec[level-1][1]*cells_per_dim_vec[level-1][2]), 1e-6); BOOST_CHECK(entity.father().level() == 0); - // Check entity.father().index() belongs to the patch_cells (corresponding LGR parents) - BOOST_CHECK_EQUAL((std::find(patch_cells.begin(), patch_cells.end(), - entity.father().index()) == patch_cells.end()) , false); const auto& child_to_parent = (*data[level]).child_to_parent_cells_[cell]; BOOST_CHECK_EQUAL( child_to_parent[0] == -1, false); BOOST_CHECK_EQUAL( child_to_parent[0] == 0, true); @@ -287,8 +291,6 @@ void refinePatch_and_check(Dune::CpGrid& coarse_grid, BOOST_CHECK_EQUAL( child_to_parent[0] == 0, true); BOOST_CHECK_EQUAL( child_to_parent[1], entity.father().index()); BOOST_CHECK( entity.father() == entity.getOrigin()); - BOOST_CHECK( (data[startIJK_vec.size() +1] -> global_cell_[entity.index()]) == - (data[0] -> global_cell_[entity.getOrigin().index()]) ); BOOST_CHECK( entity.getOrigin().level() == 0); BOOST_CHECK( std::get<0>((*data[0]).parent_to_children_cells_[child_to_parent[1]]) == entity.level()); BOOST_CHECK_EQUAL((std::find(std::get<1>((*data[0]).parent_to_children_cells_[child_to_parent[1]]).begin(), @@ -457,6 +459,12 @@ void refinePatch_and_check(Dune::CpGrid& coarse_grid, BOOST_CHECK( static_cast(startIJK_vec.size()) == coarse_grid.maxLevel()); BOOST_CHECK( (*data[data.size()-1]).parent_to_children_cells_.empty()); + auto itMin = std::min_element((data.back() -> global_cell_).begin(), (data.back()-> global_cell_).end()); + auto itMax = std::max_element((data.back() -> global_cell_).begin(), (data.back() -> global_cell_).end()); + BOOST_CHECK_EQUAL( *itMin, 0); + const auto& maxCartesianIdx = coarse_grid.logicalCartesianSize()[0]*coarse_grid.logicalCartesianSize()[1]*coarse_grid.logicalCartesianSize()[2] -1; + BOOST_CHECK_EQUAL( *itMax, maxCartesianIdx); + for (long unsigned int l = 0; l < startIJK_vec.size() +1; ++l) // level 0,1,2,... , last patch { const auto& view = coarse_grid.levelGridView(l); @@ -890,4 +898,3 @@ BOOST_AUTO_TEST_CASE(global_norefine) check_global_refine(coarse_grid, fine_grid); } - diff --git a/tests/cpgrid/inactiveCell_lgr_test.cpp b/tests/cpgrid/inactiveCell_lgr_test.cpp index d648a9c36..c73a1b311 100644 --- a/tests/cpgrid/inactiveCell_lgr_test.cpp +++ b/tests/cpgrid/inactiveCell_lgr_test.cpp @@ -173,13 +173,21 @@ void testInactiveCellsLgrs(const std::string& deckString, BOOST_CHECK( entity.level() == 0); } + if (!(data[level] -> global_cell_.empty())) + { + auto itMin = std::min_element((data[level] -> global_cell_).begin(), (data[level] -> global_cell_).end()); + auto itMax = std::max_element((data[level] -> global_cell_).begin(), (data[level] -> global_cell_).end()); + BOOST_CHECK( *itMin >= 0); + const auto& maxCartesianIdxLevel = data[level]->logical_cartesian_size_[0]*data[level]->logical_cartesian_size_[1]* data[level]->logical_cartesian_size_[2]; + BOOST_CHECK( *itMax < maxCartesianIdxLevel); + } + // LGRs for (int cell = 0; cell < data[level]-> size(0); ++cell) { Dune::cpgrid::Entity<0> entity = Dune::cpgrid::Entity<0>(*data[level], cell, true); BOOST_CHECK( entity.hasFather() == true); BOOST_CHECK( entity.getOrigin() == entity.father()); - BOOST_CHECK( entity.index() == (data[level] -> global_cell_[entity.index()])); // global_cell_ = {0,1,..., total cells -1} BOOST_CHECK( entity.getOrigin().level() == 0); BOOST_CHECK_CLOSE(entity.geometryInFather().volume(), 1./(cells_per_dim_vec[level-1][0]*cells_per_dim_vec[level-1][1]*cells_per_dim_vec[level-1][2]), 1e-6); @@ -285,6 +293,12 @@ void testInactiveCellsLgrs(const std::string& deckString, } } // end-level-for-loop + auto itMin = std::min_element((data.back() -> global_cell_).begin(), (data.back()-> global_cell_).end()); + auto itMax = std::max_element((data.back() -> global_cell_).begin(), (data.back() -> global_cell_).end()); + BOOST_CHECK( *itMin >= 0); + const auto& maxCartesianIdx = grid.logicalCartesianSize()[0]*grid.logicalCartesianSize()[1]*grid.logicalCartesianSize()[2]; + BOOST_CHECK( *itMax < maxCartesianIdx); + BOOST_CHECK( static_cast(startIJK_vec.size()) == grid.maxLevel()); BOOST_CHECK( (*data[data.size()-1]).parent_to_children_cells_.empty()); diff --git a/tests/cpgrid/lookupdataCpGrid_test.cpp b/tests/cpgrid/lookupdataCpGrid_test.cpp index 2722574d3..7acd83275 100644 --- a/tests/cpgrid/lookupdataCpGrid_test.cpp +++ b/tests/cpgrid/lookupdataCpGrid_test.cpp @@ -195,6 +195,7 @@ BOOST_AUTO_TEST_CASE(one_lgr_grid) const std::array cells_per_dim = {2,2,2}; const std::array startIJK = {1,0,1}; const std::array endIJK = {3,2,3}; // patch_dim = {3-1, 2-0, 3-1} ={2,2,2} + // Cells to be refined {13,14,17,18, 25,26,29,30} const std::string lgr_name = {"LGR1"}; grid.addLgrsUpdateLeafView({cells_per_dim}, {startIJK}, {endIJK}, {lgr_name}); @@ -353,7 +354,7 @@ void fieldProp_check(const Dune::CpGrid& grid, Opm::EclipseGrid eclGrid, std::st BOOST_AUTO_TEST_CASE(fieldProp) { - const std::string deckString = + const std::string deckString = R"( RUNSPEC DIMENS 1 1 5 / @@ -405,7 +406,7 @@ BOOST_AUTO_TEST_CASE(fieldProp) { BOOST_AUTO_TEST_CASE(fieldPropLgr) { - const std::string deckString = R"( RUNSPEC + const std::string deckString = R"( RUNSPEC DIMENS 1 1 5 /