[bugfix][periodic] Avoid using YaspGrid for inserting elements in

StructuredGridFactory when periodic transformation is present.

dune/alugrid/common/structuredgridfactory.hh

@@ -258,6 +258,16 @@ namespace Dune
       CollectiveCommunication comm( MPIHelper :: getCommunicator() );
       static_assert( dim == dimworld, "YaspGrid is used for creation of the structured grid which only supports dim == dimworld");
 
+      // if periodic transformations are active we cannot use the YaspGrid
+      // approach to insert the grid cells, otherwise the periodic elements
+      // are not inserted
+      dgf::PeriodicFaceTransformationBlock trafoBlock( input, dimworld );
+      if( trafoBlock.isactive() )
+      {
+        Dune::GridPtr< Grid > grid( input, mpiComm );
+        return SharedPtrType( grid.release() );
+      }
+
       Dune::dgf::IntervalBlock intervalBlock( input );
       if( !intervalBlock.isactive() )
       {

dune/alugrid/impl/serial/gitter_hexa_top.cc

@@ -1422,8 +1422,8 @@ namespace ALUGrid
     return;
   }
 
-  template< class A > void Periodic4Top < A >::refineImmediate (myrule_t r) {
-
+  template< class A > void Periodic4Top < A >::refineImmediate (myrule_t r)
+  {
     // Die Methode wird nur vom restore () und vom refineBalance () auf-
     // gerufen und geht davon aus, dass das betroffene Element noch nicht
     // verfeinert ist -> ist ein Blatt der Hierarchie.
@@ -1475,9 +1475,10 @@ namespace ALUGrid
       // sich direkt auf die Balancierungsregel des entsprechenden Polygonverbinders
       // projezieren l"asst (n"amlich 1:1). Deshalb unterscheidet der Aufruf nicht nach
       // der angeforderten Regel in einer 'case' Anweisung.
+      alugrid_assert( fce == 0 || fce == 1 );
 
       typedef typename myhface4_t::myrule_t myhface4rule_t;
-      int opp = fce == 0 ? 1 : 0;
+      int opp = 1 - fce; // (fce == 0) ? 1 : 0;
       if (myhface4 (opp)->refine (myhface4rule_t (r).rotate (twist (opp)), twist (opp)))
       {
         refineImmediate( r );

examples/CMakeLists.txt

@@ -29,7 +29,7 @@ foreach(variant ball euler transport)
   if("${variant}" STREQUAL "ball")
     target_compile_definitions(main_${variant} PRIVATE "ALUGRID_COUNT_GLOBALCOMM")
   endif()
-  target_compile_definitions(main_${variant} PUBLIC "ALUGRID_CUBE" "GRIDDIM=3" "WORLDDIM=3")
+  target_compile_definitions(main_${variant} PUBLIC "ALUGRID_CUBE" "GRIDDIM=2" "WORLDDIM=2")
 
   foreach(type cube simplex conform)
     string(TOUPPER ${type} type_upper)

examples/dgf/periodic2.dgf

+DGF
+
+Interval
+ 0.0  0.0    % first corner
+ 1.0  1.0    % second corner
+ 8    8      % 8 cells in each direction (x, y)
+#
+
+
+PERIODICFACETRANSFORMATION
+1 0, 0 1 + 1 0
+#
+

examples/dgf/periodic3.dgf

+DGF
+
+Interval
+ 0.0  0.0 0.0   % first corner
+ 1.0  1.0 1.0   % second corner
+ 4    4   4     % 4 cells in each direction (x, y, z)
+#
+
+
+PERIODICFACETRANSFORMATION
+1 0 0, 0 0 1 + 1 0 0
+#
+

examples/problem-ball.hh

@@ -5,6 +5,7 @@
 #include <cassert>
 #include <iostream>
 #include <sstream>
+#include <memory>
 
 #include <dune/common/fvector.hh>
 
@@ -22,21 +23,27 @@ struct BallData
 
   const static int dimDomain = DomainType::dimension;
 
-  explicit BallData (const int problem) : c_(0.5), r0_(0.3), problem_( dimDomain == 3 ? problem : 0)
+  explicit BallData (const int problem)
+    : c_(0.5), r0_(0.3), problem_( problem ) //dimDomain == 3 ? problem : 0)
   {}
 
   //! \copydoc ProblemData::gridFile
   std::string gridFile ( const std::string &path, const int mpiSize ) const
   {
     std::ostringstream dgfFileName;
-    if( problem_ == 1 )
-      dgfFileName << path << "/dgf/cube_hc_512.dgf";
-    else if ( problem_ == 2 )
-      dgfFileName << path << "/dgf/cube_hc_4096.dgf";
-    else if ( problem_ == 3 )
-      dgfFileName << path << "/dgf/cube_hc_32768.dgf";
+    if( dimDomain == 3 && problem_ < 4 )
+    {
+      if( problem_ == 1 )
+        dgfFileName << path << "/dgf/cube_hc_512.dgf";
+      else if ( problem_ == 2 )
+        dgfFileName << path << "/dgf/cube_hc_4096.dgf";
+      else if ( problem_ == 3 )
+        dgfFileName << path << "/dgf/cube_hc_32768.dgf";
+    }
     else if ( problem_ == 4 )
       dgfFileName << path << "/dgf/input" << dimDomain << ".dgf";
+    else if ( problem_ == 5 )
+      dgfFileName << path << "/dgf/periodic" << dimDomain << ".dgf";
     else
       dgfFileName << path << "/dgf/unitcube" << dimDomain << "d.dgf";
     return dgfFileName.str();
@@ -60,8 +67,19 @@ struct BallData
     DomainType xx(x);
     xx -= c_;
     DomainType y(0);
-    y[0] = std::cos(time*2.*M_PI)*r0_;
-    y[1] = std::sin(time*2.*M_PI)*r0_;
+
+    if( problem_ < 5 )
+    {
+      y[0] = std::cos(time*2.*M_PI)*r0_;
+      y[1] = std::sin(time*2.*M_PI)*r0_;
+    }
+    else
+    {
+      // lateral translation of marking zone
+      int inter = int(time);
+      y[0] = time - double(inter) -0.5;
+    }
+
     xx -= y;
     double r = xx.two_norm();
     return ( (r>0.15 && r<0.25)? 1 : 0 );
@@ -106,27 +124,9 @@ struct BallModel : public TransportModel<dimD>
    *  \param problem switch between different data settings
    */
   BallModel( unsigned int problem )
-  : problem_( 0 )
-  {
-    switch( problem )
-    {
-    case 0:
-    case 1:
-    case 2:
-    case 3:
-      problem_ = new BallData< dimDomain >( problem );
-      break;
-
-    default:
-      std::cerr << "ProblemData not defined - using problem 1!" << std::endl;
-      problem_ = new BallData< dimDomain >( problem );
-    }
-  }
-
-  /** \brief destructor */
-  ~BallModel()
+  : problem_( new BallData< dimDomain >( problem ) )
   {
-    delete problem_;
+    std::cout << "Creating BallModel with problem " << problem << std::endl;
   }
 
   /** \copydoc TransportProblem::data */
@@ -158,7 +158,7 @@ struct BallModel : public TransportModel<dimD>
     return indicator( normal,time,xGlobal, uLeft, problem().boundaryValue(xGlobal,time) );
   }
 
-  Problem *problem_;
+  std::unique_ptr< Problem > problem_;
 };
 
 #endif