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 ed176f08a92c7ab210ddebe16e3523fa6da56966..d8770b0e735c0b287c75bdeac2d7803efdf752a2 100644
--- a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalMethodsBase.h
+++ b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalMethodsBase.h
@@ -31,10 +31,12 @@ class tnlDirectEikonalMethodsBase< Meshes::Grid< 1, Real, Device, Index > >
typedef Index IndexType;
typedef Functions::MeshFunction< MeshType > MeshFunctionType;
typedef Functions::MeshFunction< MeshType, 1, bool > InterfaceMapType;
+ using MeshFunctionPointer = SharedPointer< MeshFunctionType >;
+ using InterfaceMapPointer = SharedPointer< InterfaceMapType >;
- void initInterface( const MeshFunctionType& input,
- MeshFunctionType& output,
- InterfaceMapType& interfaceMap );
+ void initInterface( const MeshFunctionPointer& input,
+ MeshFunctionPointer& output,
+ InterfaceMapPointer& interfaceMap );
template< typename MeshEntity >
__cuda_callable__ void updateCell( MeshFunctionType& u,
@@ -55,10 +57,12 @@ class tnlDirectEikonalMethodsBase< Meshes::Grid< 2, Real, Device, Index > >
typedef Index IndexType;
typedef Functions::MeshFunction< MeshType > MeshFunctionType;
typedef Functions::MeshFunction< MeshType, 2, bool > InterfaceMapType;
+ using MeshFunctionPointer = SharedPointer< MeshFunctionType >;
+ using InterfaceMapPointer = SharedPointer< InterfaceMapType >;
- void initInterface( const MeshFunctionType& input,
- MeshFunctionType& output,
- InterfaceMapType& interfaceMap );
+ void initInterface( const MeshFunctionPointer& input,
+ MeshFunctionPointer& output,
+ InterfaceMapPointer& interfaceMap );
template< typename MeshEntity >
__cuda_callable__ void updateCell( MeshFunctionType& u,
@@ -81,10 +85,12 @@ class tnlDirectEikonalMethodsBase< Meshes::Grid< 3, Real, Device, Index > >
typedef Index IndexType;
typedef Functions::MeshFunction< MeshType > MeshFunctionType;
typedef Functions::MeshFunction< MeshType, 3, bool > InterfaceMapType;
+ using MeshFunctionPointer = SharedPointer< MeshFunctionType >;
+ using InterfaceMapPointer = SharedPointer< InterfaceMapType >;
- void initInterface( const MeshFunctionType& input,
- MeshFunctionType& output,
- InterfaceMapType& interfaceMap );
+ void initInterface( const MeshFunctionPointer& input,
+ MeshFunctionPointer& output,
+ InterfaceMapPointer& interfaceMap );
template< typename MeshEntity >
__cuda_callable__ void updateCell( MeshFunctionType& u,
@@ -109,12 +115,12 @@ template < typename Real, typename Device, typename Index >
__global__ void CudaUpdateCellCaller( Functions::MeshFunction< Meshes::Grid< 2, Real, Device, Index >, 2, bool >& interfaceMap,
Functions::MeshFunction< Meshes::Grid< 2, Real, Device, Index > >& aux );
-/*template < typename Real, typename Device, typename Index >
+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 );*/
+ Functions::MeshFunction< Meshes::Grid< 2, Real, Device, Index >, 2, bool >& interfaceMap );
-__global__ void CudaInitCaller( const Functions::MeshFunction< Meshes::Grid< 2, double, TNL::Devices::Cuda, int > >& input );
+//__global__ void CudaInitCaller( const Functions::MeshFunction< Meshes::Grid< 2, double, TNL::Devices::Cuda, int > >& input );
#endif
#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 8444b6c37801f3a428742760d687cda7a0b9ed1f..37fd245128b9c549add0653c59cfa31c9db2d7df 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
@@ -16,12 +16,15 @@ template< typename Real,
typename Index >
void
tnlDirectEikonalMethodsBase< Meshes::Grid< 1, Real, Device, Index > >::
-initInterface( const MeshFunctionType& input,
- MeshFunctionType& output,
- InterfaceMapType& interfaceMap )
+initInterface( const MeshFunctionPointer& _input,
+ MeshFunctionPointer& _output,
+ InterfaceMapPointer& _interfaceMap )
{
- const MeshType& mesh = input.getMesh();
+ 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() = 1;
cell.getCoordinates().x() < mesh.getDimensions().x() - 1;
@@ -74,9 +77,9 @@ template< typename Real,
typename Index >
void
tnlDirectEikonalMethodsBase< Meshes::Grid< 2, Real, Device, Index > >::
-initInterface( const MeshFunctionType& input,
- MeshFunctionType& output,
- InterfaceMapType& interfaceMap )
+initInterface( const MeshFunctionPointer& _input,
+ MeshFunctionPointer& _output,
+ InterfaceMapPointer& _interfaceMap )
{
/*
doplnit přepočty pro cudu:
@@ -88,7 +91,7 @@ initInterface( const MeshFunctionType& input,
if( std::is_same< Device, Devices::Cuda >::value )
{
#ifdef HAVE_CUDA
- const MeshType& mesh = input.getMesh();
+ const MeshType& mesh = _input->getMesh();
const int cudaBlockSize( 16 );
int numBlocksX = Devices::Cuda::getNumberOfBlocks( mesh.getDimensions().x(), cudaBlockSize );
@@ -96,14 +99,19 @@ initInterface( const MeshFunctionType& input,
dim3 blockSize( cudaBlockSize, cudaBlockSize );
dim3 gridSize( numBlocksX, numBlocksY );
Devices::Cuda::synchronizeDevice();
- //CudaInitCaller< Real, Device, Index ><<< gridSize, blockSize >>>( input, output, interfaceMap );
- CudaInitCaller<<< gridSize, blockSize >>>( input );
+ CudaInitCaller<<< gridSize, blockSize >>>( _input.template getData< Device >(),
+ _output.template modifyData< Device >(),
+ _interfaceMap.template modifyData< Device >() );
+ //CudaInitCaller<<< gridSize, blockSize >>>( input );
cudaDeviceSynchronize();
TNL_CHECK_CUDA_DEVICE;
#endif
}
if( std::is_same< Device, Devices::Host >::value )
{
+ const MeshFunctionType& input = _input.getData();
+ MeshFunctionType& output = _output.modifyData();
+ InterfaceMapType& interfaceMap = _interfaceMap.modifyData();
const MeshType& mesh = input.getMesh();
typedef typename MeshType::Cell Cell;
Cell cell( mesh );
@@ -289,10 +297,13 @@ template< typename Real,
typename Index >
void
tnlDirectEikonalMethodsBase< Meshes::Grid< 3, Real, Device, Index > >::
-initInterface( const MeshFunctionType& input,
- MeshFunctionType& output,
- InterfaceMapType& interfaceMap )
+initInterface( const MeshFunctionPointer& _input,
+ MeshFunctionPointer& _output,
+ InterfaceMapPointer& _interfaceMap )
{
+ const MeshFunctionType& input = _input.getData();
+ MeshFunctionType& output = _output.modifyData();
+ InterfaceMapType& interfaceMap = _interfaceMap.modifyData();
const MeshType& mesh = input.getMesh();
typedef typename MeshType::Cell Cell;
Cell cell( mesh );
@@ -602,16 +613,16 @@ __cuda_callable__ void sortMinims( T1 pom[])
}
}
-/*template < typename Real, typename Device, typename Index >
+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 )
{
int i = threadIdx.x + blockDim.x*blockIdx.x;
int j = blockDim.y*blockIdx.y + threadIdx.y;
- const Meshes::Grid< 2, Real, Device, Index >& mesh = input.getMesh();
+ const Meshes::Grid< 2, Real, Device, Index >& mesh = input.template getMesh< Devices::Cuda >();
- //if( i < mesh.getDimensions().x() && j < mesh.getDimensions().y() )
+ if( i < mesh.getDimensions().x() && j < mesh.getDimensions().y() )
{
typedef typename Meshes::Grid< 2, Real, Device, Index >::Cell Cell;
Cell cell( mesh );
@@ -671,13 +682,13 @@ __global__ void CudaInitCaller( const Functions::MeshFunction< Meshes::Grid< 2,
}
}
}
-}*/
+}
-__global__ void CudaInitCaller( const Functions::MeshFunction< Meshes::Grid< 2, double, TNL::Devices::Cuda, int > >& input )
+/*__global__ void CudaInitCaller( const Functions::MeshFunction< Meshes::Grid< 2, double, TNL::Devices::Cuda, int > >& input )
{
int i = threadIdx.x + blockDim.x*blockIdx.x;
int j = blockDim.y*blockIdx.y + threadIdx.y;
//const Meshes::Grid< 2, double, TNL::Devices::Cuda, int >& mesh = input.getMesh();
-}
+}*/
diff --git a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalProblem.h b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalProblem.h
index ca3b8368ccf7a78650ba5fc985f89a52877a9a70..66513d76f9043dd1ee1fed24f79c949e95ddf54a 100644
--- a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalProblem.h
+++ b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalProblem.h
@@ -36,6 +36,8 @@ class tnlDirectEikonalProblem
typedef Functions::MeshFunction< Mesh > MeshFunctionType;
typedef Problems::PDEProblem< Mesh, RealType, DeviceType, IndexType > BaseType;
using AnisotropyType = Anisotropy;
+ using AnisotropyPointer = SharedPointer< AnisotropyType, DeviceType >;
+ using MeshFunctionPointer = SharedPointer< MeshFunctionType >;
using typename BaseType::MeshType;
using typename BaseType::DofVectorType;
@@ -76,11 +78,11 @@ class tnlDirectEikonalProblem
protected:
- MeshFunctionType u;
+ MeshFunctionPointer u;
- MeshFunctionType initialData;
+ MeshFunctionPointer initialData;
- AnisotropyType anisotropy;
+ AnisotropyPointer anisotropy;
};
diff --git a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalProblem_impl.h b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalProblem_impl.h
index b10c6949e6ffea16eec6e7a4129186c7d0d817b5..57ac7d6896f4db6dc91389af45d9d7c6cef43ec2 100644
--- a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalProblem_impl.h
+++ b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalProblem_impl.h
@@ -95,7 +95,7 @@ tnlDirectEikonalProblem< Mesh, Anisotropy, Real, Index >::
bindDofs( const MeshPointer& mesh,
const DofVectorPointer& dofs )
{
- this->u.bind( mesh, dofs );
+ this->u->bind( mesh, dofs );
}
template< typename Mesh,
@@ -110,8 +110,8 @@ setInitialCondition( const Config::ParameterContainer& parameters,
MeshDependentDataPointer& meshdependentData )
{
String inputFile = parameters.getParameter< String >( "input-file" );
- this->initialData.setMesh( mesh );
- if( !this->initialData.boundLoad( inputFile ) )
+ this->initialData->setMesh( mesh );
+ if( !this->initialData->boundLoad( inputFile ) )
return false;
return true;
}
diff --git a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod.h b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod.h
index 61104a143c34ca4ac0ee1dedb98c4dda46996e88..48f3301320c719d2c1be7d4b9c663b754d8b57cd 100644
--- a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod.h
+++ b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod.h
@@ -34,14 +34,18 @@ class FastSweepingMethod< Meshes::Grid< 1, Real, Device, Index >, Anisotropy >
typedef Meshes::Grid< 1, Real, Device, Index > MeshType;
typedef Real RealType;
- typedef TNL::Devices::Host DeviceType;
+ typedef Device DeviceType;
typedef Index IndexType;
typedef Anisotropy AnisotropyType;
typedef tnlDirectEikonalMethodsBase< Meshes::Grid< 1, Real, Device, Index > > BaseType;
using MeshPointer = SharedPointer< MeshType >;
+ using AnisotropyPointer = SharedPointer< AnisotropyType, DeviceType >;
using typename BaseType::InterfaceMapType;
using typename BaseType::MeshFunctionType;
+ using typename BaseType::InterfaceMapPointer;
+ using typename BaseType::MeshFunctionPointer;
+
FastSweepingMethod();
@@ -50,8 +54,8 @@ class FastSweepingMethod< Meshes::Grid< 1, Real, Device, Index >, Anisotropy >
void setMaxIterations( const IndexType& maxIterations );
void solve( const MeshPointer& mesh,
- const AnisotropyType& anisotropy,
- MeshFunctionType& u );
+ const AnisotropyPointer& anisotropy,
+ MeshFunctionPointer& u );
protected:
@@ -72,14 +76,17 @@ class FastSweepingMethod< Meshes::Grid< 2, Real, Device, Index >, Anisotropy >
typedef Meshes::Grid< 2, Real, Device, Index > MeshType;
typedef Real RealType;
- typedef TNL::Devices::Host DeviceType;
+ typedef Device DeviceType;
typedef Index IndexType;
typedef Anisotropy AnisotropyType;
typedef tnlDirectEikonalMethodsBase< Meshes::Grid< 2, Real, Device, Index > > BaseType;
using MeshPointer = SharedPointer< MeshType >;
+ using AnisotropyPointer = SharedPointer< AnisotropyType, DeviceType >;
using typename BaseType::InterfaceMapType;
using typename BaseType::MeshFunctionType;
+ using typename BaseType::InterfaceMapPointer;
+ using typename BaseType::MeshFunctionPointer;
FastSweepingMethod();
@@ -88,8 +95,8 @@ class FastSweepingMethod< Meshes::Grid< 2, Real, Device, Index >, Anisotropy >
void setMaxIterations( const IndexType& maxIterations );
void solve( const MeshPointer& mesh,
- const AnisotropyType& anisotropy,
- MeshFunctionType& u );
+ const AnisotropyPointer& anisotropy,
+ MeshFunctionPointer& u );
protected:
@@ -110,14 +117,17 @@ class FastSweepingMethod< Meshes::Grid< 3, Real, Device, Index >, Anisotropy >
typedef Meshes::Grid< 3, Real, Device, Index > MeshType;
typedef Real RealType;
- typedef TNL::Devices::Host DeviceType;
+ typedef Device DeviceType;
typedef Index IndexType;
typedef Anisotropy AnisotropyType;
typedef tnlDirectEikonalMethodsBase< Meshes::Grid< 3, Real, Device, Index > > BaseType;
using MeshPointer = SharedPointer< MeshType >;
+ using AnisotropyPointer = SharedPointer< AnisotropyType, DeviceType >;
using typename BaseType::InterfaceMapType;
using typename BaseType::MeshFunctionType;
+ using typename BaseType::InterfaceMapPointer;
+ using typename BaseType::MeshFunctionPointer;
FastSweepingMethod();
@@ -126,8 +136,8 @@ class FastSweepingMethod< Meshes::Grid< 3, Real, Device, Index >, Anisotropy >
void setMaxIterations( const IndexType& maxIterations );
void solve( const MeshPointer& mesh,
- const AnisotropyType& anisotropy,
- MeshFunctionType& u );
+ const AnisotropyPointer& anisotropy,
+ MeshFunctionPointer& u );
protected:
@@ -135,8 +145,6 @@ class FastSweepingMethod< Meshes::Grid< 3, Real, Device, Index >, Anisotropy >
const IndexType maxIterations;
};
-
-
#include "tnlFastSweepingMethod1D_impl.h"
#include "tnlFastSweepingMethod2D_impl.h"
#include "tnlFastSweepingMethod3D_impl.h"
diff --git a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod1D_impl.h b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod1D_impl.h
index 874945eafd7f0cb1eff4667eba6ef026ba026634..105f1ceda9b7af99e896e9aadae910330a543004 100644
--- a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod1D_impl.h
+++ b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod1D_impl.h
@@ -55,16 +55,16 @@ template< typename Real,
void
FastSweepingMethod< Meshes::Grid< 1, Real, Device, Index >, Anisotropy >::
solve( const MeshPointer& mesh,
- const AnisotropyType& anisotropy,
- MeshFunctionType& u )
+ const AnisotropyPointer& anisotropy,
+ MeshFunctionPointer& u )
{
- MeshFunctionType aux;
- InterfaceMapType interfaceMap;
- aux.setMesh( mesh );
- interfaceMap.setMesh( mesh );
+ MeshFunctionPointer aux;
+ InterfaceMapPointer interfaceMap;
+ aux->setMesh( 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" );
}
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 48640e3fc5772c5f629dcaf9c4b66b4a63e9b179..e36224d075140899ad26ae7d7e1e2348b7d0e254 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
@@ -55,22 +55,24 @@ template< typename Real,
void
FastSweepingMethod< Meshes::Grid< 2, Real, Device, Index >, Anisotropy >::
solve( const MeshPointer& mesh,
- const AnisotropyType& anisotropy,
- MeshFunctionType& u )
+ const AnisotropyPointer& anisotropy,
+ MeshFunctionPointer& u )
{
- MeshFunctionType aux;
- InterfaceMapType interfaceMap;
- aux.setMesh( mesh );
- interfaceMap.setMesh( mesh );
+ MeshFunctionPointer auxPtr;
+ InterfaceMapPointer interfaceMapPtr;
+ auxPtr->setMesh( mesh );
+ interfaceMapPtr->setMesh( mesh );
std::cout << "Initiating the interface cells ..." << std::endl;
- BaseType::initInterface( u, aux, interfaceMap );
+ BaseType::initInterface( u, auxPtr, interfaceMapPtr );
cudaDeviceSynchronize();
- aux.save( "aux-ini.tnl" );
+ auxPtr->save( "aux-ini.tnl" );
typename MeshType::Cell cell( *mesh );
IndexType iteration( 0 );
+ InterfaceMapType interfaceMap = *interfaceMapPtr;
+ MeshFunctionType aux = *auxPtr;
while( iteration < this->maxIterations )
{
if( std::is_same< DeviceType, Devices::Host >::value )
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 9844c4fcaf741cac258355f9d40cb0dae08bc7a7..1fd8c9f2a0f4c2ef68ec3185f1fc65983837cc56 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
@@ -55,20 +55,22 @@ template< typename Real,
void
FastSweepingMethod< Meshes::Grid< 3, Real, Device, Index >, Anisotropy >::
solve( const MeshPointer& mesh,
- const AnisotropyType& anisotropy,
- MeshFunctionType& u )
+ const AnisotropyPointer& anisotropy,
+ MeshFunctionPointer& u )
{
- MeshFunctionType aux;
- InterfaceMapType interfaceMap;
- aux.setMesh( mesh );
- interfaceMap.setMesh( mesh );
+ MeshFunctionPointer auxPtr;
+ InterfaceMapPointer interfaceMapPtr;
+ auxPtr->setMesh( mesh );
+ interfaceMapPtr->setMesh( mesh );
std::cout << "Initiating the interface cells ..." << std::endl;
- BaseType::initInterface( u, aux, interfaceMap );
- aux.save( "aux-ini.tnl" );
+ BaseType::initInterface( u, auxPtr, interfaceMapPtr );
+ auxPtr->save( "aux-ini.tnl" );
typename MeshType::Cell cell( *mesh );
IndexType iteration( 0 );
+ MeshFunctionType aux = *auxPtr;
+ InterfaceMapType interfaceMap = * interfaceMapPtr;
while( iteration < this->maxIterations )
{
for( cell.getCoordinates().z() = 0;