diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index b509b17e52aa8679097f9efef417f43b8a4b8638..983bc5784b130835a146fff6d3fa278f226bbbb1 100755
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -2,5 +2,5 @@ add_subdirectory( heat-equation )
add_subdirectory( navier-stokes )
#add_subdirectory( mean-curvature-flow )
#add_subdirectory( hamilton-jacobi )
-#add_subdirectory( hamilton-jacobi-parallel )
+add_subdirectory( hamilton-jacobi-parallel )
add_subdirectory( fast-sweeping )
diff --git a/examples/CMakeLists.txt~ b/examples/CMakeLists.txt~
index 0a74d246d3e11ec4cbf7a73840e91561dd99d11b..b509b17e52aa8679097f9efef417f43b8a4b8638 100755
--- a/examples/CMakeLists.txt~
+++ b/examples/CMakeLists.txt~
@@ -3,4 +3,4 @@ add_subdirectory( navier-stokes )
#add_subdirectory( mean-curvature-flow )
#add_subdirectory( hamilton-jacobi )
#add_subdirectory( hamilton-jacobi-parallel )
-#add_subdirectory( fast-sweeping )
+add_subdirectory( fast-sweeping )
diff --git a/examples/hamilton-jacobi-parallel/MainBuildConfig.h b/examples/hamilton-jacobi-parallel/MainBuildConfig.h
index 8ece05b9dc65909996a5f6f974d995c7a01cf82b..1b904c0c8b1d096a72a6ee8c3cc3cd1979d164b4 100644
--- a/examples/hamilton-jacobi-parallel/MainBuildConfig.h
+++ b/examples/hamilton-jacobi-parallel/MainBuildConfig.h
@@ -18,7 +18,7 @@
#ifndef MAINBUILDCONFIG_H_
#define MAINBUILDCONFIG_H_
-#include <solvers/tnlConfigTags.h>
+#include <solvers/tnlBuildConfigTags.h>
class MainBuildConfig
{
diff --git a/examples/hamilton-jacobi-parallel/main.h b/examples/hamilton-jacobi-parallel/main.h
index 692a2bc504ab735f69e293132d839d15e07b9743..ffda6cd5b4f5757ccf4c1430607b6a7bbf09c95a 100644
--- a/examples/hamilton-jacobi-parallel/main.h
+++ b/examples/hamilton-jacobi-parallel/main.h
@@ -17,7 +17,7 @@
#include "tnlParallelEikonalSolver.h"
#include "parallelEikonalConfig.h"
#include "MainBuildConfig.h"
-#include <solvers/tnlConfigTags.h>
+#include <solvers/tnlBuildConfigTags.h>
#include <operators/godunov-eikonal/parallelGodunovEikonal.h>
#include <mesh/tnlGrid.h>
#include <core/tnlDevice.h>
diff --git a/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver.h b/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver.h
index 2d049321c003c5464c2f47f0521e7c31ad2ae40e..ff5fa5be882ec5a542e8942c37e70e10efe022f9 100644
--- a/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver.h
+++ b/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver.h
@@ -20,8 +20,10 @@
#include <config/tnlParameterContainer.h>
#include <core/vectors/tnlVector.h>
#include <core/vectors/tnlStaticVector.h>
+#include <functions/tnlMeshFunction.h>
#include <core/tnlHost.h>
#include <mesh/tnlGrid.h>
+#include <mesh/grids/tnlGridEntity.h>
#include <limits.h>
#include <core/tnlDevice.h>
@@ -91,9 +93,13 @@ public:
VectorType runSubgrid( int boundaryCondition, VectorType u, int subGridID);
- VectorType u0, work_u;
+ tnlMeshFunction<MeshType> u0;
+ VectorType work_u;
IntVectorType subgridValues, boundaryConditions, unusedCell, calculationsCount;
MeshType mesh, subMesh;
+
+// tnlGridEntity< MeshType, 2, tnlGridEntityNoStencilStorage > Entity;
+
SchemeHost schemeHost;
SchemeDevice schemeDevice;
double delta, tau0, stopTime,cflCondition;
@@ -214,7 +220,8 @@ public:
VectorType runSubgrid( int boundaryCondition, VectorType u, int subGridID);
- VectorType u0, work_u;
+ tnlMeshFunction<MeshType> u0;
+ VectorType work_u;
IntVectorType subgridValues, boundaryConditions, unusedCell, calculationsCount;
MeshType mesh, subMesh;
SchemeHost schemeHost;
@@ -326,6 +333,32 @@ template <typename SchemeHost, typename SchemeDevice, typename Device>
__global__ void synchronize2CUDA3D(tnlParallelEikonalSolver<3, SchemeHost, SchemeDevice, Device, double, int >* cudaSolver);
#endif
-#include "tnlParallelEikonalSolver2D_impl.h"
+__device__
+double fabsMin( double x, double y)
+{
+ double fx = abs(x);
+
+ if(Min(fx,abs(y)) == fx)
+ return x;
+ else
+ return y;
+}
+
+__device__
+double atomicFabsMin(double* address, double val)
+{
+ unsigned long long int* address_as_ull =
+ (unsigned long long int*)address;
+ unsigned long long int old = *address_as_ull, assumed;
+ do {
+ assumed = old;
+ old = atomicCAS(address_as_ull, assumed,__double_as_longlong( fabsMin(__longlong_as_double(assumed),val) ));
+ } while (assumed != old);
+ return __longlong_as_double(old);
+}
+
+
+#include "tnlParallelEikonalSolver2D_impl.h"
+#include "tnlParallelEikonalSolver3D_impl.h"
#endif /* TNLPARALLELEIKONALSOLVER_H_ */
diff --git a/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver2D_impl.h b/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver2D_impl.h
index a989361a89da748c5d74c439252bbedd9b7d0ccf..aaa73cfb337d89c51858303d28fa887ca01c8dab 100644
--- a/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver2D_impl.h
+++ b/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver2D_impl.h
@@ -21,6 +21,7 @@
#include "tnlParallelEikonalSolver.h"
#include <core/mfilename.h>
+
template< typename SchemeHost, typename SchemeDevice, typename Device>
tnlParallelEikonalSolver<2,SchemeHost, SchemeDevice, Device, double, int>::tnlParallelEikonalSolver()
{
@@ -61,7 +62,7 @@ bool tnlParallelEikonalSolver<2,SchemeHost, SchemeDevice, Device, double, int>::
this->subMesh.setDimensions( this->n, this->n );
this->subMesh.setDomain( tnlStaticVector<2,double>(0.0, 0.0),
- tnlStaticVector<2,double>(this->mesh.getHx()*(double)(this->n), this->mesh.getHy()*(double)(this->n)) );
+ tnlStaticVector<2,double>(mesh.template getSpaceStepsProducts< 1, 0 >()*(double)(this->n), mesh.template getSpaceStepsProducts< 0, 1 >()*(double)(this->n)) );
this->subMesh.save("submesh.tnl");
@@ -71,7 +72,7 @@ bool tnlParallelEikonalSolver<2,SchemeHost, SchemeDevice, Device, double, int>::
//cout << this->mesh.getCellCenter(0) << endl;
this->delta = parameters.getParameter <double>("delta");
- this->delta *= this->mesh.getHx()*this->mesh.getHy();
+ this->delta *= mesh.template getSpaceStepsProducts< 1, 0 >()*mesh.template getSpaceStepsProducts< 0, 1 >();
cout << "Setting delta to " << this->delta << endl;
@@ -80,14 +81,14 @@ bool tnlParallelEikonalSolver<2,SchemeHost, SchemeDevice, Device, double, int>::
this->stopTime = parameters.getParameter <double>("stop-time");
this->cflCondition = parameters.getParameter <double>("cfl-condition");
- this -> cflCondition *= sqrt(this->mesh.getHx()*this->mesh.getHy());
+ this -> cflCondition *= sqrt(mesh.template getSpaceStepsProducts< 1, 0 >()*mesh.template getSpaceStepsProducts< 0, 1 >());
cout << "Setting CFL to " << this->cflCondition << endl;
stretchGrid();
this->stopTime /= (double)(this->gridCols);
this->stopTime *= (1.0+1.0/((double)(this->n) - 2.0));
cout << "Setting stopping time to " << this->stopTime << endl;
- //this->stopTime = 1.5*((double)(this->n))*parameters.getParameter <double>("stop-time")*this->mesh.getHx();
+ //this->stopTime = 1.5*((double)(this->n))*parameters.getParameter <double>("stop-time")*this->mesh.template getSpaceStepsProducts< 1, 0 >();
//cout << "Setting stopping time to " << this->stopTime << endl;
cout << "Initializating scheme..." << endl;
@@ -693,71 +694,71 @@ tnlParallelEikonalSolver<2,SchemeHost, SchemeDevice, Device, double, int>::runSu
if(x == 0 && (boundaryCondition & 4) && y ==0)
{
- if((u[subMesh.getCellYSuccessor( i )] - u[i])/subMesh.getHy() > 1.0)
+ if((u[subMesh.getCellYSuccessor( i )] - u[i])/subMesh.template getSpaceStepsProducts< 0, 1 >() > 1.0)
{
//cout << "x = 0; y = 0" << endl;
- u[i] = u[subMesh.getCellYSuccessor( i )] - subMesh.getHy();
+ u[i] = u[subMesh.getCellYSuccessor( i )] - subMesh.template getSpaceStepsProducts< 0, 1 >();
}
}
else if(x == 0 && (boundaryCondition & 4) && y == subMesh.getDimensions().y() - 1)
{
- if((u[subMesh.getCellYPredecessor( i )] - u[i])/subMesh.getHy() > 1.0)
+ if((u[subMesh.getCellYPredecessor( i )] - u[i])/subMesh.template getSpaceStepsProducts< 0, 1 >() > 1.0)
{
//cout << "x = 0; y = n" << endl;
- u[i] = u[subMesh.getCellYPredecessor( i )] - subMesh.getHy();
+ u[i] = u[subMesh.getCellYPredecessor( i )] - subMesh.template getSpaceStepsProducts< 0, 1 >();
}
}
else if(x == subMesh.getDimensions().x() - 1 && (boundaryCondition & 2) && y ==0)
{
- if((u[subMesh.getCellYSuccessor( i )] - u[i])/subMesh.getHy() > 1.0)
+ if((u[subMesh.getCellYSuccessor( i )] - u[i])/subMesh.template getSpaceStepsProducts< 0, 1 >() > 1.0)
{
//cout << "x = n; y = 0" << endl;
- u[i] = u[subMesh.getCellYSuccessor( i )] - subMesh.getHy();
+ u[i] = u[subMesh.getCellYSuccessor( i )] - subMesh.template getSpaceStepsProducts< 0, 1 >();
}
}
else if(x == subMesh.getDimensions().x() - 1 && (boundaryCondition & 2) && y == subMesh.getDimensions().y() - 1)
{
- if((u[subMesh.getCellYPredecessor( i )] - u[i])/subMesh.getHy() > 1.0)
+ if((u[subMesh.getCellYPredecessor( i )] - u[i])/subMesh.template getSpaceStepsProducts< 0, 1 >() > 1.0)
{
//cout << "x = n; y = n" << endl;
- u[i] = u[subMesh.getCellYPredecessor( i )] - subMesh.getHy();
+ u[i] = u[subMesh.getCellYPredecessor( i )] - subMesh.template getSpaceStepsProducts< 0, 1 >();
}
}
else if(y == 0 && (boundaryCondition & 8) && x ==0)
{
- if((u[subMesh.getCellXSuccessor( i )] - u[i])/subMesh.getHx() > 1.0)
+ if((u[subMesh.getCellXSuccessor( i )] - u[i])/subMesh.template getSpaceStepsProducts< 1, 0 >() > 1.0)
{
//cout << "y = 0; x = 0" << endl;
- u[i] = u[subMesh.getCellXSuccessor( i )] - subMesh.getHx();
+ u[i] = u[subMesh.getCellXSuccessor( i )] - subMesh.template getSpaceStepsProducts< 1, 0 >();
}
}
else if(y == 0 && (boundaryCondition & 8) && x == subMesh.getDimensions().x() - 1)
{
- if((u[subMesh.getCellXPredecessor( i )] - u[i])/subMesh.getHx() > 1.0)
+ if((u[subMesh.getCellXPredecessor( i )] - u[i])/subMesh.template getSpaceStepsProducts< 1, 0 >() > 1.0)
{
//cout << "y = 0; x = n" << endl;
- u[i] = u[subMesh.getCellXPredecessor( i )] - subMesh.getHx();
+ u[i] = u[subMesh.getCellXPredecessor( i )] - subMesh.template getSpaceStepsProducts< 1, 0 >();
}
}
else if(y == subMesh.getDimensions().y() - 1 && (boundaryCondition & 1) && x ==0)
{
- if((u[subMesh.getCellXSuccessor( i )] - u[i])/subMesh.getHx() > 1.0) {
+ if((u[subMesh.getCellXSuccessor( i )] - u[i])/subMesh.template getSpaceStepsProducts< 1, 0 >() > 1.0) {
//cout << "y = n; x = 0" << endl;
- u[i] = u[subMesh.getCellXSuccessor( i )] - subMesh.getHx();
+ u[i] = u[subMesh.getCellXSuccessor( i )] - subMesh.template getSpaceStepsProducts< 1, 0 >();
}
}
else if(y == subMesh.getDimensions().y() - 1 && (boundaryCondition & 1) && x == subMesh.getDimensions().x() - 1)
{
- if((u[subMesh.getCellXPredecessor( i )] - u[i])/subMesh.getHx() > 1.0)
+ if((u[subMesh.getCellXPredecessor( i )] - u[i])/subMesh.template getSpaceStepsProducts< 1, 0 >() > 1.0)
{
//cout << "y = n; x = n" << endl;
- u[i] = u[subMesh.getCellXPredecessor( i )] - subMesh.getHx();
+ u[i] = u[subMesh.getCellXPredecessor( i )] - subMesh.template getSpaceStepsProducts< 1, 0 >();
}
}
}*/
@@ -952,7 +953,9 @@ tnlParallelEikonalSolver<2,SchemeHost, SchemeDevice, Device, double, int>::runSu
double maxResidue( 1.0 );
//double lastResidue( 10000.0 );
- while( time < finalTime /*|| maxResidue > subMesh.getHx()*/)
+ tnlGridEntity<MeshType, 2, tnlGridEntityNoStencilStorage > Entity(subMesh);
+ tnlNeighbourGridEntityGetter<tnlGridEntity< MeshType, 2, tnlGridEntityNoStencilStorage >,2> neighbourEntities(Entity);
+ while( time < finalTime /*|| maxResidue > subMesh.template getSpaceStepsProducts< 1, 0 >()*/)
{
/****
* Compute the RHS
@@ -960,7 +963,10 @@ tnlParallelEikonalSolver<2,SchemeHost, SchemeDevice, Device, double, int>::runSu
for( int i = 0; i < fu.getSize(); i ++ )
{
- fu[ i ] = schemeHost.getValue( this->subMesh, i, this->subMesh.getCellCoordinates(i), u, time, boundaryCondition );
+ Entity.setCoordinates(tnlStaticVector<2,int>(i % subMesh.getDimensions().x(),i / subMesh.getDimensions().x()));
+ Entity.refresh();
+ neighbourEntities.refresh(subMesh,Entity.getIndex());
+ fu[ i ] = schemeHost.getValue( this->subMesh, i, tnlStaticVector<2,int>(i % subMesh.getDimensions().x(),i / subMesh.getDimensions().x()), u, time, boundaryCondition,neighbourEntities);
}
maxResidue = fu. absMax();
@@ -970,10 +976,10 @@ tnlParallelEikonalSolver<2,SchemeHost, SchemeDevice, Device, double, int>::runSu
/* if (maxResidue < 0.05)
cout << "Max < 0.05" << endl;*/
- if(currentTau > 1.0 * this->subMesh.getHx())
+ if(currentTau > 1.0 * this->subMesh.template getSpaceStepsProducts< 1, 0 >())
{
- //cout << currentTau << " >= " << 2.0 * this->subMesh.getHx() << endl;
- currentTau = 1.0 * this->subMesh.getHx();
+ //cout << currentTau << " >= " << 2.0 * this->subMesh.template getSpaceStepsProducts< 1, 0 >() << endl;
+ currentTau = 1.0 * this->subMesh.template getSpaceStepsProducts< 1, 0 >();
}
/*if(maxResidue > lastResidue)
currentTau *=(1.0/10.0);*/
@@ -1138,13 +1144,13 @@ void tnlParallelEikonalSolver<2,SchemeHost, SchemeDevice, Device, double, int>::
if(computeFU)
{
if(boundaryCondition == 4)
- u[l] = u[threadIdx.y * blockDim.x] + Sign(u[0])*this->subMesh.getHx()*(threadIdx.x) ;//+ 2*Sign(u[0])*this->subMesh.getHx()*(threadIdx.x+this->n);
+ u[l] = u[threadIdx.y * blockDim.x] + Sign(u[0])*this->subMesh.template getSpaceStepsProducts< 1, 0 >()*(threadIdx.x) ;//+ 2*Sign(u[0])*this->subMesh.template getSpaceStepsProducts< 1, 0 >()*(threadIdx.x+this->n);
else if(boundaryCondition == 2)
- u[l] = u[threadIdx.y * blockDim.x + blockDim.x - 1] + Sign(u[0])*this->subMesh.getHx()*(this->n - 1 - threadIdx.x);//+ 2*Sign(u[0])*this->subMesh.getHx()*(blockDim.x - threadIdx.x - 1+this->n);
+ u[l] = u[threadIdx.y * blockDim.x + blockDim.x - 1] + Sign(u[0])*this->subMesh.template getSpaceStepsProducts< 1, 0 >()*(this->n - 1 - threadIdx.x);//+ 2*Sign(u[0])*this->subMesh.template getSpaceStepsProducts< 1, 0 >()*(blockDim.x - threadIdx.x - 1+this->n);
else if(boundaryCondition == 8)
- u[l] = u[threadIdx.x] + Sign(u[0])*this->subMesh.getHx()*(threadIdx.y) ;//+ 2*Sign(u[0])*this->subMesh.getHx()*(threadIdx.y+this->n);
+ u[l] = u[threadIdx.x] + Sign(u[0])*this->subMesh.template getSpaceStepsProducts< 1, 0 >()*(threadIdx.y) ;//+ 2*Sign(u[0])*this->subMesh.template getSpaceStepsProducts< 1, 0 >()*(threadIdx.y+this->n);
else if(boundaryCondition == 1)
- u[l] = u[(blockDim.y - 1)* blockDim.x + threadIdx.x] + Sign(u[0])*this->subMesh.getHx()*(this->n - 1 - threadIdx.y) ;//+ Sign(u[0])*this->subMesh.getHx()*(blockDim.y - threadIdx.y - 1 +this->n);
+ u[l] = u[(blockDim.y - 1)* blockDim.x + threadIdx.x] + Sign(u[0])*this->subMesh.template getSpaceStepsProducts< 1, 0 >()*(this->n - 1 - threadIdx.y) ;//+ Sign(u[0])*this->subMesh.template getSpaceStepsProducts< 1, 0 >()*(blockDim.y - threadIdx.y - 1 +this->n);
}
}
@@ -1154,28 +1160,43 @@ void tnlParallelEikonalSolver<2,SchemeHost, SchemeDevice, Device, double, int>::
double fu = 0.0;
// if(threadIdx.x * threadIdx.y == 0)
// {
-// currentTau = this->tau0;
+// currentTau = finalTime;
// }
double finalTime = this->stopTime;
__syncthreads();
- if( time + currentTau > finalTime ) currentTau = finalTime - time;
+// if( time + currentTau > finalTime ) currentTau = finalTime - time;
+
+
+ tnlGridEntity<MeshType, 2, tnlGridEntityNoStencilStorage > Entity(subMesh);
+ tnlNeighbourGridEntityGetter<tnlGridEntity< MeshType, 2, tnlGridEntityNoStencilStorage >,2> neighbourEntities(Entity);
+ Entity.setCoordinates(tnlStaticVector<2,int>(i,j));
+ Entity.refresh();
+ neighbourEntities.refresh(subMesh,Entity.getIndex());
while( time < finalTime )
{
if(computeFU)
- fu = schemeHost.getValueDev( this->subMesh, l, tnlStaticVector<2,int>(i,j)/*this->subMesh.getCellCoordinates(l)*/, u, time, boundaryCondition);
+ fu = schemeHost.getValueDev( this->subMesh, l, tnlStaticVector<2,int>(i,j)/*this->subMesh.getCellCoordinates(l)*/, u, time, boundaryCondition, neighbourEntities);
- sharedTau[l]=cfl/fabs(fu);
+ sharedTau[l]=abs(cfl/fu);
if(l == 0)
{
- if(sharedTau[0] > 1.0 * this->subMesh.getHx()) sharedTau[0] = 1.0 * this->subMesh.getHx();
+ if(sharedTau[0] > 1.0 * this->subMesh.template getSpaceStepsProducts< 1, 0 >()) sharedTau[0] = 1.0 * this->subMesh.template getSpaceStepsProducts< 1, 0 >();
}
else if(l == blockDim.x*blockDim.y - 1)
if( time + sharedTau[l] > finalTime ) sharedTau[l] = finalTime - time;
+// if( (Sign(u[l]+sharedTau[l]*fu) != Sign(u[l])) && fu != 0.0 && fu != -0.0)
+// {
+// printf("orig: %10f", sharedTau[l]);
+// sharedTau[l]=abs(u[l]/(1.1*fu)) ;
+// printf(" new: %10f\n", sharedTau[l]);
+// }
+
+
if((blockDim.x == 16) && (l < 128)) sharedTau[l] = Min(sharedTau[l],sharedTau[l+128]);
__syncthreads();
@@ -1304,7 +1325,7 @@ void /*tnlParallelEikonalSolver<2,SchemeHost, SchemeDevice, Device, double, int>
// boundary_index = 0;
// }
- __threadfence();
+// __threadfence();
if((subgridValue == INT_MAX || fabs(u_cmp) + cudaSolver->delta < fabs(u) ) && (subgridValue_cmp != INT_MAX && subgridValue_cmp != -INT_MAX))
{
cudaSolver->unusedCell_cuda[gid] = 0;
@@ -1313,7 +1334,7 @@ void /*tnlParallelEikonalSolver<2,SchemeHost, SchemeDevice, Device, double, int>
cudaSolver->work_u_cuda[gid] = u_cmp;
u=u_cmp;
}
- //__threadfence();
+ __threadfence();
if(threadIdx.y == 0 && (blockIdx.y != 0)/* && (cudaSolver->currentStep & 1)*/)
{
u_cmp = cudaSolver->work_u_cuda[gid - blockDim.x*gridDim.x];
@@ -1327,7 +1348,7 @@ void /*tnlParallelEikonalSolver<2,SchemeHost, SchemeDevice, Device, double, int>
boundary_index = 0;
}
- __threadfence();
+// __threadfence();
if((subgridValue == INT_MAX || fabs(u_cmp) + cudaSolver->delta < fabs(u) ) && (subgridValue_cmp != INT_MAX && subgridValue_cmp != -INT_MAX))
{
cudaSolver->unusedCell_cuda[gid] = 0;
@@ -1336,6 +1357,7 @@ void /*tnlParallelEikonalSolver<2,SchemeHost, SchemeDevice, Device, double, int>
cudaSolver->work_u_cuda[gid] = u_cmp;
}
}
+ __threadfence();
__syncthreads();
if(threadIdx.x+threadIdx.y == 0)
@@ -1347,6 +1369,19 @@ void /*tnlParallelEikonalSolver<2,SchemeHost, SchemeDevice, Device, double, int>
2 * boundary[1] +
4 * boundary[2] +
8 * boundary[3]);
+
+
+ if(blockIdx.x+blockIdx.y ==0)
+ {
+ cudaSolver->currentStep = cudaSolver->currentStep + 1;
+ *(cudaSolver->runcuda) = 0;
+ }
+//
+// int stepValue = cudaSolver->currentStep + 4;
+// if( cudaSolver->getSubgridValueCUDA2D(blockIdx.y*gridDim.x + blockIdx.x) == -INT_MAX )
+// cudaSolver->setSubgridValueCUDA2D(blockIdx.y*gridDim.x + blockIdx.x, stepValue);
+//
+// atomicMax((cudaSolver->runcuda),cudaSolver->getBoundaryConditionCUDA2D(blockIdx.y*gridDim.x + blockIdx.x));
}
@@ -1467,11 +1502,11 @@ template <typename SchemeHost, typename SchemeDevice, typename Device>
__global__
void synchronize2CUDA2D(tnlParallelEikonalSolver<2,SchemeHost, SchemeDevice, Device, double, int >* cudaSolver)
{
- if(blockIdx.x+blockIdx.y ==0)
- {
- cudaSolver->currentStep = cudaSolver->currentStep + 1;
- *(cudaSolver->runcuda) = 0;
- }
+// if(blockIdx.x+blockIdx.y ==0)
+// {
+// cudaSolver->currentStep = cudaSolver->currentStep + 1;
+// *(cudaSolver->runcuda) = 0;
+// }
int stepValue = cudaSolver->currentStep + 4;
if( cudaSolver->getSubgridValueCUDA2D(blockIdx.y*gridDim.x + blockIdx.x) == -INT_MAX )
@@ -1500,7 +1535,7 @@ void /*tnlParallelEikonalSolver<2,SchemeHost, SchemeDevice, Device, double, int>
//this->subMesh_cuda.setDimensions( this->n_cuda, this->n_cuda );
//this->subMesh_cuda.setDomain( tnlStaticVector<2,double>(0.0, 0.0),
- //tnlStaticVector<2,double>(this->mesh_cuda.getHx()*(double)(this->n_cuda), this->mesh_cuda.getHy()*(double)(this->n_cuda)) );
+ //tnlStaticVector<2,double>(this->mesh_cuda.template getSpaceStepsProducts< 1, 0 >()*(double)(this->n_cuda), this->mesh_cuda.template getSpaceStepsProducts< 0, 1 >()*(double)(this->n_cuda)) );
//this->subMesh_cuda.save("submesh.tnl");
@@ -1510,7 +1545,7 @@ void /*tnlParallelEikonalSolver<2,SchemeHost, SchemeDevice, Device, double, int>
//cout << this->mesh.getCellCenter(0) << endl;
//this->delta_cuda = parameters.getParameter <double>("delta");
- //this->delta_cuda *= this->mesh_cuda.getHx()*this->mesh_cuda.getHy();
+ //this->delta_cuda *= this->mesh_cuda.template getSpaceStepsProducts< 1, 0 >()*this->mesh_cuda.template getSpaceStepsProducts< 0, 1 >();
//cout << "Setting delta to " << this->delta << endl;
@@ -1519,7 +1554,7 @@ void /*tnlParallelEikonalSolver<2,SchemeHost, SchemeDevice, Device, double, int>
//this->stopTime_cuda = parameters.getParameter <double>("stop-time");
//this->cflCondition_cuda = parameters.getParameter <double>("cfl-condition");
- //this -> cflCondition_cuda *= sqrt(this->mesh_cuda.getHx()*this->mesh_cuda.getHy());
+ //this -> cflCondition_cuda *= sqrt(this->mesh_cuda.template getSpaceStepsProducts< 1, 0 >()*this->mesh_cuda.template getSpaceStepsProducts< 0, 1 >());
//cout << "Setting CFL to " << this->cflCondition << endl;
////
////
@@ -1559,7 +1594,7 @@ void /*tnlParallelEikonalSolver<2,SchemeHost, SchemeDevice, Device, double, int>
//this->stopTime_cuda /= (double)(this->gridCols_cuda);
//this->stopTime_cuda *= (1.0+1.0/((double)(this->n_cuda) - 1.0));
//cout << "Setting stopping time to " << this->stopTime << endl;
- //this->stopTime_cuda = 1.5*((double)(this->n_cuda))*parameters.getParameter <double>("stop-time")*this->mesh_cuda.getHx();
+ //this->stopTime_cuda = 1.5*((double)(this->n_cuda))*parameters.getParameter <double>("stop-time")*this->mesh_cuda.template getSpaceStepsProducts< 1, 0 >();
//cout << "Setting stopping time to " << this->stopTime << endl;
//cout << "Initializating scheme..." << endl;
diff --git a/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver3D_impl.h b/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver3D_impl.h
index dbd085e5c25ba3eb3c10f64ec1f70067175e970a..00efb85a1a758e24420bbf600b2566dc7f408786 100644
--- a/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver3D_impl.h
+++ b/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver3D_impl.h
@@ -61,7 +61,7 @@ bool tnlParallelEikonalSolver<3,SchemeHost, SchemeDevice, Device, double, int>::
this->subMesh.setDimensions( this->n, this->n );
this->subMesh.setDomain( tnlStaticVector<3,double>(0.0, 0.0),
- tnlStaticVector<3,double>(this->mesh.getHx()*(double)(this->n), this->mesh.getHy()*(double)(this->n),this->mesh.getHz()*(double)(this->n)) );
+ tnlStaticVector<3,double>(mesh.template getSpaceStepsProducts< 1, 0, 0 >()*(double)(this->n), mesh.template getSpaceStepsProducts< 0, 1, 0 >()*(double)(this->n),mesh.template getSpaceStepsProducts< 0, 0, 1 >()*(double)(this->n)) );
this->subMesh.save("submesh.tnl");
@@ -71,7 +71,7 @@ bool tnlParallelEikonalSolver<3,SchemeHost, SchemeDevice, Device, double, int>::
//cout << this->mesh.getCellCenter(0) << endl;
this->delta = parameters.getParameter <double>("delta");
- this->delta *= this->mesh.getHx()*this->mesh.getHy();
+ this->delta *= mesh.template getSpaceStepsProducts< 1, 0, 0 >()*mesh.template getSpaceStepsProducts< 0, 1, 0 >();
cout << "Setting delta to " << this->delta << endl;
@@ -80,14 +80,14 @@ bool tnlParallelEikonalSolver<3,SchemeHost, SchemeDevice, Device, double, int>::
this->stopTime = parameters.getParameter <double>("stop-time");
this->cflCondition = parameters.getParameter <double>("cfl-condition");
- this -> cflCondition *= sqrt(this->mesh.getHx()*this->mesh.getHy());
+ this -> cflCondition *= sqrt(mesh.template getSpaceStepsProducts< 1, 0, 0 >()*mesh.template getSpaceStepsProducts< 0, 1, 0 >());
cout << "Setting CFL to " << this->cflCondition << endl;
stretchGrid();
this->stopTime /= (double)(this->gridCols);
this->stopTime *= (1.0+1.0/((double)(this->n) - 2.0));
cout << "Setting stopping time to " << this->stopTime << endl;
- //this->stopTime = 1.5*((double)(this->n))*parameters.getParameter <double>("stop-time")*this->mesh.getHx();
+ //this->stopTime = 1.5*((double)(this->n))*parameters.getParameter <double>("stop-time")*mesh.template getSpaceStepsProducts< 1, 0, 0 >();
//cout << "Setting stopping time to " << this->stopTime << endl;
cout << "Initializating scheme..." << endl;
@@ -962,6 +962,8 @@ tnlParallelEikonalSolver<3,SchemeHost, SchemeDevice, Device, double, int>::runSu
double maxResidue( 1.0 );
//double lastResidue( 10000.0 );
+ tnlGridEntity<MeshType, 3, tnlGridEntityNoStencilStorage > Entity(subMesh);
+ tnlNeighbourGridEntityGetter<tnlGridEntity< MeshType, 3, tnlGridEntityNoStencilStorage >,3> neighbourEntities(Entity);
while( time < finalTime /*|| maxResidue > subMesh.getHx()*/)
{
/****
@@ -970,7 +972,13 @@ tnlParallelEikonalSolver<3,SchemeHost, SchemeDevice, Device, double, int>::runSu
for( int i = 0; i < fu.getSize(); i ++ )
{
- fu[ i ] = schemeHost.getValue( this->subMesh, i, this->subMesh.getCellCoordinates(i), u, time, boundaryCondition );
+ tnlStaticVector<3,int> coords(i % subMesh.getDimensions().x(),
+ (i % subMesh.getDimensions().x()*subMesh.getDimensions().y())/ subMesh.getDimensions().x(),
+ i / subMesh.getDimensions().x()*subMesh.getDimensions().y());
+ Entity.setCoordinates(coords);
+ Entity.refresh();
+ neighbourEntities.refresh(subMesh,Entity.getIndex());
+ fu[ i ] = schemeHost.getValue( this->subMesh, i, coords,u, time, boundaryCondition );
}
maxResidue = fu. absMax();
@@ -1181,11 +1189,17 @@ void tnlParallelEikonalSolver<3,SchemeHost, SchemeDevice, Device, double, int>::
__syncthreads();
if( time + currentTau > finalTime ) currentTau = finalTime - time;
+ tnlGridEntity<MeshType, 3, tnlGridEntityNoStencilStorage > Entity(subMesh);
+ tnlNeighbourGridEntityGetter<tnlGridEntity< MeshType, 3, tnlGridEntityNoStencilStorage >,3> neighbourEntities(Entity);
+ Entity.setCoordinates(tnlStaticVector<3,int>(i,j,k));
+ Entity.refresh();
+ neighbourEntities.refresh(subMesh,Entity.getIndex());
+
while( time < finalTime )
{
if(computeFU)
- fu = schemeHost.getValueDev( this->subMesh, l, tnlStaticVector<3,int>(i,j)/*this->subMesh.getCellCoordinates(l)*/, u, time, boundaryCondition);
+ fu = schemeHost.getValueDev( this->subMesh, l, tnlStaticVector<3,int>(i,j,k)/*this->subMesh.getCellCoordinates(l)*/, u, time, boundaryCondition);
sharedTau[l]=cfl/fabs(fu);
diff --git a/src/operators/godunov-eikonal/godunovEikonal1D_impl.h b/src/operators/godunov-eikonal/godunovEikonal1D_impl.h
index 10ae8021b9d132488cbd87b85e1b9cd8f486c6d3..090a1f0316115e12b89d519e64775e1c5ce7056a 100644
--- a/src/operators/godunov-eikonal/godunovEikonal1D_impl.h
+++ b/src/operators/godunov-eikonal/godunovEikonal1D_impl.h
@@ -84,7 +84,7 @@ bool godunovEikonalScheme< tnlGrid< 1,MeshReal, Device, MeshIndex >, Real, Index
}
- h = originalMesh.getHx();
+ h = originalMesh.template getSpaceStepsProducts< 1, 0 >();
cout << "h = " << h << endl;
epsilon = parameters. getParameter< double >( "epsilon" );
diff --git a/src/operators/godunov-eikonal/godunovEikonal2D_impl.h b/src/operators/godunov-eikonal/godunovEikonal2D_impl.h
index dc819d4c7c997330cf93bfc47fbdac3339482fed..a10fb8f9c4133fd31a8f67db9057c9761e9b1c13 100644
--- a/src/operators/godunov-eikonal/godunovEikonal2D_impl.h
+++ b/src/operators/godunov-eikonal/godunovEikonal2D_impl.h
@@ -84,8 +84,8 @@ bool godunovEikonalScheme< tnlGrid< 2,MeshReal, Device, MeshIndex >, Real, Index
- hx = originalMesh.getHx();
- hy = originalMesh.getHy();
+ hx = originalMesh.template getSpaceStepsProducts< 1, 0 >();
+ hy = originalMesh.template getSpaceStepsProducts< 0, 1 >();
epsilon = parameters. getParameter< double >( "epsilon" );
diff --git a/src/operators/godunov-eikonal/godunovEikonal3D_impl.h b/src/operators/godunov-eikonal/godunovEikonal3D_impl.h
index 765e501c0b6c791df2dc6729f8cd5435bd3a76f5..4f8a31550cf8b247d3c52c2542a75b5ef3148765 100644
--- a/src/operators/godunov-eikonal/godunovEikonal3D_impl.h
+++ b/src/operators/godunov-eikonal/godunovEikonal3D_impl.h
@@ -80,8 +80,8 @@ bool godunovEikonalScheme< tnlGrid< 3,MeshReal, Device, MeshIndex >, Real, Index
}
- hx = originalMesh.getHx();
- hy = originalMesh.getHy();
+ hx = originalMesh.template getSpaceStepsProducts< 1, 0 >();
+ hy = originalMesh.template getSpaceStepsProducts< 0, 1 >();
hz = originalMesh.getHz();
epsilon = parameters. getParameter< double >( "epsilon" );
diff --git a/src/operators/godunov-eikonal/parallelGodunovEikonal.h b/src/operators/godunov-eikonal/parallelGodunovEikonal.h
index aff70555eada007bd9f4ddfff996a2ed2c3661e3..b15caf715a9f2bd008b010cf9807e7f76bfdbc22 100644
--- a/src/operators/godunov-eikonal/parallelGodunovEikonal.h
+++ b/src/operators/godunov-eikonal/parallelGodunovEikonal.h
@@ -26,6 +26,7 @@
#include <core/mfilename.h>
#include <mesh/tnlGrid.h>
#include <core/tnlCuda.h>
+#include <mesh/grids/tnlGridEntity.h>
template< typename Mesh,
@@ -147,8 +148,8 @@ public:
const CoordinatesType& coordinates,
const Vector& u,
const RealType& time,
- const IndexType boundaryCondition ) const;
-
+ const IndexType boundaryCondition,
+ const tnlNeighbourGridEntityGetter<tnlGridEntity< MeshType, 2, tnlGridEntityNoStencilStorage >,2> neighbourEntities) const;
#ifdef HAVE_CUDA
__device__
@@ -158,7 +159,8 @@ public:
const CoordinatesType& coordinates,
const RealType* u,
const RealType& time,
- const IndexType boundaryCondition) const;
+ const IndexType boundaryCondition,
+ const tnlNeighbourGridEntityGetter<tnlGridEntity< MeshType, 2, tnlGridEntityNoStencilStorage >,2> neighbourEntities) const;
#ifdef HAVE_CUDA
__device__ __host__
#endif
diff --git a/src/operators/godunov-eikonal/parallelGodunovEikonal1D_impl.h b/src/operators/godunov-eikonal/parallelGodunovEikonal1D_impl.h
index 6ab773ec88d8968a22c590b4fcf8f4fb964e74f3..cb0670518470254a201d9d025b2fa2b271d3f76c 100644
--- a/src/operators/godunov-eikonal/parallelGodunovEikonal1D_impl.h
+++ b/src/operators/godunov-eikonal/parallelGodunovEikonal1D_impl.h
@@ -84,7 +84,7 @@ bool parallelGodunovEikonalScheme< tnlGrid< 1,MeshReal, Device, MeshIndex >, Rea
}
- h = originalMesh.getHx();
+ h = originalMesh.template getSpaceStepsProducts< 1, 0 >();
cout << "h = " << h << endl;
epsilon = parameters. getParameter< double >( "epsilon" );
diff --git a/src/operators/godunov-eikonal/parallelGodunovEikonal2D_impl.h b/src/operators/godunov-eikonal/parallelGodunovEikonal2D_impl.h
index 1b4906b185d13692c4a25099a746c55651655364..5ac0dd0d9b05ab21ade02f4526c9b65f01337f5c 100644
--- a/src/operators/godunov-eikonal/parallelGodunovEikonal2D_impl.h
+++ b/src/operators/godunov-eikonal/parallelGodunovEikonal2D_impl.h
@@ -97,9 +97,9 @@ bool parallelGodunovEikonalScheme< tnlGrid< 2,MeshReal, Device, MeshIndex >, Rea
- hx = originalMesh.getHx();
+ hx = originalMesh.template getSpaceStepsProducts< 1, 0 >();
ihx = 1.0/hx;
- hy = originalMesh.getHy();
+ hy = originalMesh.template getSpaceStepsProducts< 0, 1 >();
ihy = 1.0/hy;
this->epsilon = parameters. getParameter< double >( "epsilon" );
@@ -144,7 +144,8 @@ Real parallelGodunovEikonalScheme< tnlGrid< 2, MeshReal, Device, MeshIndex >, Re
const CoordinatesType& coordinates,
const Vector& u,
const Real& time,
- const IndexType boundaryCondition ) const
+ const IndexType boundaryCondition,
+ const tnlNeighbourGridEntityGetter<tnlGridEntity< MeshType, 2, tnlGridEntityNoStencilStorage >,2> neighbourEntities ) const
{
@@ -181,24 +182,24 @@ Real parallelGodunovEikonalScheme< tnlGrid< 2, MeshReal, Device, MeshIndex >, Re
if(coordinates.x() == mesh.getDimensions().x() - 1)
- xf += u[mesh.template getCellNextToCell<-1,0>( cellIndex )];
+ xf += u[neighbourEntities.template getEntityIndex< -1, 0 >()];
else
- xf += u[mesh.template getCellNextToCell<1,0>( cellIndex )];
+ xf += u[neighbourEntities.template getEntityIndex< 1, 0 >()];
if(coordinates.x() == 0)
- xb -= u[mesh.template getCellNextToCell<1,0>( cellIndex )];
+ xb -= u[neighbourEntities.template getEntityIndex< 1, 0 >()];
else
- xb -= u[mesh.template getCellNextToCell<-1,0>( cellIndex )];
+ xb -= u[neighbourEntities.template getEntityIndex< -1, 0 >()];
if(coordinates.y() == mesh.getDimensions().y() - 1)
- yf += u[mesh.template getCellNextToCell<0,-1>( cellIndex )];
+ yf += u[neighbourEntities.template getEntityIndex< 0, -1 >()];
else
- yf += u[mesh.template getCellNextToCell<0,1>( cellIndex )];
+ yf += u[neighbourEntities.template getEntityIndex< 0, 1 >()];
if(coordinates.y() == 0)
- yb -= u[mesh.template getCellNextToCell<0,1>( cellIndex )];
+ yb -= u[neighbourEntities.template getEntityIndex< 0, 1 >()];
else
- yb -= u[mesh.template getCellNextToCell<0,-1>( cellIndex )];
+ yb -= u[neighbourEntities.template getEntityIndex< 0, -1 >()];
//xb *= ihx;
@@ -268,9 +269,9 @@ Real parallelGodunovEikonalScheme< tnlGrid< 2, MeshReal, Device, MeshIndex >, Re
// else *//*if(boundaryCondition & 4)
// xf = 0.0;
// else /**/if(coordinates.x() == mesh.getDimensions().x() - 1)
-// xf = negativePart((u[mesh.template getCellNextToCell<-1,0>( cellIndex )] - u[cellIndex])*ihx);
+// xf = negativePart((u[neighbourEntities.template getEntityIndex< -1, 0 >()] - u[cellIndex])*ihx);
// else
-// xf = negativePart((u[mesh.template getCellNextToCell<1,0>( cellIndex )] - u[cellIndex])*ihx);
+// xf = negativePart((u[neighbourEntities.template getEntityIndex< 1, 0 >()] - u[cellIndex])*ihx);
//
// /**/ /* if(boundaryCondition & 4)
// xb = (u[cellIndex] - u[mesh.getCellXPredecessor( cellIndex )])*ihx;
@@ -279,25 +280,25 @@ Real parallelGodunovEikonalScheme< tnlGrid< 2, MeshReal, Device, MeshIndex >, Re
// else /**/if(coordinates.x() == 0)
// xb = positivePart((u[cellIndex] - u[mesh.template getCellNextToCell<+1,0>( cellIndex )])*ihx);
// else
-// xb = positivePart((u[cellIndex] - u[mesh.template getCellNextToCell<-1,0>( cellIndex )])*ihx);
+// xb = positivePart((u[cellIndex] - u[neighbourEntities.template getEntityIndex< -1, 0 >()])*ihx);
//
// /**/ /* if(boundaryCondition & 1)
// yf = (u[mesh.getCellYSuccessor( cellIndex )] - u[cellIndex])*ihy;
// else *//*if(boundaryCondition & 8)
// yf = 0.0;
// else /**/if(coordinates.y() == mesh.getDimensions().y() - 1)
-// yf = negativePart((u[mesh.template getCellNextToCell<0,-1>( cellIndex )] - u[cellIndex])*ihy);
+// yf = negativePart((u[neighbourEntities.template getEntityIndex< 0, -1 >()] - u[cellIndex])*ihy);
// else
-// yf = negativePart((u[mesh.template getCellNextToCell<0,1>( cellIndex )] - u[cellIndex])*ihy);
+// yf = negativePart((u[neighbourEntities.template getEntityIndex< 0, 1 >()] - u[cellIndex])*ihy);
//
// /**/ /* if(boundaryCondition & 8)
// yb = (u[cellIndex] - u[mesh.getCellYPredecessor( cellIndex )])*ihy;
// else *//*if(boundaryCondition & 1)
// yb = 0.0;
// else /**/if(coordinates.y() == 0)
-// yb = positivePart((u[cellIndex] - u[mesh.template getCellNextToCell<0,1>( cellIndex )])*ihy);
+// yb = positivePart((u[cellIndex] - u[neighbourEntities.template getEntityIndex< 0, 1 >()])*ihy);
// else
-// yb = positivePart((u[cellIndex] - u[mesh.template getCellNextToCell<0,-1>( cellIndex )])*ihy);
+// yb = positivePart((u[cellIndex] - u[neighbourEntities.template getEntityIndex< 0, -1 >()])*ihy);
//
// if(xb - xf > 0.0)
// xf = 0.0;
@@ -322,36 +323,36 @@ Real parallelGodunovEikonalScheme< tnlGrid< 2, MeshReal, Device, MeshIndex >, Re
// else*//* if(boundaryCondition & 4)
// xf = 0.0;
// else /**/if(coordinates.x() == mesh.getDimensions().x() - 1)
-// xf = positivePart((u[mesh.template getCellNextToCell<-1,0>( cellIndex )] - u[cellIndex])*ihx);
+// xf = positivePart((u[neighbourEntities.template getEntityIndex< -1, 0 >()] - u[cellIndex])*ihx);
// else
-// xf = positivePart((u[mesh.template getCellNextToCell<1,0>( cellIndex )] - u[cellIndex])*ihx);
+// xf = positivePart((u[neighbourEntities.template getEntityIndex< 1, 0 >()] - u[cellIndex])*ihx);
//
// /**/ /* if(boundaryCondition & 4)
// xb = (u[cellIndex] - u[mesh.getCellXPredecessor( cellIndex )])*ihx;
// else*//* if(boundaryCondition & 2)
// xb = 0.0;
// else /**/if(coordinates.x() == 0)
-// xb = negativePart((u[cellIndex] - u[mesh.template getCellNextToCell<1,0>( cellIndex )])*ihx);
+// xb = negativePart((u[cellIndex] - u[neighbourEntities.template getEntityIndex< 1, 0 >()])*ihx);
// else
-// xb = negativePart((u[cellIndex] - u[mesh.template getCellNextToCell<-1,0>( cellIndex )])*ihx);
+// xb = negativePart((u[cellIndex] - u[neighbourEntities.template getEntityIndex< -1, 0 >()])*ihx);
//
// /**/ /* if(boundaryCondition & 1)
// yf = (u[mesh.getCellYSuccessor( cellIndex )] - u[cellIndex])*ihy;
// else *//*if(boundaryCondition & 8)
// yf = 0.0;
// else /**/if(coordinates.y() == mesh.getDimensions().y() - 1)
-// yf = positivePart((u[mesh.template getCellNextToCell<0,-1>( cellIndex )] - u[cellIndex])*ihy);
+// yf = positivePart((u[neighbourEntities.template getEntityIndex< 0, -1 >()] - u[cellIndex])*ihy);
// else
-// yf = positivePart((u[mesh.template getCellNextToCell<0,1>( cellIndex )] - u[cellIndex])*ihy);
+// yf = positivePart((u[neighbourEntities.template getEntityIndex< 0, 1 >()] - u[cellIndex])*ihy);
//
// /**/ /* if(boundaryCondition & 8)
// yb = (u[cellIndex] - u[mesh.getCellYPredecessor( cellIndex )])*ihy;
// else*//* if(boundaryCondition & 1)
// yb = 0.0;
// else /**/if(coordinates.y() == 0)
-// yb = negativePart((u[cellIndex] - u[mesh.template getCellNextToCell<0,1>( cellIndex )])*ihy);
+// yb = negativePart((u[cellIndex] - u[neighbourEntities.template getEntityIndex< 0, 1 >()])*ihy);
// else
-// yb = negativePart((u[cellIndex] - u[mesh.template getCellNextToCell<0,-1>( cellIndex )])*ihy);
+// yb = negativePart((u[cellIndex] - u[neighbourEntities.template getEntityIndex< 0, -1 >()])*ihy);
//
//
// if(xb - xf > 0.0)
@@ -399,37 +400,42 @@ Real parallelGodunovEikonalScheme< tnlGrid< 2, MeshReal, Device, MeshIndex >, Re
const CoordinatesType& coordinates,
const Real* u,
const Real& time,
- const IndexType boundaryCondition) const
+ const IndexType boundaryCondition,
+ const tnlNeighbourGridEntityGetter<tnlGridEntity< MeshType, 2, tnlGridEntityNoStencilStorage >,2> neighbourEntities) const
{
- RealType signui = sign(u[cellIndex],this->epsilon);
+ RealType signui;
+ if(boundaryCondition == 0)
+ signui = sign(u[cellIndex],/*(boundaryCondition != 0) * */this->epsilon);
+ else
+ signui = Sign(u[cellIndex]);
RealType xb = u[cellIndex];
RealType xf = -u[cellIndex];
RealType yb = u[cellIndex];
RealType yf = -u[cellIndex];
- RealType a,b,c;
+ RealType a,b/*,c*/;
if(coordinates.x() == mesh.getDimensions().x() - 1)
- xf += u[mesh.template getCellNextToCell<-1,0>( cellIndex )];
+ xf += u[neighbourEntities.template getEntityIndex< -1, 0 >()];
else
- xf += u[mesh.template getCellNextToCell<1,0>( cellIndex )];
+ xf += u[neighbourEntities.template getEntityIndex< 1, 0 >()];
if(coordinates.x() == 0)
- xb -= u[mesh.template getCellNextToCell<1,0>( cellIndex )];
+ xb -= u[neighbourEntities.template getEntityIndex< 1, 0 >()];
else
- xb -= u[mesh.template getCellNextToCell<-1,0>( cellIndex )];
+ xb -= u[neighbourEntities.template getEntityIndex< -1, 0 >()];
if(coordinates.y() == mesh.getDimensions().y() - 1)
- yf += u[mesh.template getCellNextToCell<0,-1>( cellIndex )];
+ yf += u[neighbourEntities.template getEntityIndex< 0, -1 >()];
else
- yf += u[mesh.template getCellNextToCell<0,1>( cellIndex )];
+ yf += u[neighbourEntities.template getEntityIndex< 0, 1 >()];
if(coordinates.y() == 0)
- yb -= u[mesh.template getCellNextToCell<0,1>( cellIndex )];
+ yb -= u[neighbourEntities.template getEntityIndex< 0, 1 >()];
else
- yb -= u[mesh.template getCellNextToCell<0,-1>( cellIndex )];
+ yb -= u[neighbourEntities.template getEntityIndex< 0, -1 >()];
if(signui > 0.0)
@@ -466,12 +472,12 @@ Real parallelGodunovEikonalScheme< tnlGrid< 2, MeshReal, Device, MeshIndex >, Re
else
b = yf;
- c =(1.0 - sqrt(a*a+b*b)*ihx );
+// c =(1.0 - sqrt(a*a+b*b)*ihx );
- if(c > 0.0 )
- return Sign(u[cellIndex])*c;
- else
- return signui*c;
+// if(c > 0.0 )
+// return Sign(u[cellIndex])*c;
+// else
+ return signui*(1.0 - sqrt(a*a+b*b)*ihx );
}
diff --git a/src/operators/godunov-eikonal/parallelGodunovEikonal3D_impl.h b/src/operators/godunov-eikonal/parallelGodunovEikonal3D_impl.h
index 7b1346dca7dcd19b704b2cbfeb53cf97589c831a..25bde7a7772a0bd7b0124fa8d0e133c1a3efbd0f 100644
--- a/src/operators/godunov-eikonal/parallelGodunovEikonal3D_impl.h
+++ b/src/operators/godunov-eikonal/parallelGodunovEikonal3D_impl.h
@@ -80,8 +80,8 @@ bool parallelGodunovEikonalScheme< tnlGrid< 3,MeshReal, Device, MeshIndex >, Rea
}
- hx = originalMesh.getHx();
- hy = originalMesh.getHy();
+ hx = originalMesh.template getSpaceStepsProducts< 1, 0 >();
+ hy = originalMesh.template getSpaceStepsProducts< 0, 1 >();
hz = originalMesh.getHz();
ihx = 1.0/hx;
ihy = 1.0/hy;