WIP building levels aware of neighbor-cell-ref

Author: aritorto

opm/grid/cpgrid/CpGrid.cpp

@@ -2249,7 +2249,8 @@ bool CpGrid::refineAndUpdateGrid(bool throwOnFailure,
                                     cornerInMarkedElemWithEquivRefinedCorner,
                                     cells_per_dim_vec,
                                     preAdaptMaxLevel,
-                                    cellRefinements);
+                                    cellRefinements,
+                                    withoutFaults);
 
     for (int level = 0; level < levels; ++level) {
         const int refinedLevelGridIdx = level + preAdaptMaxLevel +1;

opm/grid/cpgrid/LgrFaultHelpers.cpp

@@ -718,10 +718,12 @@ void addVertices(const Dune::cpgrid::CpGridData&
     const auto& oldGridVertices = *(oldGrid.getGeometry().geomVector(std::integral_constant<int,3>()));
     Dune::cpgrid::EntityVariableBase<Dune::cpgrid::Geometry<0,3>>& newGridVertices = *(newGridGeomData.geometries.geomVector(std::integral_constant<int,3>()));
 
+    int newVerticesCount = 0;
     for (const auto& [parentFace, boundInfo] : parentCellFace_to_boundaryFaceInfo) {
-        newGridVertices.resize(numVertices + boundInfo.foundNewVertices.size());
+        newVerticesCount += boundInfo.foundNewVertices.size();
     }
-    
+    newGridVertices.resize(numVertices + newVerticesCount);
+      
     // add the "old" vertices (the parent-cell-face-aware are includeed here)
     for (int i = 0; i < oldGrid.size(3); ++i) {
         newGridVertices[i] = oldGridVertices.get(i);
@@ -732,7 +734,6 @@ void addVertices(const Dune::cpgrid::CpGridData&
             newVertex_to_newGridVertexIdx[vertex] = numVertices; 
             newGridVertices[numVertices] = Dune::cpgrid::Geometry<0, 3>(vertex);
             newGridBoundaryVertex_to_parentFace[numVertices] = parentFace;
-        
             ++numVertices;
         }
     }
@@ -965,8 +966,8 @@ void addFaces(int
 
     int newGrid_num_points = 0;
     int newGrid_face_count = 0;
-    
-    for (int i = 0; i < oldGrid.numFaces(); ++i) {
+
+  for (int i = 0; i < oldGrid.numFaces(); ++i) {
 
         const auto faceToCellSize = oldGrid.faceToCell(i).size();
 
@@ -1071,15 +1072,124 @@ void addFaces(int
                 }
             }
         }
+  }
+  
+
+    
+    
+  /*  for (int i = 0; i < oldGrid.numFaces(); ++i) {
+        
+        const auto faceToCellSize = oldGrid.faceToCell(i).size();
+        std::vector<int> newFaceIndices{};
+
+        if (faceToCellSize == 2) { // inner face->strore it
+            const auto face =  Dune::cpgrid::EntityRep<1>(i, true);
+               
+            old_to_new_faceIdxList[i] = std::vector<int>{newGrid_face_count};
+            new_to_old_faceIdx[newGrid_face_count] = i;
+
+            newGridFaces[newGrid_face_count] = (*oldGrid.getGeometry().geomVector(std::integral_constant<int,1>()))[face];
+            mutable_newGridFaceTags[newGrid_face_count] = oldGrid.faceTag(i);
+            mutable_newGridFaceNormals[newGrid_face_count] = oldGrid.faceNormals(i);
+
+            std::vector<int> faceToPoint{};
+            faceToPoint.reserve(oldGrid.faceToPoint(i).size());
+            for (const auto& vertexIdx : oldGrid.faceToPoint(i)) {
+                faceToPoint.push_back(vertexIdx);
+            }
+            newGrid_num_points += faceToPoint.size();
+            aux_newGrid_face_to_point[newGrid_face_count] = faceToPoint;
+
+            const auto& faceToCell = oldGrid.faceToCell(i);
+            newGridGeomData.face_to_cell.appendRow(faceToCell.begin(), faceToCell.end());
+
+            ++newGrid_face_count;
+        }
     }
+    for (int i = 0; i < oldGrid.numFaces(); ++i) {
+        const auto faceToCellSize = oldGrid.faceToCell(i).size();
+        if (faceToCellSize==2)
+            continue;
+
+        std::vector<int> newFaceIndices{};
+        
+        for (const auto& [parentFaceIdx, boundaryFaceInfo] : boundaryFacesInfo) {
+            const auto& newFacesInfo = boundaryFaceInfo.vanishedCellRefFace_to_newRefinedFaces[i];
+            for (const auto& [grid2_faceIdx, newFaceToCoord] : newFacesInfo) {
+
+                const auto overlapIt = boundaryFaceInfo.overlapFaces.find(grid2_faceIdx);
+                assert(overlapIt != boundaryFaceInfo.overlapFaces.end());
+                const auto& overlapFacesInfo = overlapIt->second;
+
+                newFaceIndices.push_back(newGrid_face_count);
+                old_to_new_faceIdxList[i].push_back(newGrid_face_count);
+                new_to_old_faceIdx[newGrid_face_count] = i;
+                newGridBoundaryFace_to_parentFace[newGrid_face_count] = parentFaceIdx;
+
+                const auto [faceCenter, faceArea, faceNormal] = computeFaceCenterAreaNormal(newFaceToCoord);
+
+                newGridFaces[newGrid_face_count] = Dune::cpgrid::Geometry<2,3>(faceCenter, faceArea);
+                mutable_newGridFaceTags[newGrid_face_count] = oldGrid.faceTag(i); // shared tag
+                mutable_newGridFaceNormals[newGrid_face_count] = faceNormal;
+
+                std::vector<int> faceToPoint{};
+                faceToPoint.reserve(newFaceToCoord.size());
+
+                for (const auto& vertex : newFaceToCoord) {
+                    int vertexIdx = -1; // invalid to be rewritten
+                    // check if vertex (is new) and has already been stored 
+                    auto it = newVertex_to_newGridVertexIdx.find(vertex);
+                    if (it != newVertex_to_newGridVertexIdx.end()) {
+                        vertexIdx = it->second;
+                    }
+                    else { // otherwise, vertex existed already in gridData1 (before "correction")
+                        auto iit = boundaryFaceInfo.boundaryVertex_to_vertexIdx.find(vertex);
+                        assert( iit != boundaryFaceInfo.boundaryVertex_to_vertexIdx.end());
+                        vertexIdx = iit->second;
+                    }
+                    assert(vertexIdx>=0);
+                    faceToPoint.push_back(vertexIdx);
+                }
+                newGrid_num_points += faceToPoint.size();
+                aux_newGrid_face_to_point[newGrid_face_count] = faceToPoint;
+
+                const auto& face_to_cell = oldGrid.faceToCell(i);
+                newGridGeomData.face_to_cell.appendRow(face_to_cell.begin(), face_to_cell.end());
+                ++newGrid_face_count;
+            }
+            if (newFacesInfo.empty()) { std::cout<< "empty for face i: " << i << std::endl;
+                const auto face =  Dune::cpgrid::EntityRep<1>(i, true);
+               
+                old_to_new_faceIdxList[i] = std::vector<int>{newGrid_face_count};
+                new_to_old_faceIdx[newGrid_face_count] = i;
+
+                newGridFaces[newGrid_face_count] = (*oldGrid.getGeometry().geomVector(std::integral_constant<int,1>()))[face];
+                mutable_newGridFaceTags[newGrid_face_count] = oldGrid.faceTag(i);
+                mutable_newGridFaceNormals[newGrid_face_count] = oldGrid.faceNormals(i);
+
+                std::vector<int> faceToPoint{};
+                faceToPoint.reserve(oldGrid.faceToPoint(i).size());
+                for (const auto& vertexIdx : oldGrid.faceToPoint(i)) {
+                    faceToPoint.push_back(vertexIdx);
+                }
+                newGrid_num_points += faceToPoint.size();
+                aux_newGrid_face_to_point[newGrid_face_count] = faceToPoint;
+
+                const auto& faceToCell = oldGrid.faceToCell(i);
+                newGridGeomData.face_to_cell.appendRow(faceToCell.begin(), faceToCell.end());
+
+                ++newGrid_face_count;
+            }
+            
+        }
+        }*/
     newGridGeomData.face_to_point.reserve(newGrid_face_count, newGrid_num_points);
     for (int face = 0; face < newGrid_face_count; ++face) {
         newGridGeomData.face_to_point.appendRow(aux_newGrid_face_to_point[face].begin(),
                                                 aux_newGrid_face_to_point[face].end());
     }
 }
 
-
 void addCells(GeomData& geomData,
               const Dune::cpgrid::CpGridData& gridData)
 {
@@ -1155,7 +1265,7 @@ void makeCellRefinementParentFaceAware(bool
     int numVertices = cellRefGrid.size(3);
     const auto& parentCellToPoint = parentGrid.cellToPoint(parentCellElem.index());
     std::map<int,int> correctedCellRefDataVertexIdx_to_parentGridVertexIdx{};
-
+     
     addVertices(numVertices,
                 cellRefGrid,
                 parentGrid_vertex_to_vertexIdx,
@@ -1327,21 +1437,21 @@ makeCellRefinementNeighborsAware(const Dune::cpgrid::Entity<0>& parentCell,
     int numFaces = cellRefGrid.numFaces();
 
     for (const auto& parentFace : parentCellToFace) {
-        
+      
         const auto faceToCell = parentGrid.faceToCell(parentFace.index());
 
-        if (faceToCell.size() == 1 || (faceInMarkedElemAndRefinedFaces[parentFace.index()].size()<=1)){
+        if (faceToCell.size() == 1) { // || (faceInMarkedElemAndRefinedFaces[parentFace.index()].size()<=1)){
             // face at parent grid boundary does need to be corrected (cell refinement is already parentCellFaces aware).
             continue;
         }
-
+        
         assert(faceToCell.size() == 2);
         assert(faceInMarkedElemAndRefinedFaces[parentFace.index()].size() == 2);
 
         const auto& [p1, refinedFaces1] = faceInMarkedElemAndRefinedFaces[parentFace.index()][0];
         const auto& [p2, refinedFaces2] = faceInMarkedElemAndRefinedFaces[parentFace.index()][1];
 
-        int cellRefIdx = (p1 == parentCell.index())? p1 : p2;
+        // int cellRefIdx = (p1 == parentCell.index())? p1 : p2;
         int neighborCellRefIdx = (p1 == parentCell.index())? p2 : p1;
 
         const auto& cellRefFaces = (p1 == parentCell.index())? refinedFaces1 : refinedFaces2;
@@ -1392,8 +1502,6 @@ makeCellRefinementNeighborsAware(const Dune::cpgrid::Entity<0>& parentCell,
                 boundaryInfo.newVertex_to_newGridVertexIdx,
                 cellRefGridBoundaryInfo.boundaryRefinedVertex_to_parentFace,
                 cellRefGridBoundaryInfo.boundaryRefinedVertexCoincidesWithParentVertex);
-   
-    std::cout<< numFaces << " numfaces before neighbor aware" << std::endl;
 
     addFaces(numFaces,
              cellRefGrid,
@@ -1407,7 +1515,7 @@ makeCellRefinementNeighborsAware(const Dune::cpgrid::Entity<0>& parentCell,
 
     addCells(neighborAwareGeomData, cellRefGrid);
 
-    std::cout<< neighborAwareCellRef.numFaces() << " numFaces after addition " << std::endl;
+ 
 
     return {neighborAwareCellRef_ptr, boundaryInfo};
 }
@@ -1418,6 +1526,8 @@ void makeCellRefinementsNeighborsAware(std::vector<std::shared_ptr<Dune::cpgrid:
                                        std::vector<std::vector<std::pair<int, std::vector<int>>>>& faceInMarkedElemAndRefinedFaces,
                                        std::vector<BoundaryInfo>& boundaryInfoCellRefs) // needs to be updated
 {
+    const auto cellRefsBeforeCorrection = cellRefs; // copy on propose
+    
     for (int elemIdx = 0; elemIdx < parentGrid.size(0); ++elemIdx) {
         auto& cellRef = cellRefs[elemIdx];
         if (cellRef == nullptr)
@@ -1426,7 +1536,7 @@ void makeCellRefinementsNeighborsAware(std::vector<std::shared_ptr<Dune::cpgrid:
         const auto& parentCell = Dune::cpgrid::Entity<0>(parentGrid, elemIdx, true);
 
         auto [correctedCellRef, correctedBoundaryInfo ] =  makeCellRefinementNeighborsAware(parentCell,
-                                                                          cellRefs, // to get neighbor information
+                                                                         cellRefsBeforeCorrection, // to get neighbor information
                                                                           parentGrid,
                                                                           *cellRef,
                                                                           cellRefsBoundaryInfo[elemIdx],
@@ -1468,7 +1578,7 @@ bool areClose(const Dune::FieldVector<double,3>& v,
 
 
 std::optional<int> getNeighborCellRefVertexIdx(const Dune::cpgrid::CpGridData& neighborCellRef,
-                                const std::vector<int>& faceList,
+                                               const std::vector<int>& faceList,
                                                const Dune::FieldVector<double,3>& vertex)
 {
     for (const auto& faceIdx : faceList) {
@@ -1493,7 +1603,7 @@ std::optional<int> getNeighborCellRefFaceIdx(const Dune::cpgrid::CpGridData& nei
         faceCoords.insert(Dune::cpgrid::Entity<3>(cellRef, vertexIdx, true).geometry().center());
     }
 
-    /* for (const auto& faceIdx : neighborFaceList) {
+    for (const auto& faceIdx : neighborFaceList) {
         std::set<Dune::FieldVector<double,3>, FieldVectorLess> neighborFaceCoords{};
         for (const auto& vIdx : neighborCellRef.faceToPoint(faceIdx)){
             neighborFaceCoords.insert(Dune::cpgrid::Entity<3>(neighborCellRef, vIdx, true).geometry().center());
@@ -1510,16 +1620,14 @@ std::optional<int> getNeighborCellRefFaceIdx(const Dune::cpgrid::CpGridData& nei
         }
             return faceIdx;
     }
-    return std::nullopt;*/
+    return std::nullopt;
 
 
-    for (const auto& faceIdx : neighborFaceList) {
+    /*  for (const auto& faceIdx : neighborFaceList) {
     std::set<Dune::FieldVector<double,3>, FieldVectorLess> neighborFaceCoords;
 
     for (const auto& vIdx : neighborCellRef.faceToPoint(faceIdx)) {
-        // std::cout<< vIdx << " adding vertex index in neighborList "<< std::endl;
-        neighborFaceCoords.insert(
-            Dune::cpgrid::Entity<3>(neighborCellRef, vIdx, true).geometry().center());
+        neighborFaceCoords.insert(Dune::cpgrid::Entity<3>(neighborCellRef, vIdx, true).geometry().center());
     }
 
     if (faceCoords.size() != neighborFaceCoords.size())
@@ -1547,7 +1655,7 @@ std::optional<int> getNeighborCellRefFaceIdx(const Dune::cpgrid::CpGridData& nei
         return faceIdx;
 }
 
-return std::nullopt;
+return std::nullopt;*/
 }