Fixing setup of grid subdomain overlaps for periodic boundary conditions.

There is no overlap if the periodic boundary conditions are applied within one MPI process.

src/TNL/Meshes/DistributedMeshes/SubdomainOverlapsGetter.h

@@ -46,15 +46,49 @@ class SubdomainOverlapsGetter< Mesh< MeshConfig, Device >, Communicator >
                                const SubdomainOverlapsType& periodicBoundariesOverlapSize = 0 );
 };
 
-template< int Dimension,
-          typename Real,
+// TODO: Specializations by the grid dimension can be avoided when the MPI directions are 
+// rewritten in a dimension independent way
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename Communicator >
+class SubdomainOverlapsGetter< Grid< 1, Real, Device, Index >, Communicator >
+{
+   public:
+      
+      static const int Dimension = 1;
+      using MeshType = Grid< Dimension, Real, Device, Index >;
+      using DeviceType = Device;
+      using IndexType = Index;
+      using DistributedMeshType = DistributedMesh< MeshType >;
+      using SubdomainOverlapsType = typename DistributedMeshType::SubdomainOverlapsType;
+      using CoordinatesType = typename DistributedMeshType::CoordinatesType;
+      using CommunicatorType = Communicator;
+      
+      // Computes subdomain overlaps
+      /* SubdomainOverlapsType is a touple of the same size as the mesh dimensions. 
+       * lower.x() is overlap of the subdomain at boundary where x = 0,
+       * upper.x() is overlap of the subdomain at boundary where x = grid.getDimensions().x() - 1,
+       */
+      static void getOverlaps( const DistributedMeshType* distributedMesh,
+                               SubdomainOverlapsType& lower,
+                               SubdomainOverlapsType& upper,
+                               IndexType subdomainOverlapSize,
+                               const SubdomainOverlapsType& periodicBoundariesOverlapSize = 0 );
+   
+};
+
+
+template< typename Real,
           typename Device,
           typename Index,
           typename Communicator >
-class SubdomainOverlapsGetter< Grid< Dimension, Real, Device, Index >, Communicator >
+class SubdomainOverlapsGetter< Grid< 2, Real, Device, Index >, Communicator >
 {
    public:
       
+      static const int Dimension = 2;
       using MeshType = Grid< Dimension, Real, Device, Index >;
       using DeviceType = Device;
       using IndexType = Index;
@@ -66,9 +100,9 @@ class SubdomainOverlapsGetter< Grid< Dimension, Real, Device, Index >, Communica
       // Computes subdomain overlaps
       /* SubdomainOverlapsType is a touple of the same size as the mesh dimensions. 
        * lower.x() is overlap of the subdomain at boundary where x = 0,
-       * lower.y() is overlap of the subdomain at boundary where y = 0 etc.
+       * lower.y() is overlap of the subdomain at boundary where y = 0,
        * upper.x() is overlap of the subdomain at boundary where x = grid.getDimensions().x() - 1,
-       * upper.y() is overlap of the subdomain at boundary where y = grid.getDimensions().y() - 1 etc.
+       * upper.y() is overlap of the subdomain at boundary where y = grid.getDimensions().y() - 1.
        */
       static void getOverlaps( const DistributedMeshType* distributedMesh,
                                SubdomainOverlapsType& lower,
@@ -78,6 +112,41 @@ class SubdomainOverlapsGetter< Grid< Dimension, Real, Device, Index >, Communica
    
 };
 
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename Communicator >
+class SubdomainOverlapsGetter< Grid< 3, Real, Device, Index >, Communicator >
+{
+   public:
+      
+      static const int Dimension = 3;
+      using MeshType = Grid< Dimension, Real, Device, Index >;
+      using DeviceType = Device;
+      using IndexType = Index;
+      using DistributedMeshType = DistributedMesh< MeshType >;
+      using SubdomainOverlapsType = typename DistributedMeshType::SubdomainOverlapsType;
+      using CoordinatesType = typename DistributedMeshType::CoordinatesType;
+      using CommunicatorType = Communicator;
+      
+      // Computes subdomain overlaps
+      /* SubdomainOverlapsType is a touple of the same size as the mesh dimensions. 
+       * lower.x() is overlap of the subdomain at boundary where x = 0,
+       * lower.y() is overlap of the subdomain at boundary where y = 0,
+       * lower.z() is overlap of the subdomain at boundary where z = 0,
+       * upper.x() is overlap of the subdomain at boundary where x = grid.getDimensions().x() - 1,
+       * upper.y() is overlap of the subdomain at boundary where y = grid.getDimensions().y() - 1,
+       * upper.z() is overlap of the subdomain at boundary where z = grid.getDimensions().z() - 1,
+       */
+      static void getOverlaps( const DistributedMeshType* distributedMesh,
+                               SubdomainOverlapsType& lower,
+                               SubdomainOverlapsType& upper,
+                               IndexType subdomainOverlapSize,
+                               const SubdomainOverlapsType& periodicBoundariesOverlapSize = 0 );
+   
+};
+
+
       } // namespace DistributedMeshes
    } // namespace Meshes
 } // namespace TNL

