diff --git a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnl-run-fsm-eoc-test b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnl-run-fsm-eoc-test
index 57a40528b4b33551af8ca5f34f69cb0d00fcc30f..0dc50246f5ced2ec65616f7eecc65ebf860a2cdb 100755
--- a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnl-run-fsm-eoc-test
+++ b/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 \
diff --git a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalMethodsBase.h b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalMethodsBase.h
index 1182ac5745612aafc9e8fa1bdacec9f1883fd681..a4a46396cdc127e0ccbc9420b81c959011d18ad3 100644
--- a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalMethodsBase.h
+++ b/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"
diff --git a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalMethodsBase_impl.h b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalMethodsBase_impl.h
index 461c70f55619556fd364c4e137d861884d67840c..4cb8327a16034a061feeb6a94f6faec15caf61e3 100644
--- a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalMethodsBase_impl.h
+++ b/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
diff --git a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod2D_impl.h b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod2D_impl.h
index 13e5c824912ef387c47ddaca1c5eb98acfe30fe6..bb1d4a018faf433decb517664941c20eda746feb 100644
--- a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod2D_impl.h
+++ b/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");
}
diff --git a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod3D_impl.h b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod3D_impl.h
index 6fad87c7ea63f50a13e7b9e38c65994998e2b441..2c281617f4b3bf0c41fb480b3e7a32b2e348483f 100644
--- a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod3D_impl.h
+++ b/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");
}