diff --git a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/HamiltonJacobiProblem_impl.h b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/HamiltonJacobiProblem_impl.h index b529b09fd7a492f046fb9990cdccc49c1d6704fc..3cc6388497c062865a812cca629c57d45d715d09 100644 --- a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/HamiltonJacobiProblem_impl.h +++ b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/HamiltonJacobiProblem_impl.h @@ -121,14 +121,14 @@ setInitialCondition( const Config::ParameterContainer& parameters, DofVectorType& dofs, MeshDependentDataType& meshDependentData ) { - this->bindDofs( mesh, dofs ); - const String& initialConditionFile = parameters.getParameter< String >( "initial-condition" ); - if( ! this->solution.boundLoad( initialConditionFile ) ) - { - std::cerr << "I am not able to load the initial condition from the file " << initialConditionFile << "." <<std::endl; - return false; - } - return true; + this->bindDofs( mesh, dofs ); + const String& initialConditionFile = parameters.getParameter< String >( "initial-condition" ); + if( ! this->solution.boundLoad( initialConditionFile ) ) + { + std::cerr << "I am not able to load the initial condition from the file " << initialConditionFile << "." <<std::endl; + return false; + } + return true; } template< typename Mesh, 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 f712ce2cc9ff3de5701a8550722b8904d33c4229..76d98a2de79882b1a5c2b89b6c73dc30b34e59b8 100644 --- a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalMethodsBase.h +++ b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalMethodsBase.h @@ -127,9 +127,6 @@ class tnlDirectEikonalMethodsBase< Meshes::Grid< 3, Real, Device, Index > > const Real velocity = 1.0 ); }; -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 > __cuda_callable__ void sortMinims( T1 pom[] ); 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 8f7937541d4708f519cd73a1cb50f0b582fbad48..786c6a48862d6107904107d9b75ab4802d05557f 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 @@ -102,7 +102,7 @@ updateBlocks( InterfaceMapType interfaceMap, ArrayContainer BlockIterHost, int numThreadsPerBlock/*, Real **sArray*/ ) { #pragma omp parallel for schedule( dynamic ) - for( int i = 0; i < BlockIterHost.getSize(); i++ ) + for( IndexType i = 0; i < BlockIterHost.getSize(); i++ ) { if( BlockIterHost[ i ] ) { @@ -142,14 +142,14 @@ updateBlocks( InterfaceMapType interfaceMap, if( sArray == nullptr ) std::cout << "Error while allocating memory for sArray." << std::endl; - for( int thri = 0; thri < sizeSArray; thri++ ){ - for( int thrj = 0; thrj < sizeSArray; thrj++ ) + for( IndexType thri = 0; thri < sizeSArray; thri++ ){ + for( IndexType thrj = 0; thrj < sizeSArray; thrj++ ) sArray[ thri * sizeSArray + thrj ] = std::numeric_limits< Real >::max(); } //printf("numThreadsPerBlock = %d\n", numThreadsPerBlock); - for( int thrj = 0; thrj < numThreadsPerBlock + 1; thrj++ ) + for( IndexType thrj = 0; thrj < numThreadsPerBlock + 1; thrj++ ) { if( dimX > (blIdx+1) * numThreadsPerBlock && thrj+1 < ykolik ) sArray[ ( thrj+1 )* sizeSArray +xkolik] = aux[ blIdy*numThreadsPerBlock*dimX - dimX + blIdx*numThreadsPerBlock - 1 + (thrj+1)*dimX + xkolik ]; @@ -165,14 +165,14 @@ updateBlocks( InterfaceMapType interfaceMap, sArray[thrj+1] = aux[ blIdy*numThreadsPerBlock*dimX - dimX + blIdx*numThreadsPerBlock - 1 + thrj+1 ]; } - for( int k = 0; k < numThreadsPerBlock; k++ ){ - for( int l = 0; l < numThreadsPerBlock; l++ ) + for( IndexType k = 0; k < numThreadsPerBlock; k++ ){ + for( IndexType l = 0; l < numThreadsPerBlock; l++ ) if( blIdy * numThreadsPerBlock + k < dimY && blIdx * numThreadsPerBlock + l < dimX ) sArray[(k+1) * sizeSArray + l+1] = aux[ blIdy * numThreadsPerBlock * dimX + numThreadsPerBlock * blIdx + k*dimX + l ]; } - for( int k = 0; k < numThreadsPerBlock; k++ ){ - for( int l = 0; l < numThreadsPerBlock; l++ ){ + for( IndexType k = 0; k < numThreadsPerBlock; k++ ){ + for( IndexType l = 0; l < numThreadsPerBlock; l++ ){ if( blIdy * numThreadsPerBlock + k < dimY && blIdx * numThreadsPerBlock + l < dimX ){ //std::cout << "proslo i = " << k * numThreadsPerBlock + l << std::endl; if( ! interfaceMap[ blIdy * numThreadsPerBlock * dimX + numThreadsPerBlock * blIdx + k*dimX + l ] ) @@ -191,8 +191,8 @@ updateBlocks( InterfaceMapType interfaceMap, std::cout << std::endl; }*/ - for( int k = 0; k < numThreadsPerBlock; k++ ) - for( int l = numThreadsPerBlock-1; l >-1; l-- ) { + for( IndexType k = 0; k < numThreadsPerBlock; k++ ) + for( IndexType l = numThreadsPerBlock-1; l >-1; l-- ) { if( blIdy * numThreadsPerBlock + k < dimY && blIdx * numThreadsPerBlock + l < dimX ) { if( ! interfaceMap[ blIdy * numThreadsPerBlock * dimX + numThreadsPerBlock * blIdx + k*dimX + l ] ) @@ -209,8 +209,8 @@ updateBlocks( InterfaceMapType interfaceMap, std::cout << std::endl; }*/ - for( int k = numThreadsPerBlock-1; k > -1; k-- ) - for( int l = 0; l < numThreadsPerBlock; l++ ) { + for( IndexType k = numThreadsPerBlock-1; k > -1; k-- ) + for( IndexType l = 0; l < numThreadsPerBlock; l++ ) { if( blIdy * numThreadsPerBlock + k < dimY && blIdx * numThreadsPerBlock + l < dimX ) { if( ! interfaceMap[ blIdy * numThreadsPerBlock * dimX + numThreadsPerBlock * blIdx + k*dimX + l ] ) @@ -227,8 +227,8 @@ updateBlocks( InterfaceMapType interfaceMap, std::cout << std::endl; }*/ - for( int k = numThreadsPerBlock-1; k > -1; k-- ){ - for( int l = numThreadsPerBlock-1; l >-1; l-- ) { + for( IndexType k = numThreadsPerBlock-1; k > -1; k-- ){ + for( IndexType l = numThreadsPerBlock-1; l >-1; l-- ) { if( blIdy * numThreadsPerBlock + k < dimY && blIdx * numThreadsPerBlock + l < dimX ) { if( ! interfaceMap[ blIdy * numThreadsPerBlock * dimX + numThreadsPerBlock * blIdx + k*dimX + l ] ) @@ -252,8 +252,8 @@ updateBlocks( InterfaceMapType interfaceMap, } - for( int k = 0; k < numThreadsPerBlock; k++ ){ - for( int l = 0; l < numThreadsPerBlock; l++ ) { + for( IndexType k = 0; k < numThreadsPerBlock; k++ ){ + for( IndexType l = 0; l < numThreadsPerBlock; l++ ) { if( blIdy * numThreadsPerBlock + k < dimY && blIdx * numThreadsPerBlock + l < dimX ) helpFunc[ blIdy * numThreadsPerBlock * dimX + numThreadsPerBlock * blIdx + k*dimX + l ] = sArray[ (k + 1)* sizeSArray + l + 1 ]; //std::cout<< sArray[k+1][l+1]; @@ -276,7 +276,7 @@ updateBlocks( const InterfaceMapType interfaceMap, ArrayContainer BlockIterHost, int numThreadsPerBlock/*, Real **sArray*/ ) { //#pragma omp parallel for schedule( dynamic ) - for( int i = 0; i < BlockIterHost.getSize(); i++ ) + for( IndexType i = 0; i < BlockIterHost.getSize(); i++ ) { if( BlockIterHost[ i ] ) { @@ -323,15 +323,15 @@ updateBlocks( const InterfaceMapType interfaceMap, if( sArray == nullptr ) std::cout << "Error while allocating memory for sArray." << std::endl; - for( int k = 0; k < sizeSArray; k++ ) - for( int l = 0; l < sizeSArray; l++ ) - for( int m = 0; m < sizeSArray; m++ ){ + for( IndexType k = 0; k < sizeSArray; k++ ) + for( IndexType l = 0; l < sizeSArray; l++ ) + for( IndexType m = 0; m < sizeSArray; m++ ){ sArray[ m * sizeSArray * sizeSArray + k * sizeSArray + l ] = std::numeric_limits< Real >::max(); } - for( int thrk = 0; thrk < numThreadsPerBlock; thrk++ ) - for( int thrj = 0; thrj < numThreadsPerBlock; thrj++ ) + for( IndexType thrk = 0; thrk < numThreadsPerBlock; thrk++ ) + for( IndexType thrj = 0; thrj < numThreadsPerBlock; thrj++ ) { if( blIdx != 0 && thrj+1 < ykolik && thrk+1 < zkolik ) sArray[(thrk+1 )* sizeSArray * sizeSArray + (thrj+1)*sizeSArray + 0] = @@ -358,9 +358,9 @@ updateBlocks( const InterfaceMapType interfaceMap, aux[ (blIdz+1)*numThreadsPerBlock * dimX * dimY + blIdy * numThreadsPerBlock*dimX + blIdx*numThreadsPerBlock + thrj * dimX + thrk ]; } - for( int m = 0; m < numThreadsPerBlock; m++ ){ - for( int k = 0; k < numThreadsPerBlock; k++ ){ - for( int l = 0; l < numThreadsPerBlock; l++ ){ + for( IndexType m = 0; m < numThreadsPerBlock; m++ ){ + for( IndexType k = 0; k < numThreadsPerBlock; k++ ){ + for( IndexType l = 0; l < numThreadsPerBlock; l++ ){ if( blIdy * numThreadsPerBlock + k < dimY && blIdx * numThreadsPerBlock + l < dimX && blIdz * numThreadsPerBlock + m < dimZ ) sArray[(m+1) * sizeSArray * sizeSArray + (k+1) *sizeSArray + l+1] = aux[ blIdz * numThreadsPerBlock * dimX * dimY + blIdy * numThreadsPerBlock * dimX + blIdx*numThreadsPerBlock + m*dimX*dimY + k*dimX + l ]; @@ -379,9 +379,9 @@ updateBlocks( const InterfaceMapType interfaceMap, s = "helpFunc-"+ std::to_string(numWhile) + ".tnl"; helpFunc.save( s );*/ - for( int m = 0; m < numThreadsPerBlock; m++ ){ - for( int k = 0; k < numThreadsPerBlock; k++ ){ - for( int l = 0; l < numThreadsPerBlock; l++ ){ + for( IndexType m = 0; m < numThreadsPerBlock; m++ ){ + for( IndexType k = 0; k < numThreadsPerBlock; k++ ){ + for( IndexType l = 0; l < numThreadsPerBlock; l++ ){ if( blIdy * numThreadsPerBlock + k < dimY && blIdx * numThreadsPerBlock + l < dimX && blIdz * numThreadsPerBlock + m < dimZ ){ //std::cout << "proslo i = " << k * numThreadsPerBlock + l << std::endl; if( ! interfaceMap[ blIdz * numThreadsPerBlock * dimX * dimY + blIdy * numThreadsPerBlock * dimX + blIdx*numThreadsPerBlock + m*dimX*dimY + k*dimX + l ] ) @@ -404,9 +404,9 @@ updateBlocks( const InterfaceMapType interfaceMap, s = "helpFunc-"+ std::to_string(numWhile) + ".tnl"; helpFunc.save( s );*/ - for( int m = numThreadsPerBlock-1; m >-1; m-- ){ - for( int k = 0; k < numThreadsPerBlock; k++ ){ - for( int l = 0; l <numThreadsPerBlock; l++ ){ + for( IndexType m = numThreadsPerBlock-1; m >-1; m-- ){ + for( IndexType k = 0; k < numThreadsPerBlock; k++ ){ + for( IndexType l = 0; l <numThreadsPerBlock; l++ ){ if( blIdy * numThreadsPerBlock + k < dimY && blIdx * numThreadsPerBlock + l < dimX && blIdz * numThreadsPerBlock + m < dimZ ) { if( ! interfaceMap[ blIdz * numThreadsPerBlock * dimX * dimY + blIdy * numThreadsPerBlock * dimX + blIdx*numThreadsPerBlock + m*dimX*dimY + k*dimX + l ] ) @@ -427,9 +427,9 @@ updateBlocks( const InterfaceMapType interfaceMap, s = "helpFunc-"+ std::to_string(numWhile) + ".tnl"; helpFunc.save( s );*/ - for( int m = 0; m < numThreadsPerBlock; m++ ){ - for( int k = 0; k < numThreadsPerBlock; k++ ){ - for( int l = numThreadsPerBlock-1; l >-1; l-- ){ + for( IndexType m = 0; m < numThreadsPerBlock; m++ ){ + for( IndexType k = 0; k < numThreadsPerBlock; k++ ){ + for( IndexType l = numThreadsPerBlock-1; l >-1; l-- ){ if( blIdy * numThreadsPerBlock + k < dimY && blIdx * numThreadsPerBlock + l < dimX && blIdz * numThreadsPerBlock + m < dimZ ) { if( ! interfaceMap[ blIdz * numThreadsPerBlock * dimX * dimY + blIdy * numThreadsPerBlock * dimX + blIdx*numThreadsPerBlock + m*dimX*dimY + k*dimX + l ] ) @@ -450,9 +450,9 @@ updateBlocks( const InterfaceMapType interfaceMap, s = "helpFunc-"+ std::to_string(numWhile) + ".tnl"; helpFunc.save( s ); */ - for( int m = numThreadsPerBlock-1; m >-1; m-- ){ - for( int k = 0; k < numThreadsPerBlock; k++ ){ - for( int l = numThreadsPerBlock-1; l >-1; l-- ){ + for( IndexType m = numThreadsPerBlock-1; m >-1; m-- ){ + for( IndexType k = 0; k < numThreadsPerBlock; k++ ){ + for( IndexType l = numThreadsPerBlock-1; l >-1; l-- ){ if( blIdy * numThreadsPerBlock + k < dimY && blIdx * numThreadsPerBlock + l < dimX && blIdz * numThreadsPerBlock + m < dimZ ) { if( ! interfaceMap[ blIdz * numThreadsPerBlock * dimX * dimY + blIdy * numThreadsPerBlock * dimX + blIdx*numThreadsPerBlock + m*dimX*dimY + k*dimX + l ] ) @@ -473,9 +473,9 @@ updateBlocks( const InterfaceMapType interfaceMap, s = "helpFunc-"+ std::to_string(numWhile) + ".tnl"; helpFunc.save( s );*/ - for( int m = 0; m < numThreadsPerBlock; m++ ){ - for( int k = numThreadsPerBlock-1; k > -1; k-- ){ - for( int l = 0; l <numThreadsPerBlock; l++ ){ + for( IndexType m = 0; m < numThreadsPerBlock; m++ ){ + for( IndexType k = numThreadsPerBlock-1; k > -1; k-- ){ + for( IndexType l = 0; l <numThreadsPerBlock; l++ ){ if( blIdy * numThreadsPerBlock + k < dimY && blIdx * numThreadsPerBlock + l < dimX && blIdz * numThreadsPerBlock + m < dimZ ) { if( ! interfaceMap[blIdz * numThreadsPerBlock * dimX * dimY + blIdy * numThreadsPerBlock * dimX + blIdx*numThreadsPerBlock + m*dimX*dimY + k*dimX + l ] ) @@ -496,9 +496,9 @@ updateBlocks( const InterfaceMapType interfaceMap, s = "helpFunc-"+ std::to_string(numWhile) + ".tnl"; helpFunc.save( s );*/ - for( int m = numThreadsPerBlock-1; m >-1; m-- ){ - for( int k = numThreadsPerBlock-1; k > -1; k-- ){ - for( int l = 0; l <numThreadsPerBlock; l++ ){ + for( IndexType m = numThreadsPerBlock-1; m >-1; m-- ){ + for( IndexType k = numThreadsPerBlock-1; k > -1; k-- ){ + for( IndexType l = 0; l <numThreadsPerBlock; l++ ){ if( blIdy * numThreadsPerBlock + k < dimY && blIdx * numThreadsPerBlock + l < dimX && blIdz * numThreadsPerBlock + m < dimZ ) { if( ! interfaceMap[ blIdz * numThreadsPerBlock * dimX * dimY + blIdy * numThreadsPerBlock * dimX + blIdx*numThreadsPerBlock + m*dimX*dimY + k*dimX + l ] ) @@ -519,9 +519,9 @@ updateBlocks( const InterfaceMapType interfaceMap, s = "helpFunc-"+ std::to_string(numWhile) + ".tnl"; helpFunc.save( s );*/ - for( int m = 0; m < numThreadsPerBlock; m++ ){ - for( int k = numThreadsPerBlock-1; k > -1; k-- ){ - for( int l = numThreadsPerBlock-1; l >-1; l-- ){ + for( IndexType m = 0; m < numThreadsPerBlock; m++ ){ + for( IndexType k = numThreadsPerBlock-1; k > -1; k-- ){ + for( IndexType l = numThreadsPerBlock-1; l >-1; l-- ){ if( blIdy * numThreadsPerBlock + k < dimY && blIdx * numThreadsPerBlock + l < dimX && blIdz * numThreadsPerBlock + m < dimZ ) { if( ! interfaceMap[ blIdz * numThreadsPerBlock * dimX * dimY + blIdy * numThreadsPerBlock * dimX + blIdx*numThreadsPerBlock + m*dimX*dimY + k*dimX + l ] ) @@ -543,9 +543,9 @@ updateBlocks( const InterfaceMapType interfaceMap, helpFunc.save( s );*/ - for( int m = numThreadsPerBlock-1; m >-1; m-- ){ - for( int k = numThreadsPerBlock-1; k > -1; k-- ){ - for( int l = numThreadsPerBlock-1; l >-1; l-- ){ + for( IndexType m = numThreadsPerBlock-1; m >-1; m-- ){ + for( IndexType k = numThreadsPerBlock-1; k > -1; k-- ){ + for( IndexType l = numThreadsPerBlock-1; l >-1; l-- ){ if( blIdy * numThreadsPerBlock + k < dimY && blIdx * numThreadsPerBlock + l < dimX && blIdz * numThreadsPerBlock + m < dimZ ) { if( ! interfaceMap[ blIdz * numThreadsPerBlock * dimX * dimY + blIdy * numThreadsPerBlock * dimX + blIdx*numThreadsPerBlock + m*dimX*dimY + k*dimX + l ] ) @@ -571,9 +571,9 @@ updateBlocks( const InterfaceMapType interfaceMap, } - for( int k = 0; k < numThreadsPerBlock; k++ ){ - for( int l = 0; l < numThreadsPerBlock; l++ ) { - for( int m = 0; m < numThreadsPerBlock; m++ ){ + for( IndexType k = 0; k < numThreadsPerBlock; k++ ){ + for( IndexType l = 0; l < numThreadsPerBlock; l++ ) { + for( IndexType m = 0; m < numThreadsPerBlock; m++ ){ if( blIdy * numThreadsPerBlock + k < dimY && blIdx * numThreadsPerBlock + l < dimX && blIdz * numThreadsPerBlock + m < dimZ ){ helpFunc[ blIdz * numThreadsPerBlock * dimX * dimY + blIdy * numThreadsPerBlock * dimX + blIdx*numThreadsPerBlock + m*dimX*dimY + k*dimX + l ] = sArray[ (m+1) * sizeSArray * sizeSArray + (k+1) *sizeSArray + l+1 ]; @@ -1224,26 +1224,6 @@ updateCell( MeshFunctionType& u, } } -template < typename T1, typename T2 > -T1 meet2DCondition( T1 a, T1 b, const T2 ha, const T2 hb, const T1 value, double v) -{ - T1 tmp; - if( fabs( a ) != std::numeric_limits< T1 >::max && - fabs( b ) != std::numeric_limits< T1 >::max && - fabs( a - b ) < ha/v )//TNL::sqrt( (ha * ha + hb * hb)/2 )/v ) - { - tmp = ( ha * ha * b + hb * hb * a + - TNL::sign( value ) * ha * hb * TNL::sqrt( ( ha * ha + hb * hb )/( v * v ) - - ( a - b ) * ( a - b ) ) )/( ha * ha + hb * hb ); - } - else - { - tmp = std::numeric_limits< T1 >::max; - } - - return tmp; -} - template < typename T1 > __cuda_callable__ void sortMinims( T1 pom[] ) { @@ -1274,7 +1254,7 @@ __cuda_callable__ void sortMinims( T1 pom[] ) tmp[3] = pom[5]; tmp[4] = pom[4]; tmp[5] = pom[3]; } - for( int i = 0; i < 6; i++ ) + for( unsigned int i = 0; i < 6; i++ ) { pom[ i ] = tmp[ i ]; } 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 0dde93802513e8cd3ca576b35cb8093d0d54b070..61465bee9f9784130a547dfe7aaf8bdac6ebdc2d 100644 --- a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalProblem.h +++ b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalProblem.h @@ -18,6 +18,8 @@ #include <TNL/Pointers/SharedPointer.h> #include "tnlFastSweepingMethod.h" +#include <TNL/Meshes/DistributedMeshes/DistributedGridIO.h> + template< typename Mesh, typename Communicator, typename Anisotropy, @@ -46,6 +48,8 @@ class tnlDirectEikonalProblem using MeshPointer = Pointers::SharedPointer< MeshType >; using DofVectorPointer = Pointers::SharedPointer< DofVectorType >; + typedef Communicator CommunicatorType; + static constexpr bool isTimeDependent() { return false; }; static String getType(); @@ -79,6 +83,8 @@ class tnlDirectEikonalProblem AnisotropyPointer anisotropy; -}; + Meshes::DistributedMeshes::DistrGridIOTypes distributedIOType; + + }; #include "tnlDirectEikonalProblem_impl.h" 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 3d5fa15dc23d82eb886d344bc7cf4fff7f4e34d3..7c65ef94d4818c2fe6b5b703dfe390f23ea37d21 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 @@ -113,11 +113,26 @@ tnlDirectEikonalProblem< Mesh, Communicator, Anisotropy, Real, Index >:: setInitialCondition( const Config::ParameterContainer& parameters, DofVectorPointer& dofs ) { - String inputFile = parameters.getParameter< String >( "input-file" ); - this->initialData->setMesh( this->getMesh() ); - if( !this->initialData->boundLoad( inputFile ) ) - return false; - return true; + this->bindDofs( dofs ); + String inputFile = parameters.getParameter< String >( "input-file" ); + this->initialData->setMesh( this->getMesh() ); + std::cout<<"setInitialCondition" <<std::endl; + if(CommunicatorType::isDistributed()) + { + std::cout<<"Nodes Distribution: " << u->getMesh().getDistributedMesh()->printProcessDistr() << std::endl; + if(distributedIOType==Meshes::DistributedMeshes::MpiIO) + Meshes::DistributedMeshes::DistributedGridIO<MeshFunctionType,Meshes::DistributedMeshes::MpiIO> ::load(inputFile, *u ); + if(distributedIOType==Meshes::DistributedMeshes::LocalCopy) + Meshes::DistributedMeshes::DistributedGridIO<MeshFunctionType,Meshes::DistributedMeshes::LocalCopy> ::load(inputFile, *u ); + u->template synchronize<CommunicatorType>(); + } + else + { + if( !this->initialData->boundLoad( inputFile ) ) + std::cerr << "I am not able to load the initial condition from the file " << inputFile << "." << std::endl; + return false; + } + return true; } 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 d5ce1efe164bc835866dc7fc95909b3defdf846f..001f4c69a45e98b245b8e470f489db65c84956bf 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 @@ -83,13 +83,34 @@ solve( const MeshPointer& mesh, MeshFunctionType aux = *auxPtr; -//#ifdef HAVE_MPI - bool a = Communicators::MpiCommunicator::IsInitialized(); +#ifdef HAVE_MPI + int i = Communicators::MpiCommunicator::GetRank( Communicators::MpiCommunicator::AllGroup ); + //printf( "Hello world from rank: %d ", i ); + //Communicators::MpiCommunicator::Request r = Communicators::MpiCommunicator::ISend( auxPtr, 0, 0, Communicators::MpiCommunicator::AllGroup ); + if( i == 1 ) + /*for( int k = 0; k < 16*16; k++ ) + aux[ k ] = 10;*/ + printf( "1: mesh x: %d\n", mesh->getDimensions().x() ); + printf( "1: mesh y: %d\n", mesh->getDimensions().y() ); + //aux.save("aux_proc1.tnl"); + if( i == 0 ) + printf( "0: mesh x: %d\n", mesh->getDimensions().x() ); + printf( "0: mesh y: %d\n", mesh->getDimensions().y() ); + //aux.save("aux_proc0.tnl"); + /*for( int k = 0; k < mesh->getDimensions().x()*mesh->getDimensions().y(); k++ ) + aux[ k ] = 10; + for( int k = 0; k < mesh->getDimensions().x(); k++ ){ + for( int l = 0; l < mesh->getDimensions().y(); l++ ) + printf("%f.2\t",aux[ k * 16 + l ] ); + printf("\n"); + }*/ + + /*bool a = Communicators::MpiCommunicator::IsInitialized(); if( a ) printf("Je Init\n"); else - printf("Neni Init\n"); -//#endif + printf("Neni Init\n");*/ +#endif while( iteration < this->maxIterations ) {