Refactoring

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

+/* 
+ * File:   tnlDirectEikonalMethodBase1D_impl.h
+ * Author: Fencl
+ *
+ * Created on March 15, 2019
+ */
+
+#pragma once
+
+template< typename Real,
+        typename Device,
+        typename Index >
+void
+tnlDirectEikonalMethodsBase< Meshes::Grid< 1, Real, Device, Index > >::
+initInterface( const MeshFunctionPointer& _input,
+        MeshFunctionPointer& _output,
+        InterfaceMapPointer& _interfaceMap  )
+{
+  if( std::is_same< Device, Devices::Cuda >::value )
+  {
+#ifdef HAVE_CUDA
+    const MeshType& mesh = _input->getMesh();
+    
+    const int cudaBlockSize( 16 );
+    int numBlocksX = Devices::Cuda::getNumberOfBlocks( mesh.getDimensions().x(), cudaBlockSize );
+    dim3 blockSize( cudaBlockSize );
+    dim3 gridSize( numBlocksX );
+    Devices::Cuda::synchronizeDevice();
+    CudaInitCaller<<< gridSize, blockSize >>>( _input.template getData< Device >(),
+            _output.template modifyData< Device >(),
+            _interfaceMap.template modifyData< Device >() );
+    cudaDeviceSynchronize();
+    TNL_CHECK_CUDA_DEVICE;
+#endif
+  }
+  if( std::is_same< Device, Devices::Host >::value )
+  {
+    const MeshType& mesh = _input->getMesh();
+    typedef typename MeshType::Cell Cell;
+    const MeshFunctionType& input = _input.getData();
+    MeshFunctionType& output = _output.modifyData();
+    InterfaceMapType& interfaceMap = _interfaceMap.modifyData();
+    Cell cell( mesh );
+    for( cell.getCoordinates().x() = 0;
+            cell.getCoordinates().x() < mesh.getDimensions().x();
+            cell.getCoordinates().x() ++ )
+    {
+      cell.refresh();
+      output[ cell.getIndex() ] =
+              input( cell ) >= 0 ? std::numeric_limits< RealType >::max() :
+                -std::numeric_limits< RealType >::max();
+      interfaceMap[ cell.getIndex() ] = false;
+    }
+    
+    
+    const RealType& h = mesh.getSpaceSteps().x();
+    for( cell.getCoordinates().x() = 0;
+            cell.getCoordinates().x() < mesh.getDimensions().x() - 1;
+            cell.getCoordinates().x() ++ )
+    {
+      cell.refresh();
+      const RealType& c = input( cell );      
+      if( ! cell.isBoundaryEntity()  )
+      {
+        const auto& neighbors = cell.getNeighborEntities();
+        Real pom = 0;
+        //const IndexType& c = cell.getIndex();
+        const IndexType e = neighbors.template getEntityIndex<  1 >();
+        if( c * input[ e ] <= 0 )
+        {
+          pom = TNL::sign( c )*( h * c )/( c - input[ e ]);
+          if( TNL::abs( output[ cell.getIndex() ] ) > TNL::abs( pom ) )
+            output[ cell.getIndex() ] = pom;
+          
+          pom = pom - TNL::sign( c )*h; //output[ e ] = (hx * c)/( c - input[ e ]) - hx;
+          if( TNL::abs( output[ e ] ) > TNL::abs( pom ) )
+            output[ e ] = pom; 
+          
+          interfaceMap[ cell.getIndex() ] = true;
+          interfaceMap[ e ] = true;
+        }
+      }
+    }
+  }
+}
+
+template< typename Real,
+        typename Device,
+        typename Index >
+template< typename MeshEntity >
+void
+tnlDirectEikonalMethodsBase< Meshes::Grid< 1, Real, Device, Index > >::
+updateCell( MeshFunctionType& u,
+        const MeshEntity& cell, 
+        const RealType v )
+{
+  const auto& neighborEntities = cell.template getNeighborEntities< 1 >();
+  const MeshType& mesh = cell.getMesh();
+  const RealType& h = mesh.getSpaceSteps().x();
+  const RealType value = u( cell );
+  RealType a, tmp = std::numeric_limits< RealType >::max();
+  
+  if( cell.getCoordinates().x() == 0 )
+    a = u[ neighborEntities.template getEntityIndex< 1 >() ];
+  else if( cell.getCoordinates().x() == mesh.getDimensions().x() - 1 )
+    a = u[ neighborEntities.template getEntityIndex< -1 >() ];
+  else
+  {
+    a = TNL::argAbsMin( u[ neighborEntities.template getEntityIndex< -1 >() ],
+            u[ neighborEntities.template getEntityIndex<  1 >() ] );
+  }
+  
+  if( fabs( a ) == std::numeric_limits< RealType >::max() )
+    return;
+  
+  tmp = a + TNL::sign( value ) * h/v;
+  
+  u[ cell.getIndex() ] = argAbsMin( value, tmp );
+}
+
+template< typename Real,
+        typename Device,
+        typename Index >
+__cuda_callable__
+bool
+tnlDirectEikonalMethodsBase< Meshes::Grid< 1, Real, Device, Index > >::
+updateCell( volatile Real sArray[18], int thri, const Real h, const Real v )
+{
+  const RealType value = sArray[ thri ];
+  RealType a, tmp = std::numeric_limits< RealType >::max();
+  
+  a = TNL::argAbsMin( sArray[ thri+1 ],
+          sArray[ thri-1 ] );
+  
+  if( fabs( a ) == std::numeric_limits< RealType >::max() )
+    return false;
+  
+  tmp = a + TNL::sign( value ) * h/v;
+  
+  
+  sArray[ thri ] = argAbsMin( value, tmp );
+  
+  tmp = value - sArray[ thri ];
+  if ( fabs( tmp ) >  0.001*h )
+    return true;
+  else
+    return false;
+}
+
+#ifdef HAVE_CUDA
+template < typename Real, typename Device, typename Index >
+__global__ void CudaInitCaller( const Functions::MeshFunction< Meshes::Grid< 1, Real, Device, Index > >& input, 
+        Functions::MeshFunction< Meshes::Grid< 1, Real, Device, Index > >& output,
+        Functions::MeshFunction< Meshes::Grid< 1, Real, Device, Index >, 1, bool >& interfaceMap )
+{
+  int i = threadIdx.x + blockDim.x*blockIdx.x;
+  const Meshes::Grid< 1, Real, Device, Index >& mesh = input.template getMesh< Devices::Cuda >();
+  
+  if( i < mesh.getDimensions().x()  )
+  {
+    typedef typename Meshes::Grid< 1, Real, Device, Index >::Cell Cell;
+    Cell cell( mesh );
+    cell.getCoordinates().x() = i;
+    cell.refresh();
+    const Index cind = cell.getIndex();
+    
+    
+    output[ cind ] =
+            input( cell ) >= 0 ? std::numeric_limits< Real >::max() :
+              - std::numeric_limits< Real >::max();
+    interfaceMap[ cind ] = false; 
+    
+    const Real& h = mesh.getSpaceSteps().x();
+    cell.refresh();
+    const Real& c = input( cell );
+    if( ! cell.isBoundaryEntity()  )
+    {
+      auto neighbors = cell.getNeighborEntities();
+      Real pom = 0;
+      const Index e = neighbors.template getEntityIndex< 1 >();
+      const Index w = neighbors.template getEntityIndex< -1 >();
+      if( c * input[ e ] <= 0 )
+      {
+        pom = TNL::sign( c )*( h * c )/( c - input[ e ]);
+        if( TNL::abs( output[ cind ] ) > TNL::abs( pom ) )
+          output[ cind ] = pom;                       
+        
+        interfaceMap[ cind ] = true;
+      }
+      if( c * input[ w ] <= 0 )
+      {
+        pom = TNL::sign( c )*( h * c )/( c - input[ w ]);
+        if( TNL::abs( output[ cind ] ) > TNL::abs( pom ) ) 
+          output[ cind ] = pom;
+        
+        interfaceMap[ cind ] = true;
+      }
+    }
+  }
+  
+}
+#endif
+
+

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

