@@ -877,6 +877,34 @@ void addFaces(int
877877 }
878878}
879879
880+ std::vector<int > modifyFaceToPoint (const std::vector<Dune::FieldVector<double ,3 >>& newFaceToCoord,
881+ const GridModificationMapping& modificationMap,
882+ const BoundaryFaceInfo& boundaryFaceInfo)
883+ {
884+ std::vector<int > faceToPoint;
885+
886+ for (const auto & vertex : newFaceToCoord) {
887+ int vertexIdx = -1 ;
888+ // Check if the vertex is new (e.g., added when modifying the grid)
889+ auto newIt = modificationMap.newGridVertexCoordinates_to_newGridVertexIdx .find (vertex);
890+ if (newIt != modificationMap.newGridVertexCoordinates_to_newGridVertexIdx .end ()) {
891+ vertexIdx = newIt->second ;
892+ }
893+ else { // Vertex already existed in the original grid ()
894+ auto oldIt = boundaryFaceInfo.boundaryVertex_to_vertexIdx .find (vertex);
895+ assert (oldIt != boundaryFaceInfo.boundaryVertex_to_vertexIdx .end ());
896+
897+ vertexIdx = oldIt->second ;
898+ }
899+
900+ assert (vertexIdx >= 0 );
901+ faceToPoint.push_back (vertexIdx);
902+ }
903+
904+ return faceToPoint;
905+ }
906+
907+
880908void populateFaceData (int oldGridNumFaces,
881909 const Dune::cpgrid::CpGridData& oldGrid,
882910 const std::map<int ,BoundaryFaceInfo>& boundaryFacesInfo,
@@ -904,7 +932,7 @@ void populateFaceData(int oldGridNumFaces,
904932
905933 modificationMap.newGridFaceIdx_to_oldGridFaceIdx .resize (upperBoundFaceSize, invalidIdx);
906934
907- // intentional copy
935+ // intentional copy
908936 const auto hasParent = newGridFaceIdx_to_parentFaceIdx;
909937 newGridFaceIdx_to_parentFaceIdx.resize (upperBoundFaceSize, invalidIdx);
910938
@@ -914,107 +942,73 @@ void populateFaceData(int oldGridNumFaces,
914942 int newGrid_num_points = 0 ;
915943 int newGrid_face_count = 0 ;
916944
945+ auto registerFace = [&](int oldFaceIdx, int newFaceIdx, int parentFaceIdx,
946+ auto && geometry, auto && normal, auto tag)
947+ {
948+ modificationMap.oldGridFaceIdx_to_newGridFaceIdxList [oldFaceIdx].push_back (newFaceIdx);
949+ modificationMap.newGridFaceIdx_to_oldGridFaceIdx [newFaceIdx] = oldFaceIdx;
950+
951+ if (parentFaceIdx != -1 ) {
952+ newGridFaceIdx_to_parentFaceIdx[newFaceIdx] = parentFaceIdx;
953+ }
954+
955+ newGridFaces[newFaceIdx] = geometry;
956+ mutable_newGridFaceTags[newFaceIdx] = tag;
957+ mutable_newGridFaceNormals[newFaceIdx] = normal;
958+
959+ const auto & faceToCell = oldGrid.faceToCell (oldFaceIdx);
960+ newGridGeomData.face_to_cell .appendRow (faceToCell.begin (), faceToCell.end ());
961+
962+ ++newGrid_face_count;
963+ };
964+
917965 for (int i = 0 ; i < oldGrid.numFaces (); ++i) {
918966
919967 std::vector<int > newFaceIndices{};
920-
968+
969+ const auto face = Dune::cpgrid::EntityRep<1 >(i, true );
921970 int parentFaceIdx = hasParent[i];
922971 auto itP = boundaryFacesInfo.find (parentFaceIdx);
923972
924973 if ((parentFaceIdx == -1 ) || (itP == boundaryFacesInfo.end ())) {
925-
926- const auto face = Dune::cpgrid::EntityRep<1 >(i, true );
927-
928- modificationMap.oldGridFaceIdx_to_newGridFaceIdxList [i].push_back (newGrid_face_count);
929- modificationMap.newGridFaceIdx_to_oldGridFaceIdx [newGrid_face_count] = i;
930-
931- if (parentFaceIdx!= -1 ) {
932- newGridFaceIdx_to_parentFaceIdx[newGrid_face_count] = parentFaceIdx;
933- }
934-
935- newGridFaces[newGrid_face_count] = (*oldGrid.getGeometry ().geomVector (std::integral_constant<int ,1 >()))[face];
936- mutable_newGridFaceTags[newGrid_face_count] = oldGrid.faceTag (i);
937- mutable_newGridFaceNormals[newGrid_face_count] = oldGrid.faceNormals (i);
938974
939- std::vector<int > faceToPoint{};
940- faceToPoint.reserve (oldGrid.faceToPoint (i).size ());
941- for (const auto & vertexIdx : oldGrid.faceToPoint (i)) {
942- faceToPoint.push_back (vertexIdx);
943- }
975+ const auto faceToPoint = oldGrid.faceToPoint (i);
976+ aux_newGrid_face_to_point[newGrid_face_count] = std::vector<int >(faceToPoint.begin (), faceToPoint.end ());
944977 newGrid_num_points += faceToPoint.size ();
945- aux_newGrid_face_to_point[newGrid_face_count] = faceToPoint;
946978
947- const auto & faceToCell = oldGrid.faceToCell (i);
948- newGridGeomData.face_to_cell .appendRow (faceToCell.begin (), faceToCell.end ());
949-
950- ++newGrid_face_count;
979+ registerFace (i, newGrid_face_count, parentFaceIdx,
980+ (*oldGrid.getGeometry ().geomVector (std::integral_constant<int ,1 >()))[face],
981+ oldGrid.faceNormals (i), oldGrid.faceTag (i));
951982 }
952983 else {
953984 const auto & boundaryFaceInfo = itP->second ;
954985 const auto & newFacesInfo = boundaryFaceInfo.vanishedCellRefFace_to_newRefinedFaces [i];
955986
956- if ( boundaryFaceInfo.faceFullyContainedInNeighbor [i] || newFacesInfo.empty ()){
957- const auto face = Dune::cpgrid::EntityRep<1 >(i, true );
987+ if ( boundaryFaceInfo.faceFullyContainedInNeighbor [i] || newFacesInfo.empty ()){
958988
959- modificationMap.oldGridFaceIdx_to_newGridFaceIdxList [i].push_back (newGrid_face_count);
960- modificationMap.newGridFaceIdx_to_oldGridFaceIdx [newGrid_face_count] = i;
961- newGridFaceIdx_to_parentFaceIdx[newGrid_face_count] = parentFaceIdx;
962-
963- newGridFaces[newGrid_face_count] = (*oldGrid.getGeometry ().geomVector (std::integral_constant<int ,1 >()))[face];
964- mutable_newGridFaceTags[newGrid_face_count] = oldGrid.faceTag (i);
965- mutable_newGridFaceNormals[newGrid_face_count] = oldGrid.faceNormals (i);
966-
967- std::vector<int > faceToPoint{};
968- faceToPoint.reserve (oldGrid.faceToPoint (i).size ());
969- for (const auto & vertexIdx : oldGrid.faceToPoint (i)) {
970- faceToPoint.push_back (vertexIdx);
971- }
989+ const auto faceToPoint = oldGrid.faceToPoint (i);
990+ aux_newGrid_face_to_point[newGrid_face_count] = std::vector<int >(faceToPoint.begin (), faceToPoint.end ());
972991 newGrid_num_points += faceToPoint.size ();
973- aux_newGrid_face_to_point[newGrid_face_count] = faceToPoint;
974992
975- const auto & faceToCell = oldGrid.faceToCell (i);
976- newGridGeomData.face_to_cell .appendRow (faceToCell.begin (), faceToCell.end ());
977-
978- ++newGrid_face_count;
993+ registerFace (i, newGrid_face_count, parentFaceIdx,
994+ (*oldGrid.getGeometry ().geomVector (std::integral_constant<int ,1 >()))[face],
995+ oldGrid.faceNormals (i), oldGrid.faceTag (i));
979996 }
980997 else {
981998 for (const auto & [grid2_faceIdx, newFaceToCoord] : newFacesInfo) {
982-
999+
9831000 newFaceIndices.push_back (newGrid_face_count);
984- modificationMap.oldGridFaceIdx_to_newGridFaceIdxList [i].push_back (newGrid_face_count);
985- modificationMap.newGridFaceIdx_to_oldGridFaceIdx [newGrid_face_count] = i;
986- newGridFaceIdx_to_parentFaceIdx[newGrid_face_count] = parentFaceIdx;
987-
988- const auto [faceCenter, faceArea, faceNormal] = computeFaceCenterAreaNormal (newFaceToCoord);
989-
990- newGridFaces[newGrid_face_count] = Dune::cpgrid::Geometry<2 ,3 >(faceCenter, faceArea);
991- mutable_newGridFaceTags[newGrid_face_count] = oldGrid.faceTag (i); // shared tag
992- mutable_newGridFaceNormals[newGrid_face_count] = faceNormal;
993-
994- std::vector<int > faceToPoint{};
995- faceToPoint.reserve (newFaceToCoord.size ());
996-
997- for (const auto & vertex : newFaceToCoord) {
998- int vertexIdx = -1 ; // invalid to be rewritten
999- // check if vertex (is new) and has already been stored
1000- auto it = modificationMap.newGridVertexCoordinates_to_newGridVertexIdx .find (vertex);
1001- if (it != modificationMap.newGridVertexCoordinates_to_newGridVertexIdx .end () ) {
1002- vertexIdx = it->second ;
1003- }
1004- else { // otherwise, vertex existed already in gridData1 (before "correction")
1005- auto iit = boundaryFaceInfo.boundaryVertex_to_vertexIdx .find (vertex);
1006- assert ( iit != boundaryFaceInfo.boundaryVertex_to_vertexIdx .end ());
1007- vertexIdx = iit->second ;
1008- }
1009- assert (vertexIdx>=0 );
1010- faceToPoint.push_back (vertexIdx);
1011- }
1012- newGrid_num_points += faceToPoint.size ();
1001+
1002+ const auto faceToPoint = modifyFaceToPoint (newFaceToCoord,
1003+ modificationMap,
1004+ boundaryFaceInfo);
10131005 aux_newGrid_face_to_point[newGrid_face_count] = faceToPoint;
1006+ newGrid_num_points += faceToPoint.size ();
10141007
1015- const auto & face_to_cell = oldGrid.faceToCell (i);
1016- newGridGeomData.face_to_cell .appendRow (face_to_cell.begin (), face_to_cell.end ());
1017- ++newGrid_face_count;
1008+ const auto [faceCenter, faceArea, faceNormal] = computeFaceCenterAreaNormal (newFaceToCoord);
1009+ registerFace (i, newGrid_face_count, parentFaceIdx,
1010+ Dune::cpgrid::Geometry<2 ,3 >(faceCenter, faceArea),
1011+ faceNormal, oldGrid.faceTag (i));
10181012 }
10191013 }
10201014 }
@@ -1041,8 +1035,7 @@ void populateCellData(const Dune::cpgrid::CpGridData& sourceGrid,
10411035
10421036 const auto element = Dune::cpgrid::Entity<0 >(sourceGrid, i, true );
10431037
1044- cells[i] = Dune::cpgrid::Geometry<3 ,3 >(element.geometry ().center (),
1045- element.geometry ().volume (),
1038+ cells[i] = Dune::cpgrid::Geometry<3 ,3 >(element.geometry ().center (), element.geometry ().volume (),
10461039 geomData.geometries .geomVector (std::integral_constant<int ,3 >()),
10471040 indices_storage_ptr);
10481041 }
0 commit comments