[!52] Fix periodic boundaries.

Merge branch 'bugfix/periodic-boundaries' into 'master'

ref:extensions/dune-alugrid This adds a suggested fix for issue [#43].

See merge request [extensions/dune-alugrid!52]

  [#43]: gitlab.dune-project.org/NoneNone/issues/43
  [extensions/dune-alugrid!52]: gitlab.dune-project.org/extensions/dune-alugrid/merge_requests/52

dune/alugrid/3d/gridfactory.cc

@@ -975,6 +975,17 @@ namespace Dune
       const FaceMapIterator fend = faceMap.end();
       for( FaceMapIterator fit = faceMap.begin(); fit != fend; ++fit )
       {
+        //for dimension == 2 we do not want to search
+        // the artificially introduced faces
+        if(dimension == 2)
+        {
+          if(elementType == hexa)
+            if(fit->second.second > 3)
+              continue;
+          if(elementType == tetra)
+            if(fit->second.second > 2)
+              continue;
+        }
         FaceType key2;
         generateFace( fit->second, key2 );
 
@@ -1010,6 +1021,9 @@ namespace Dune
   {
     typedef typename FaceMap::iterator FaceIterator;
     FaceMap faceMap;
+    // list of face that should be removed
+    std::vector< FaceType > toBeDeletedFaces;
+    toBeDeletedFaces.reserve( faceMap.size() / 10 + 1);
 
     const unsigned int numElements = elements_.size();
     for( unsigned int n = 0; n < numElements; ++n )
@@ -1051,15 +1065,35 @@ namespace Dune
       }
 
       reinsertBoundary( faceMap, pos, bndIt->second );
-      faceMap.erase( pos );
+      toBeDeletedFaces.push_back( key );
     }
 
     //the search for the periodic neighbour also deletes the
     //found faces from the boundaryIds_ - thus it has to be done
     //after the recreation of the Ids_, because of correctElementOrientation
-    for(auto it = faceMap.begin(); it!=faceMap.end(); ++it)
+    if( !faceTransformations_.empty() )
+    {
+      for(auto it = faceMap.begin(); it!=faceMap.end(); ++it)
+      {
+        //for dimension == 2 we do not want to search
+        // the artificially introduced faces
+        if(dimension == 2)
+        {
+          if(elementType == hexa)
+            if(it->second.second > 3)
+              continue;
+          if(elementType == tetra)
+            if(it->second.second > 2)
+              continue;
+        }
+        searchPeriodicNeighbor( faceMap, it, defaultId );
+      }
+    }
+
+    // erase faces that are boundries
+    for( const auto& key : toBeDeletedFaces )
     {
-      searchPeriodicNeighbor( faceMap, it, defaultId );
+      faceMap.erase( key );
     }
 
     // communicate unidentified boundaries and find process borders)

dune/alugrid/test/dgf/CMakeLists.txt

@@ -4,6 +4,7 @@ set(DGFS
   cube-testgrid-2-2.dgf
   cube-testgrid-2-3.dgf
   grid2d_str1d.dgf
+  periodic.dgf
   simplex-testgrid-2-2.dgf
   simplex-testgrid-2-3.dgf
   simplex-testgrid-3-3.dgf

dune/alugrid/test/dgf/periodic.dgf

+DGF
+
+INTERVAL
+0 0 0            % lower left corner
+1 1 1            % upper right corner
+8 8 8            % number of cells in each direction
+#
+
+PERIODICFACETRANSFORMATION
+1 0 0, 0 1 0, 0 0 1 + 1 0 0  % make periodic in x-direction
+1 0 0, 0 1 0, 0 0 1 + 0 1 0  % make periodic in y-direction
+#
+
+BOUNDARYDOMAIN
+DEFAULT 1
+#
+
+GRIDPARAMETER
+OVERLAP 0
+#
+
+#

dune/alugrid/test/test-alugrid.cc

@@ -476,7 +476,20 @@ void checkGrid( GridType& grid )
   }
 }
 
+template <class GridType>
+void checkForPeriodicBoundaries( GridType& grid )
+{
+  for (const auto& element : elements(grid.leafGridView()))
+  {
+    for (const auto& intersection : intersections(grid.leafGridView(), element))
+    {
+      if (intersection.neighbor() && intersection.boundary())
+        return;
+    }
+  }
 
+  DUNE_THROW( Dune::InvalidStateException, "No periodic boundaries found!" );
+}
 
 template <class GridType>
 void checkALUSerial(GridType & grid, int mxl = 2)
@@ -853,25 +866,42 @@ int main (int argc , char **argv) {
           filename = "./dgf/simplex-testgrid-3-3.dgf";
 
         typedef Dune::ALUGrid< 3, 3, Dune::cube, Dune::nonconforming > GridType;
-        Dune::GridPtr< GridType > gridPtr( filename );
-        gridPtr.loadBalance();
-        GridType & grid = *gridPtr;
-
-        checkCapabilities< false >( grid );
-
         {
-          std::cout << "Check serial grid" << std::endl;
-          checkALUSerial(grid,
-                         (mysize == 1) ? 1 : 0 );
+          Dune::GridPtr< GridType > gridPtr( filename );
+          gridPtr.loadBalance();
+          GridType & grid = *gridPtr;
+
+          checkCapabilities< false >( grid );
+
+          {
+            std::cout << "Check serial grid" << std::endl;
+            checkALUSerial(grid,
+                           (mysize == 1) ? 1 : 0 );
+          }
+
+          // perform parallel check only when more then one proc
+          if(mysize > 1)
+          {
+            if (myrank == 0) std::cout << "Check conform grid" << std::endl;
+            checkALUParallel(grid,1,0);
+            if (myrank == 0) std::cout << "Check non-conform grid" << std::endl;
+            checkALUParallel(grid,0,2);
+          }
         }
 
-        // perform parallel check only when more then one proc
-        if(mysize > 1)
+        // check periodic capabilities
         {
-          if (myrank == 0) std::cout << "Check conform grid" << std::endl;
-          checkALUParallel(grid,1,0);
-          if (myrank == 0) std::cout << "Check non-conform grid" << std::endl;
-          checkALUParallel(grid,0,2);
+          filename = "./dgf/periodic.dgf";
+          Dune::GridPtr< GridType > gridPtr( filename );
+          gridPtr.loadBalance();
+          GridType & grid = *gridPtr;
+
+          {
+            std::cout << "Check periodic grid" << std::endl;
+            checkALUSerial(grid,
+                           (mysize == 1) ? 1 : 0 );
+            checkForPeriodicBoundaries( grid );
+          }
         }
       }