tnl-run-fsm-eoc-test: Change of the final file for computeError

tnlDirectEikonalMethodsBase_impl.h: 2D method repaired, 3D method ready for debug
tnlDirectEikonalMethod2D_impl.h: Added save of aux-* as results from computation
tnlDirectEikonalMethod3D_impl.h: Implemented all sweaps and saving results from computation

src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnl-run-fsm-eoc-test

@@ -2,7 +2,7 @@
 
 device="host"
 dimensions="2D 3D"
-dimensions="2D"
+#dimensions="3D"
 sizes1D="16 32 64 128 256 512 1024 2048 4096"
 #sizes1D="256"
 sizes2D="16 32 64 128 256 512 1024"
@@ -12,7 +12,7 @@ testFunctions="paraboloid"
 snapshotPeriod=0.1
 finalTime=1.5
 solverName="tnl-direct-eikonal-solver"
-#solverName="gdb --args tnl-direct-eikonal-solver-dbg"
+#solverName="gdb --args tnl-direct-eikonal-solver-dbg --catch-exceptions no"
 #
 
 setupTestFunction()
@@ -59,7 +59,7 @@ setupGrid()
 setInitialCondition()
 {
    testFunction=$1
-   tnl-init --test-function ${testFunction}-sdf \
+   tnl-init --test-function ${testFunction} \
             --output-file initial-u.tnl \
             --amplitude ${amplitude} \
             --wave-length ${waveLength} \
@@ -78,7 +78,7 @@ setInitialCondition()
             --radius ${radius}
    
     tnl-init --test-function ${testFunction}-sdf \
-            --output-file final-u.tnl \
+            --output-file exact-u.tnl \
             --amplitude ${amplitude} \
             --wave-length ${waveLength} \
             --wave-length-x ${waveLengthX} \
@@ -108,7 +108,6 @@ solve()
                  --time-step 10 \
                  --time-step-order 2 \
                  --discrete-solver ${discreteSolver} \
-                 --merson-adaptivity 1.0e-5 \
                  --min-iterations 20 \
                  --convergence-residue 1.0e-12 \
                  --snapshot-period ${snapshotPeriod} \
@@ -118,7 +117,7 @@ solve()
 computeError()
 {
    tnl-diff --mesh mesh.tnl \
-            --input-files final-u.tnl u-*.tnl \
+            --input-files aux-final.tnl exact-u.tnl \
             --mode sequence \
             --snapshot-period ${snapshotPeriod} \
             --output-file errors.txt \

src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalMethodsBase.h

@@ -99,6 +99,6 @@ template < typename T1, typename T2 >
 T1 meet2DCondition( T1 a, T1 b, const T2 ha, const T2 hb, const T1 value, double v = 1);
 
 template < typename T1 >
-T1 sortMinims( T1 a, T1 b, T1 c, T1 ha, T1 hb, T1 hc);
+void sortMinims( T1 pom[] );
 
 #include "tnlDirectEikonalMethodsBase_impl.h"

src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalMethodsBase_impl.h

@@ -86,12 +86,13 @@ initInterface( const MeshFunctionType& input,
         for( cell.getCoordinates().x() = 0;
              cell.getCoordinates().x() < mesh.getDimensions().x();
              cell.getCoordinates().x() ++ )
-        {
-            output[ cell.getIndex() ] =
+            {
+                cell.refresh();
+                output[ cell.getIndex() ] =
                 input( cell ) >= 0 ? TypeInfo< RealType >::getMaxValue() :
                                    - TypeInfo< RealType >::getMaxValue();
-            interfaceMap[ cell.getIndex() ] = false;
-        }
+                interfaceMap[ cell.getIndex() ] = false;
+            }
    
    const RealType& hx = mesh.getSpaceSteps().x();
    const RealType& hy = mesh.getSpaceSteps().y();     
@@ -117,15 +118,15 @@ initInterface( const MeshFunctionType& input,
                     pom = ( hy * c )/( c - input[ n ]);
                     if( output[ cell.getIndex() ] > pom ) 
                         output[ cell.getIndex() ] = pom;
-                    
-                    output[ n ] = ( hy * c )/( c - input[ n ]) - hy;
+                    if( output[ n ] < pom-hx )
+                        output[ n ] = pom - hy; //( hy * c )/( c - input[ n ]) - hy;
                 }else
                 {
                     pom = - ( hy * c )/( c - input[ n ]);
                     if( output[ cell.getIndex() ] < pom )
                         output[ cell.getIndex() ] = pom;
-                    
-                    output[ n ] = hy - ( hy * c )/( c - input[ n ]);
+                    if( output[ n ] > hy + pom )
+                        output[ n ] = hy + pom; //hy - ( hy * c )/( c - input[ n ]);
                 }
                 interfaceMap[ cell.getIndex() ] = true;
                 interfaceMap[ n ] = true;
@@ -139,7 +140,7 @@ initInterface( const MeshFunctionType& input,
                         output[ cell.getIndex() ] = pom;
                    
                     pom = pom - hx; //output[ e ] = (hx * c)/( c - input[ e ]) - hx;
-                    if( output[ e ] >= pom )
+                    if( output[ e ] < pom )
                         output[ e ] = pom;                         
                 }else
                 {
@@ -148,7 +149,7 @@ initInterface( const MeshFunctionType& input,
                         output[ cell.getIndex() ] = pom;
                     
                     pom = pom + hx; //output[ e ] = hx - (hx * c)/( c - input[ e ]);
-                    if( output[ e ] < pom )
+                    if( output[ e ] > pom )
                         output[ e ] = pom;
                 }
                 interfaceMap[ cell.getIndex() ] = true;
@@ -157,7 +158,7 @@ initInterface( const MeshFunctionType& input,
          }
          /*output[ cell.getIndex() ] =
             c > 0 ? TypeInfo< RealType >::getMaxValue() :
-                   -TypeInfo< RealType >::getMaxValue();*/  
+                   -TypeInfo< RealType >::getMaxValue(); */
          //interfaceMap[ cell.getIndex() ] = false; //is on line 90
       }
 }
@@ -178,7 +179,7 @@ updateCell( MeshFunctionType& u,
    const RealType& hx = mesh.getSpaceSteps().x();
    const RealType& hy = mesh.getSpaceSteps().y();
    const RealType value = u( cell );
-   Real a, b, tmp;
+   Real a, b, tmp = TypeInfo< Real >::getMaxValue();
    
    if( cell.getCoordinates().x() == 0 )
       a = u[ neighborEntities.template getEntityIndex< 1,  0 >() ];
@@ -211,10 +212,10 @@ updateCell( MeshFunctionType& u,
                                  a + TNL::sign( value ) * hx;
    }*/
    tmp = meet2DCondition( a, b, hx, hy, value, v);
-   if( tmp == 0 )
+   if( tmp == TypeInfo< Real >::getMaxValue() )
       tmp = 
-          fabs( a ) >= fabs( b ) ? fabs( b ) + TNL::sign( value ) * hy :
-                    fabs( a ) + TNL::sign( value ) * hx;
+          fabs( a ) >= fabs( b ) ? b + TNL::sign( value ) * hy :
+                                   a + TNL::sign( value ) * hx;
             
    
    u[ cell.getIndex() ] = argAbsMin( value, tmp );
@@ -242,6 +243,7 @@ initInterface( const MeshFunctionType& input,
                  cell.getCoordinates().x() < mesh.getDimensions().x();
                  cell.getCoordinates().x() ++ )
             {
+                cell.refresh();
                 output[ cell.getIndex() ] =
                 input( cell ) >= 0 ? TypeInfo< RealType >::getMaxValue() :
                                    - TypeInfo< RealType >::getMaxValue();
@@ -288,16 +290,19 @@ initInterface( const MeshFunctionType& input,
                    {
                    pom = ( hy * c )/( c - input[ n ]);
                    if( output[ cell.getIndex() ] > pom ) 
-                   output[ cell.getIndex() ] = pom;
+                       output[ cell.getIndex() ] = pom;
                     
-                   output[ n ] = ( hy * c )/( c - input[ n ]) - hy;
+                   if ( output[ n ] < pom - hy)
+                        output[ n ] = pom - hy; // ( hy * c )/( c - input[ n ]) - hy;
+                   
                    }else
                    {
                      pom = - ( hy * c )/( c - input[ n ]);
                      if( output[ cell.getIndex() ] < pom )
                          output[ cell.getIndex() ] = pom;
-                    
-                     output[ n ] = hy - ( hy * c )/( c - input[ n ]);
+                     if( output[ n ] > hy + pom )
+                         output[ n ] = hy + pom; //hy - ( hy * c )/( c - input[ n ]);
+                   
                    }
                interfaceMap[ cell.getIndex() ] = true;
                interfaceMap[ n ] = true;
@@ -308,11 +313,12 @@ initInterface( const MeshFunctionType& input,
                    {
                        pom = ( hx * c )/( c - input[ e ]);
                        if( output[ cell.getIndex() ] > pom )
-                       output[ cell.getIndex() ] = pom;
+                           output[ cell.getIndex() ] = pom;
                    
                        pom = pom - hx; //output[ e ] = (hx * c)/( c - input[ e ]) - hx;
-                       if( output[ e ] >= pom )
-                           output[ e ] = pom;                         
+                       if( output[ e ] < pom )
+                           output[ e ] = pom;      
+                       
                    }else
                    {
                        pom = - (hx * c)/( c - input[ e ]);
@@ -320,7 +326,7 @@ initInterface( const MeshFunctionType& input,
                            output[ cell.getIndex() ] = pom;
                     
                        pom = pom + hx; //output[ e ] = hx - (hx * c)/( c - input[ e ]);
-                       if( output[ e ] < pom )
+                       if( output[ e ] > pom )
                            output[ e ] = pom;
                    }
                interfaceMap[ cell.getIndex() ] = true;
@@ -335,8 +341,9 @@ initInterface( const MeshFunctionType& input,
                            output[ cell.getIndex() ] = pom;
                    
                        pom = pom - hz; //output[ e ] = (hx * c)/( c - input[ e ]) - hx;
-                       if( output[ t ] >= pom )
-                           output[ t ] = pom;                         
+                       if( output[ t ] < pom )
+                           output[ t ] = pom; 
+                       
                    }else
                    {
                        pom = - (hz * c)/( c - input[ t ]);
@@ -344,8 +351,9 @@ initInterface( const MeshFunctionType& input,
                            output[ cell.getIndex() ] = pom;
                     
                        pom = pom + hz; //output[ e ] = hx - (hx * c)/( c - input[ e ]);
-                       if( output[ t ] < pom )
+                       if( output[ t ] > pom )
                            output[ t ] = pom;
+                       
                    }
                interfaceMap[ cell.getIndex() ] = true;
                interfaceMap[ t ] = true;
@@ -368,24 +376,15 @@ updateCell( MeshFunctionType& u,
             const MeshEntity& cell, 
             double v )
 {
-   const auto& neighborEntities = cell.template getNeighborEntities< 2 >();
+   const auto& neighborEntities = cell.template getNeighborEntities< 3 >();
    const MeshType& mesh = cell.getMesh();
   
    const RealType& hx = mesh.getSpaceSteps().x();
    const RealType& hy = mesh.getSpaceSteps().y();
    const RealType& hz = mesh.getSpaceSteps().z();
    const RealType value = u( cell );
-   Real a, b, c, tmp = 0;
+   Real a, b, c, tmp = TypeInfo< Real >::getMaxValue();
    
-   if( cell.getCoordinates().x() == 0 )
-      a = u[ neighborEntities.template getEntityIndex< 1,  0 >() ];
-   else if( cell.getCoordinates().x() == mesh.getDimensions().x() - 1 )
-      a = u[ neighborEntities.template getEntityIndex< -1,  0 >() ];
-   else
-   {
-      a = TNL::argAbsMin( u[ neighborEntities.template getEntityIndex< -1,  0 >() ],
-                     u[ neighborEntities.template getEntityIndex<  1,  0 >() ] );
-   }
    
    if( cell.getCoordinates().x() == 0 )
       a = u[ neighborEntities.template getEntityIndex< 1, 0, 0 >() ];
@@ -396,17 +395,17 @@ updateCell( MeshFunctionType& u,
       a = TNL::argAbsMin( u[ neighborEntities.template getEntityIndex< -1, 0, 0 >() ],
                         u[ neighborEntities.template getEntityIndex< 1, 0, 0 >() ] );
    }
-   if( cell.getCoordinates().x() == 0 )
+   if( cell.getCoordinates().y() == 0 )
       b = u[ neighborEntities.template getEntityIndex< 0, 1, 0 >() ];
-   else if( cell.getCoordinates().x() == mesh.getDimensions().x() - 1 )
+   else if( cell.getCoordinates().y() == mesh.getDimensions().y() - 1 )
       b = u[ neighborEntities.template getEntityIndex< 0, -1, 0 >() ];
    else
    {
       b = TNL::argAbsMin( u[ neighborEntities.template getEntityIndex< 0, -1, 0 >() ],
                         u[ neighborEntities.template getEntityIndex< 0, 1, 0 >() ] );
-   }if( cell.getCoordinates().x() == 0 )
+   }if( cell.getCoordinates().z() == 0 )
       c = u[ neighborEntities.template getEntityIndex< 0, 0, 1 >() ];
-   else if( cell.getCoordinates().x() == mesh.getDimensions().x() - 1 )
+   else if( cell.getCoordinates().z() == mesh.getDimensions().z() - 1 )
       c = u[ neighborEntities.template getEntityIndex< 0, 0, -1 >() ];
    else
    {
@@ -447,11 +446,12 @@ updateCell( MeshFunctionType& u,
                 sign( value ) * hy * hz * sqrt( ( hy * hy + hz * hz )/v - 
                 ( b - c ) * ( b - c ) ) )/( hy * hy + hz * hz );
        }*/
-       Real pom = sortMinims( a, b, c, (Real)hx, (Real)hy, (Real)hz);
-       tmp = meet2DCondition( pom[0], pom[1], pom[3], pom[4], value, v);
+       Real pom[6] = { a, b, c, (Real)hx, (Real)hy, (Real)hz};
+       sortMinims( pom );
+       tmp = meet2DCondition( pom[0], pom[1], pom[2], pom[3], value, v);
        
-       if( tmp == 0 )
-           tmp = pom[ 0 ] + TNL::sign( value ) * pom[ 3 ]; 
+       if( tmp == TypeInfo< Real >::getMaxValue() )
+           tmp = pom[ 0 ] + TNL::sign( value ) * pom[ 2 ]; 
         
    }
    else
@@ -477,33 +477,64 @@ T1 meet2DCondition( T1 a, T1 b, const T2 ha, const T2 hb, const T1 value, double
             ( a - b ) * ( a - b ) ) )/( ha * ha + hb * hb );
    }
    else
-       tmp = 0;
+   {
+       tmp = TypeInfo< T1 >::getMaxValue();
+   }
    
    return tmp;
 }
 
 template < typename T1 >
-T1 sortMinims( T1 a, T1 b, T1 c, T1 ha, T1 hb, T1 hc)
+void sortMinims( T1 pom[])
+{
+    T1 tmp[6]; 
+    if( pom[0] <= pom[1] && pom[1] <= pom[2]){
+        tmp[0] = pom[0]; tmp[1] = pom[1]; tmp[2] = pom[3]; tmp[3] = pom[4];
+    }
+    else if( pom[0] <= pom[2] && pom[2] <= pom[1] ){
+        tmp[0] = pom[0]; tmp[1] = pom[2]; tmp[2] = pom[3]; tmp[3] = pom[5];
+    }
+    else if( pom[1]<=pom[0] && pom[0] <= pom[2] ){
+        tmp[0] = pom[1]; tmp[1] = pom[0]; tmp[2] = pom[4]; tmp[3] = pom[3];
+    }
+    else if( pom[1] <= pom[2] && pom[2] <= pom[0] ){
+        tmp[0] = pom[1]; tmp[1] = pom[2]; tmp[2] = pom[4]; tmp[3] = pom[5];
+    }
+    else if( pom[2] <= pom[0] && pom[0] <= pom[1] ){
+        tmp[0] = pom[2]; tmp[1] = pom[0]; tmp[2] = pom[5]; tmp[3] = pom[3];
+    }
+    else if( pom[2] <= pom[1] && pom[1] <= pom[0] ){
+        tmp[0] = pom[2]; tmp[1] = pom[1]; tmp[2] = pom[5]; tmp[3] = pom[4];
+    }
+    for( int i = 0; i < 6; i++ )
+        pom[ i ] = tmp[ i ];
+        
+}
+
+/*template < typename T1 >
+T1* sortMinims( T1 pom[], T1 a, T1 b, T1 c, T1 ha, T1 hb, T1 hc)
 {
-    T1 tmp[6];
+    T1 tmp[6] = pom; 
     if( a <= b && b <= c){
-        tmp = {a,b,c,ha,hb,hc};
+        tmp[0] = a; tmp[1] = b; tmp[2] = ha; tmp[3] = hb;
     }
     else if( a <= c && c <= b ){
-        tmp = {a,c,b,ha,hc,hb};
+        tmp[0] = a; tmp[1] = c; tmp[2] = ha; tmp[3] = hc;
     }
     else if( b<=a && a<=c ){
-        tmp = {b,a,c,hb,ha,hc};
+        tmp[0] = b; tmp[1] = a; tmp[2] = hb; tmp[3] = ha;
     }
     else if( b <= c && c <= a ){
-        tmp = {b,c,a,hb,hc,ha};
+        tmp[0] = b; tmp[1] = c; tmp[2] = hb; tmp[3] = hc;
     }
     else if( c <= a && a <= b ){
-        tmp = {c,a,b,hc,ha,hb};
+        tmp[0] = c; tmp[1] = a; tmp[2] = hc; tmp[3] = ha;
     }
     else if( c <=b && b <= a ){
-        tmp = {c,b,a,hc,hb,ha};
+        tmp[0] = c; tmp[1] = b; tmp[2] = hc; tmp[3] = hb;
     }
+    for( int i = 0; i < 6; i++ )
+        pom[ i ] = tmp[ i ];
         
-    return tmp;
-}
\ No newline at end of file
+    return pom;
+}*/
\ No newline at end of file

src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod2D_impl.h

@@ -64,7 +64,7 @@ solve( const MeshPointer& mesh,
    interfaceMap.setMesh( mesh );
    std::cout << "Initiating the interface cells ..." << std::endl;
    BaseType::initInterface( u, aux, interfaceMap );
-   //aux.save( "aux-ini.tnl" );
+   aux.save( "aux-ini.tnl" );
 
    typename MeshType::Cell cell( *mesh );
    
@@ -85,7 +85,7 @@ solve( const MeshPointer& mesh,
                   this->updateCell( aux, cell );
             }
       }
-      //aux.save( "aux-1.tnl" );
+      aux.save( "aux-1.tnl" );
 
       for( cell.getCoordinates().y() = 0;
            cell.getCoordinates().y() < mesh->getDimensions().y();
@@ -101,7 +101,7 @@ solve( const MeshPointer& mesh,
                   this->updateCell( aux, cell );
             }
       }
-      //aux.save( "aux-2.tnl" );
+      aux.save( "aux-2.tnl" );
 
       for( cell.getCoordinates().y() = mesh->getDimensions().y() - 1;
            cell.getCoordinates().y() >= 0 ;
@@ -117,7 +117,7 @@ solve( const MeshPointer& mesh,
                   this->updateCell( aux, cell );
             }
          }
-      //aux.save( "aux-3.tnl" );
+      aux.save( "aux-3.tnl" );
 
 
       for( cell.getCoordinates().y() = mesh->getDimensions().y() - 1;
@@ -134,8 +134,9 @@ solve( const MeshPointer& mesh,
                   this->updateCell( aux, cell );
             }
          }   
-      //aux.save( "aux-4.tnl" );
+      aux.save( "aux-4.tnl" );
       iteration++;
    }
+   aux.save("aux-final.tnl");
 }
 

src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod3D_impl.h

@@ -65,5 +65,180 @@ solve( const MeshPointer& mesh,
    std::cout << "Initiating the interface cells ..." << std::endl;
    BaseType::initInterface( u, aux, interfaceMap );
    aux.save( "aux-ini.tnl" );   
+   
+   typename MeshType::Cell cell( *mesh );
+   
+   IndexType iteration( 0 );
+   while( iteration < this->maxIterations )
+   {
+      for( cell.getCoordinates().z() = 0;
+           cell.getCoordinates().z() < mesh->getDimensions().z();
+           cell.getCoordinates().z()++ )
+      {
+         for( cell.getCoordinates().y() = 0;
+              cell.getCoordinates().y() < mesh->getDimensions().y();
+              cell.getCoordinates().y()++ )
+         {
+            for( cell.getCoordinates().x() = 0;
+                 cell.getCoordinates().x() < mesh->getDimensions().x();
+                 cell.getCoordinates().x()++ )
+            {
+               cell.refresh();
+               if( ! interfaceMap( cell ) )
+                  this->updateCell( aux, cell );
+            }
+         }
+      }
+      aux.save( "aux-1.tnl" );
+
+      for( cell.getCoordinates().z() = 0;
+           cell.getCoordinates().z() < mesh->getDimensions().z();
+           cell.getCoordinates().z()++ )
+      {
+         for( cell.getCoordinates().y() = 0;
+              cell.getCoordinates().y() < mesh->getDimensions().y();
+              cell.getCoordinates().y()++ )
+         {
+            for( cell.getCoordinates().x() = mesh->getDimensions().x() - 1;
+                 cell.getCoordinates().x() >= 0 ;
+                 cell.getCoordinates().x()-- )		
+            {
+               //std::cerr << "2 -> ";
+               cell.refresh();
+               if( ! interfaceMap( cell ) )            
+                  this->updateCell( aux, cell );
+            }
+         }
+      }
+      aux.save( "aux-2.tnl" );
+      for( cell.getCoordinates().z() = 0;
+           cell.getCoordinates().z() < mesh->getDimensions().z();
+           cell.getCoordinates().z()++ )
+      {
+         for( cell.getCoordinates().y() = mesh->getDimensions().y() - 1;
+              cell.getCoordinates().y() >= 0 ;
+              cell.getCoordinates().y()-- )
+         {
+            for( cell.getCoordinates().x() = 0;
+                 cell.getCoordinates().x() < mesh->getDimensions().x();
+                 cell.getCoordinates().x()++ )
+            {
+               //std::cerr << "3 -> ";
+               cell.refresh();
+               if( ! interfaceMap( cell ) )            
+                  this->updateCell( aux, cell );
+            }
+         }
+      }
+      aux.save( "aux-3.tnl" );
+      
+      for( cell.getCoordinates().z() = 0;
+           cell.getCoordinates().z() < mesh->getDimensions().z();
+           cell.getCoordinates().z()++ )
+      {
+         for( cell.getCoordinates().y() = mesh->getDimensions().y() - 1;
+              cell.getCoordinates().y() >= 0;
+              cell.getCoordinates().y()-- )
+         {
+            for( cell.getCoordinates().x() = mesh->getDimensions().x() - 1;
+                 cell.getCoordinates().x() >= 0 ;
+                 cell.getCoordinates().x()-- )		
+            {
+               //std::cerr << "4 -> ";
+               cell.refresh();
+               if( ! interfaceMap( cell ) )            
+                  this->updateCell( aux, cell );
+            }
+         }
+      }     
+      aux.save( "aux-4.tnl" );
+      
+      for( cell.getCoordinates().z() = mesh->getDimensions().z() - 1;
+           cell.getCoordinates().z() >= 0;
+           cell.getCoordinates().z()-- )
+      {
+         for( cell.getCoordinates().y() = 0;
+              cell.getCoordinates().y() < mesh->getDimensions().y();
+              cell.getCoordinates().y()++ )
+         {
+            for( cell.getCoordinates().x() = 0;
+                 cell.getCoordinates().x() < mesh->getDimensions().x();
+                 cell.getCoordinates().x()++ )
+            {
+               //std::cerr << "5 -> ";
+               cell.refresh();
+               if( ! interfaceMap( cell ) )
+                  this->updateCell( aux, cell );
+            }
+         }
+      }
+      aux.save( "aux-5.tnl" );
+
+      for( cell.getCoordinates().z() = mesh->getDimensions().z() - 1;
+           cell.getCoordinates().z() >= 0;
+           cell.getCoordinates().z()-- )
+      {
+         for( cell.getCoordinates().y() = 0;
+              cell.getCoordinates().y() < mesh->getDimensions().y();
+              cell.getCoordinates().y()++ )
+         {
+            for( cell.getCoordinates().x() = mesh->getDimensions().x() - 1;
+                 cell.getCoordinates().x() >= 0 ;
+                 cell.getCoordinates().x()-- )		
+            {
+               //std::cerr << "6 -> ";
+               cell.refresh();
+               if( ! interfaceMap( cell ) )            
+                  this->updateCell( aux, cell );
+            }
+         }
+      }
+      aux.save( "aux-6.tnl" );
+      
+      for( cell.getCoordinates().z() = mesh->getDimensions().z() - 1;
+           cell.getCoordinates().z() >= 0;
+           cell.getCoordinates().z()-- )
+      {
+         for( cell.getCoordinates().y() = mesh->getDimensions().y() - 1;
+              cell.getCoordinates().y() >= 0 ;
+              cell.getCoordinates().y()-- )
+         {
+            for( cell.getCoordinates().x() = 0;
+                 cell.getCoordinates().x() < mesh->getDimensions().x();
+                 cell.getCoordinates().x()++ )
+            {
+               //std::cerr << "7 -> ";
+               cell.refresh();
+               if( ! interfaceMap( cell ) )            
+                  this->updateCell( aux, cell );
+            }
+         }
+      }
+      aux.save( "aux-7.tnl" );
+
+      for( cell.getCoordinates().z() = mesh->getDimensions().z() - 1;
+           cell.getCoordinates().z() >= 0;
+           cell.getCoordinates().z()-- )
+      {
+         for( cell.getCoordinates().y() = mesh->getDimensions().y() - 1;
+              cell.getCoordinates().y() >= 0;
+              cell.getCoordinates().y()-- )
+         {
+            for( cell.getCoordinates().x() = mesh->getDimensions().x() - 1;
+                 cell.getCoordinates().x() >= 0 ;
+                 cell.getCoordinates().x()-- )		
+            {
+               //std::cerr << "8 -> ";
+               cell.refresh();
+               if( ! interfaceMap( cell ) )            
+                  this->updateCell( aux, cell );
+            }
+         }
+      }
+      aux.save( "aux-8.tnl" );
+      iteration++;
+      
+   }
+   aux.save("aux-final.tnl");
 }