src/TNL/Meshes/DistributedMeshes/SubdomainOverlapsGetter.hpp

@@ -17,13 +17,15 @@ namespace TNL {
    namespace Meshes {
       namespace DistributedMeshes {
 
-template< int Dimension,
-          typename Real,
+/*
+ * TODO: This could work when the MPI directions are rewritten
+         
+template< typename Real,
           typename Device,
           typename Index,
           typename Communicator >
 void
-SubdomainOverlapsGetter< Grid< Dimension, Real, Device, Index >, Communicator >::
+SubdomainOverlapsGetter< Grid< 1, Real, Device, Index >, Communicator >::
 getOverlaps( const DistributedMeshType* distributedMesh,
              SubdomainOverlapsType& lower,
              SubdomainOverlapsType& upper,
@@ -33,7 +35,7 @@ getOverlaps( const DistributedMeshType* distributedMesh,
    if( ! CommunicatorType::isDistributed() )
       return;
    TNL_ASSERT_TRUE( distributedMesh != NULL, "" );
-   
+
    const CoordinatesType& subdomainCoordinates = distributedMesh->getSubdomainCoordinates();
    int rank = CommunicatorType::GetRank( CommunicatorType::AllGroup );
    
@@ -53,6 +55,129 @@ getOverlaps( const DistributedMeshType* distributedMesh,
          upper[ i ] = periodicBoundariesOverlapSize[ i ];
    }
 }
+ 
+*/
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename Communicator >
+void
+SubdomainOverlapsGetter< Grid< 1, Real, Device, Index >, Communicator >::
+getOverlaps( const DistributedMeshType* distributedMesh,
+             SubdomainOverlapsType& lower,
+             SubdomainOverlapsType& upper,
+             IndexType subdomainOverlapSize,
+             const SubdomainOverlapsType& periodicBoundariesOverlapSize )
+{
+   if( ! CommunicatorType::isDistributed() )
+      return;
+   TNL_ASSERT_TRUE( distributedMesh != NULL, "" );
+
+   const CoordinatesType& subdomainCoordinates = distributedMesh->getSubdomainCoordinates();
+   int rank = CommunicatorType::GetRank( CommunicatorType::AllGroup );
+   
+   if( subdomainCoordinates[ 0 ] > 0 )
+      lower[ 0 ] = subdomainOverlapSize;
+   else if( distributedMesh->getPeriodicNeighbors()[ ZzYzXm ] != rank )
+      lower[ 0 ] = periodicBoundariesOverlapSize[ 0 ];
+
+   if( subdomainCoordinates[ 0 ] < distributedMesh->getDomainDecomposition()[ 0 ] - 1 )
+      upper[ 0 ] = subdomainOverlapSize;
+   else if( distributedMesh->getPeriodicNeighbors()[ ZzYzXp ] != rank )
+      upper[ 0 ] = periodicBoundariesOverlapSize[ 0 ];
+}
+
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename Communicator >
+void
+SubdomainOverlapsGetter< Grid< 2, Real, Device, Index >, Communicator >::
+getOverlaps( const DistributedMeshType* distributedMesh,
+             SubdomainOverlapsType& lower,
+             SubdomainOverlapsType& upper,
+             IndexType subdomainOverlapSize,
+             const SubdomainOverlapsType& periodicBoundariesOverlapSize )
+{
+   if( ! CommunicatorType::isDistributed() )
+      return;
+   TNL_ASSERT_TRUE( distributedMesh != NULL, "" );
+
+   const CoordinatesType& subdomainCoordinates = distributedMesh->getSubdomainCoordinates();
+   int rank = CommunicatorType::GetRank( CommunicatorType::AllGroup );
+   
+   if( subdomainCoordinates[ 0 ] > 0 )
+      lower[ 0 ] = subdomainOverlapSize;
+   else if( distributedMesh->getPeriodicNeighbors()[ ZzYzXm ] != rank )
+      lower[ 0 ] = periodicBoundariesOverlapSize[ 0 ];
+
+   if( subdomainCoordinates[ 0 ] < distributedMesh->getDomainDecomposition()[ 0 ] - 1 )
+      upper[ 0 ] = subdomainOverlapSize;
+   else if( distributedMesh->getPeriodicNeighbors()[ ZzYzXp ] != rank )
+      upper[ 0 ] = periodicBoundariesOverlapSize[ 0 ];
+   
+   if( subdomainCoordinates[ 1 ] > 0 )
+      lower[ 1 ] = subdomainOverlapSize;
+   else if( distributedMesh->getPeriodicNeighbors()[ ZzYmXz ] != rank )
+      lower[ 1 ] = periodicBoundariesOverlapSize[ 1 ];
+
+   if( subdomainCoordinates[ 1 ] < distributedMesh->getDomainDecomposition()[ 1 ] - 1 )
+      upper[ 1 ] = subdomainOverlapSize;
+   else if( distributedMesh->getPeriodicNeighbors()[ ZzYpXz ] != rank )
+      upper[ 1 ] = periodicBoundariesOverlapSize[ 1 ];
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename Communicator >
+void
+SubdomainOverlapsGetter< Grid< 3, Real, Device, Index >, Communicator >::
+getOverlaps( const DistributedMeshType* distributedMesh,
+             SubdomainOverlapsType& lower,
+             SubdomainOverlapsType& upper,
+             IndexType subdomainOverlapSize,
+             const SubdomainOverlapsType& periodicBoundariesOverlapSize )
+{
+   if( ! CommunicatorType::isDistributed() )
+      return;
+   TNL_ASSERT_TRUE( distributedMesh != NULL, "" );
+
+   const CoordinatesType& subdomainCoordinates = distributedMesh->getSubdomainCoordinates();
+   int rank = CommunicatorType::GetRank( CommunicatorType::AllGroup );
+   
+   if( subdomainCoordinates[ 0 ] > 0 )
+      lower[ 0 ] = subdomainOverlapSize;
+   else if( distributedMesh->getPeriodicNeighbors()[ ZzYzXm ] != rank )
+      lower[ 0 ] = periodicBoundariesOverlapSize[ 0 ];
+
+   if( subdomainCoordinates[ 0 ] < distributedMesh->getDomainDecomposition()[ 0 ] - 1 )
+      upper[ 0 ] = subdomainOverlapSize;
+   else if( distributedMesh->getPeriodicNeighbors()[ ZzYzXp ] != rank )
+      upper[ 0 ] = periodicBoundariesOverlapSize[ 0 ];
+   
+   if( subdomainCoordinates[ 1 ] > 0 )
+      lower[ 1 ] = subdomainOverlapSize;
+   else if( distributedMesh->getPeriodicNeighbors()[ ZzYmXz ] != rank )
+      lower[ 1 ] = periodicBoundariesOverlapSize[ 1 ];
+
+   if( subdomainCoordinates[ 1 ] < distributedMesh->getDomainDecomposition()[ 1 ] - 1 )
+      upper[ 1 ] = subdomainOverlapSize;
+   else if( distributedMesh->getPeriodicNeighbors()[ ZzYpXz ] != rank )
+      upper[ 1 ] = periodicBoundariesOverlapSize[ 1 ];
+   
+   if( subdomainCoordinates[ 2 ] > 0 )
+      lower[ 2 ] = subdomainOverlapSize;
+   else if( distributedMesh->getPeriodicNeighbors()[ ZmYzXz ] != rank )
+      lower[ 2 ] = periodicBoundariesOverlapSize[ 2 ];
+
+   if( subdomainCoordinates[ 2 ] < distributedMesh->getDomainDecomposition()[ 2 ] - 1 )
+      upper[ 2 ] = subdomainOverlapSize;
+   else if( distributedMesh->getPeriodicNeighbors()[ ZpYzXz ] != rank )
+      upper[ 2 ] = periodicBoundariesOverlapSize[ 2 ];
+}
 
       } // namespace DistributedMeshes
    } // namespace Meshes