@@ -7,9 +7,9 @@
 
 #pragma once 
 
-#include <TNL/Meshes/Grid.h>
-#include <TNL/Functions/MeshFunction.h>
-#include <TNL/Devices/Cuda.h>
+//#include <TNL/Meshes/Grid.h>
+//#include <TNL/Functions/MeshFunction.h>
+//#include <TNL/Devices/Cuda.h>
 
 using namespace TNL;
 
@@ -63,25 +63,32 @@ class tnlDirectEikonalMethodsBase< Meshes::Grid< 2, Real, Device, Index > >
     typedef Functions::MeshFunction< MeshType, 2, bool > InterfaceMapType;
     typedef TNL::Containers::Array< int, Device, IndexType > ArrayContainer;
     typedef Containers::StaticVector< 2, Index > StaticVector;
+    
+    using MeshPointer = Pointers::SharedPointer<  MeshType >;
     using MeshFunctionPointer = Pointers::SharedPointer< MeshFunctionType >;
     using InterfaceMapPointer = Pointers::SharedPointer< InterfaceMapType >;
     
-    
+    // CALLER FOR HOST AND CUDA
     void initInterface( const MeshFunctionPointer& input,
             MeshFunctionPointer& output,
             InterfaceMapPointer& interfaceMap,
             StaticVector vLower, StaticVector vUpper );
-    
+        
+    // FOR HOST
     template< typename MeshEntity >
     __cuda_callable__ bool updateCell( MeshFunctionType& u,
             const MeshEntity& cell,
             const RealType velocity = 1.0 );
     
+    // FOR CUDA
     template< int sizeSArray >
-    __cuda_callable__ bool updateCell( volatile Real *sArray,
-            int thri, int thrj, const Real hx, const Real hy,
-            const Real velocity = 1.0 );
-    
+    __cuda_callable__ bool updateCell( volatile RealType *sArray,
+            int thri, int thrj, const RealType hx, const RealType hy,
+            const RealType velocity = 1.0 );
+        
+// FOR OPENMP WILL BE REMOVED
+    void getNeighbours( ArrayContainer BlockIterHost, int numBlockX, int numBlockY  );
+        
     template< int sizeSArray >
     void updateBlocks( const InterfaceMapType& interfaceMap,
             MeshFunctionType& aux,
@@ -108,16 +115,27 @@ class tnlDirectEikonalMethodsBase< Meshes::Grid< 3, Real, Device, Index > >
     using MeshFunctionPointer = Pointers::SharedPointer< MeshFunctionType >;
     using InterfaceMapPointer = Pointers::SharedPointer< InterfaceMapType >;      
     
+     // CALLER FOR HOST AND CUDA
     void initInterface( const MeshFunctionPointer& input,
             MeshFunctionPointer& output,
             InterfaceMapPointer& interfaceMap,
             StaticVector vLower, StaticVector vUpper );
     
+    // FOR HOST
     template< typename MeshEntity >
     __cuda_callable__ bool updateCell( MeshFunctionType& u,
             const MeshEntity& cell,
             const RealType velocity = 1.0);
     
+    // FOR CUDA
+    template< int sizeSArray >
+    __cuda_callable__ bool updateCell( volatile Real *sArray,
+            int thri, int thrj, int thrk, const RealType hx, const RealType hy, const RealType hz,
+            const RealType velocity = 1.0 );
+    
+    // OPENMP WILL BE REMOVED
+    void getNeighbours( ArrayContainer BlockIterHost, int numBlockX, int numBlockY, int numBlockZ );
+    
     template< int sizeSArray >
     void updateBlocks( const InterfaceMapType& interfaceMap,
             const MeshFunctionType& aux,
@@ -126,16 +144,15 @@ class tnlDirectEikonalMethodsBase< Meshes::Grid< 3, Real, Device, Index > >
     
     void getNeighbours( ArrayContainer& BlockIterHost, int numBlockX, int numBlockY, int numBlockZ );
     
-    template< int sizeSArray >
-    __cuda_callable__ bool updateCell3D( volatile Real *sArray,
-            int thri, int thrj, int thrk, const Real hx, const Real hy, const Real hz,
-            const Real velocity = 1.0 );
+    __cuda_callable__ RealType getNewValue( RealType valuesAndSteps[],
+           const RealType originalValue, const RealType v );
 };
 
 template < typename T1 >
 __cuda_callable__ void sortMinims( T1 pom[] );
 
 #ifdef HAVE_CUDA
+// 1D
 template < typename Real, typename Device, typename Index >
 __global__ void CudaInitCaller( const Functions::MeshFunction< Meshes::Grid< 1, Real, Device, Index > >& input, 
         Functions::MeshFunction< Meshes::Grid< 1, Real, Device, Index > >& output,
@@ -147,21 +164,25 @@ __global__ void CudaUpdateCellCaller( tnlDirectEikonalMethodsBase< Meshes::Grid<
         Functions::MeshFunction< Meshes::Grid< 1, Real, Device, Index > >& aux,
         bool *BlockIterDevice );
 
+
+
+
+// 2D
+template < typename Real, typename Device, typename Index >
+__global__ void CudaInitCaller( const Functions::MeshFunction< Meshes::Grid< 2, Real, Device, Index > >& input, 
+        Functions::MeshFunction< Meshes::Grid< 2, Real, Device, Index > >& output,
+        Functions::MeshFunction< Meshes::Grid< 2, Real, Device, Index >, 2, bool >& interfaceMap,
+        const Containers::StaticVector< 2, Index > vecLowerOverlas,
+        const Containers::StaticVector< 2, Index > vecUpperOerlaps );
+
 template < int sizeSArray, typename Real, typename Device, typename Index >
 __global__ void CudaUpdateCellCaller( tnlDirectEikonalMethodsBase< Meshes::Grid< 2, Real, Device, Index > > ptr,
         const Functions::MeshFunction< Meshes::Grid< 2, Real, Device, Index >, 2, bool >& interfaceMap,
         const Functions::MeshFunction< Meshes::Grid< 2, Real, Device, Index > >& aux,
         Functions::MeshFunction< Meshes::Grid< 2, Real, Device, Index > >& helpFunc,
-        TNL::Containers::Array< int, Devices::Cuda, Index > BlockIterDevice,
-        Containers::StaticVector< 2, Index > vLower, Containers::StaticVector< 2, Index > vUpper, int k,int oddEvenBlock =0);
-
-template< typename Real, typename Device, typename Index >
-__global__ void DeepCopy( const Functions::MeshFunction< Meshes::Grid< 2, Real, Device, Index > >& aux,
-        Functions::MeshFunction< Meshes::Grid< 2, Real, Device, Index > >& helpFunc, int copy, int k );
-
-template< typename Real, typename Device, typename Index >
-__global__ void DeepCopy( const Functions::MeshFunction< Meshes::Grid< 3, Real, Device, Index > >& aux,
-        Functions::MeshFunction< Meshes::Grid< 3, Real, Device, Index > >& helpFunc, int copy, int k );
+        TNL::Containers::Array< int, Devices::Cuda, Index > blockCalculationIndicator,
+        const Containers::StaticVector< 2, Index > vecLowerOverlaps, 
+        const Containers::StaticVector< 2, Index > vecUpperOverlaps, int oddEvenBlock =0);
 
 template < typename Index >
 __global__ void CudaParallelReduc( TNL::Containers::ArrayView< int, Devices::Cuda, Index > BlockIterDevice,
@@ -171,17 +192,13 @@ template < typename Index >
 __global__ void GetNeighbours( TNL::Containers::ArrayView< int, Devices::Cuda, Index > BlockIterDevice,
         TNL::Containers::ArrayView< int, Devices::Cuda, Index > BlockIterPom, int numBlockX, int numBlockY );
 
-template < typename Real, typename Device, typename Index >
-__global__ void CudaInitCaller( const Functions::MeshFunction< Meshes::Grid< 2, Real, Device, Index > >& input, 
-        Functions::MeshFunction< Meshes::Grid< 2, Real, Device, Index > >& output,
-        Functions::MeshFunction< Meshes::Grid< 2, Real, Device, Index >, 2, bool >& interfaceMap,
-        Containers::StaticVector< 2, Index > vLower, Containers::StaticVector< 2, Index > vUpper );
 
+// 3D
 template < typename Real, typename Device, typename Index >
 __global__ void CudaInitCaller3d( const Functions::MeshFunction< Meshes::Grid< 3, Real, Device, Index > >& input, 
         Functions::MeshFunction< Meshes::Grid< 3, Real, Device, Index > >& output,
         Functions::MeshFunction< Meshes::Grid< 3, Real, Device, Index >, 3, bool >& interfaceMap,
-        Containers::StaticVector< 3, Index > vLower, Containers::StaticVector< 3, Index > vUpper );
+        Containers::StaticVector< 3, Index > vecLowerOverlaps, Containers::StaticVector< 3, Index > vecUpperOverlaps );
 
 template < int sizeSArray, typename Real, typename Device, typename Index >
 __global__ void CudaUpdateCellCaller( tnlDirectEikonalMethodsBase< Meshes::Grid< 3, Real, Device, Index > > ptr,
@@ -196,4 +213,6 @@ __global__ void GetNeighbours3D( TNL::Containers::ArrayView< int, Devices::Cuda,
         int numBlockX, int numBlockY, int numBlockZ );
 #endif
 
-#include "tnlDirectEikonalMethodsBase_impl.h"
+#include "tnlDirectEikonalMethodBase1D_impl.h"
+#include "tnlDirectEikonalMethodBase2D_impl.h"
+#include "tnlDirectEikonalMethodBase3D_impl.h"

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

@@ -10,11 +10,10 @@
 
 #pragma once
 
-#include <TNL/Meshes/Grid.h>
-#include <TNL/Functions/Analytic/Constant.h>
-#include <TNL/Pointers/SharedPointer.h>
+//#include <TNL/Meshes/Grid.h>
+//#include <TNL/Functions/Analytic/Constant.h>
+//#include <TNL/Pointers/SharedPointer.h>
 #include "tnlDirectEikonalMethodsBase.h"
-#define ForDebug false // false <=> off
 
 
 template< typename Mesh,
@@ -88,6 +87,7 @@ class FastSweepingMethod< Meshes::Grid< 2, Real, Device, Index >, Communicator,
     typedef Anisotropy AnisotropyType;
     typedef tnlDirectEikonalMethodsBase< Meshes::Grid< 2, Real, Device, Index > > BaseType;
     typedef Communicator CommunicatorType;
+    typedef Containers::StaticVector< 2, Index > StaticVector;
     
     using MeshPointer = Pointers::SharedPointer<  MeshType >;
     using AnisotropyPointer = Pointers::SharedPointer< AnisotropyType, DeviceType >;
@@ -113,6 +113,15 @@ class FastSweepingMethod< Meshes::Grid< 2, Real, Device, Index >, Communicator,
     protected:
       
       const IndexType maxIterations;
+    
+      void setOverlaps( StaticVector& vecLowerOverlaps, StaticVector& vecUpperOverlaps,
+              const MeshPointer& mesh);
+      
+      bool goThroughSweep( const StaticVector boundsFrom, const StaticVector boundsTo, 
+              MeshFunctionType& aux, const InterfaceMapType& interfaceMap,
+              const AnisotropyPointer& anisotropy );
+      
+      void getInfoFromNeighbours( int& calculated, int& calculateAgain, const MeshPointer& mesh );
 };
 
 template< typename Real,
@@ -134,6 +143,7 @@ class FastSweepingMethod< Meshes::Grid< 3, Real, Device, Index >, Communicator,
     typedef Anisotropy AnisotropyType;
     typedef tnlDirectEikonalMethodsBase< Meshes::Grid< 3, Real, Device, Index > > BaseType;
     typedef Communicator CommunicatorType;
+    typedef Containers::StaticVector< 3, Index > StaticVector;
     
     using MeshPointer = Pointers::SharedPointer<  MeshType >;
     using AnisotropyPointer = Pointers::SharedPointer< AnisotropyType, DeviceType >;
@@ -161,6 +171,15 @@ class FastSweepingMethod< Meshes::Grid< 3, Real, Device, Index >, Communicator,
     protected:
       
       const IndexType maxIterations;
+      
+      void setOverlaps( StaticVector& vecLowerOverlaps, StaticVector& vecUpperOverlaps,
+              const MeshPointer& mesh);
+      
+      bool goThroughSweep( const StaticVector boundsFrom, const StaticVector boundsTo, 
+              MeshFunctionType& aux, const InterfaceMapType& interfaceMap,
+              const AnisotropyPointer& anisotropy );
+      
+      void getInfoFromNeighbours( int& calculated, int& calculateAgain, const MeshPointer& mesh );
 };