MeshUtils.hpp

Go to the documentation of this file.

/*
 * ------------------------------------------------------------------------------------------------------------
 * SPDX-License-Identifier: LGPL-2.1-only
 *
 * Copyright (c) 2018-2019 Lawrence Livermore National Security LLC
 * Copyright (c) 2018-2019 The Board of Trustees of the Leland Stanford Junior University
 * Copyright (c) 2018-2019 Total, S.A
 * Copyright (c) 2019-     GEOS/GEOSX Contributors
 * All right reserved
 *
 * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
 * ------------------------------------------------------------------------------------------------------------
 */
 
#ifndef GEOS_LINEARALGEBRA_MULTISCALE_MESHUTILS_HPP
#define GEOS_LINEARALGEBRA_MULTISCALE_MESHUTILS_HPP
 
#include "common/DataTypes.hpp"
#include "linearAlgebra/multiscale/mesh/MeshObjectManager.hpp"
#include "mesh/ObjectManagerBase.hpp"
 
namespace geos
{
 
class DomainPartition;
 
namespace multiscale
{
namespace meshUtils
{
 
class ScopedDataRegistrar
{
public:
 
  template< typename T >
  ScopedDataRegistrar( ObjectManagerBase & manager, string key, T & data )
    : m_manager( manager ),
    m_key( std::move( key ) )
  {
    m_manager.registerWrapper( m_key, &data );
  }
 
  ScopedDataRegistrar( ScopedDataRegistrar const & ) = delete;
  ScopedDataRegistrar( ScopedDataRegistrar && ) = delete;
 
  ~ScopedDataRegistrar()
  {
    m_manager.deregisterWrapper( m_key );
  }
 
  void sync( DomainPartition & domain ) const;
 
private:
 
  ObjectManagerBase & m_manager;
  string const m_key;
};
 
template< typename T, typename U = T >
void mapIndexArray( arrayView1d< T const > const & src,
                    arrayView1d< U const > const & map,
                    array1d< U > & dst )
{
  dst.clear();
  dst.reserve( src.size() );
  for( T const & val : src )
  {
    U const newVal = map[val];
    if( newVal >= 0 )
    {
      dst.emplace_back( newVal );
    }
  }
}
 
template< typename T, typename U = T >
void mapIndexArrayUnique( arrayView1d< T const > const & src,
                          arrayView1d< U const > const & map,
                          array1d< U > & dst )
{
  SortedArray< U > values;
  values.reserve( src.size() );
  for( T const & val : src )
  {
    U const newVal = map[val];
    if( newVal >= 0 )
    {
      values.insert( newVal );
    }
  }
  dst.clear();
  dst.reserve( values.size() );
  for( U const & val : values )
  {
    dst.emplace_back( val );
  }
}
 
template< typename T, typename U = T >
void mapIndexSet( SortedArrayView< T const > const & src,
                  arrayView1d< U const > const & map,
                  SortedArray< U > & dst )
{
  dst.reserve( src.size() );
  for( T const & val : src )
  {
    U const newVal = map[val];
    if( newVal >= 0 )
    {
      dst.insert( newVal );
    }
  }
}
 
template< typename POLICY, typename T, int NDIM, int USD, typename INDEX, typename FUNC >
void fillArrayByDstIndex( arrayView1d< INDEX const > const & map,
                          ArrayView< T, NDIM, USD > const & dst,
                          FUNC && src )
{
  forAll< POLICY >( map.size(), [=]( INDEX const srcIdx )
  {
    INDEX const dstIdx = map[srcIdx];
    if( dstIdx >= 0 )
    {
      LvArray::forValuesInSliceWithIndices( dst[dstIdx], [&]( T & dstValue, auto ... indices )
      {
        dstValue = src( srcIdx, indices ... );
      } );
    }
  } );
}
 
template< typename POLICY, typename T, int NDIM, int USD, typename INDEX, typename FUNC >
void fillArrayBySrcIndex( arrayView1d< INDEX const > const & map,
                          ArrayView< T, NDIM, USD > const & dst,
                          FUNC && src )
{
  GEOS_ASSERT_EQ( dst.size(), map.size() );
  forAll< POLICY >( map.size(), [=]( INDEX const dstIndex )
  {
    INDEX const srcIndex = map[dstIndex];
    if( srcIndex >= 0 )
    {
      LvArray::forValuesInSliceWithIndices( dst[dstIndex], [&]( T & dstValue, auto ... indices )
      {
        dstValue = src( srcIndex, indices ... );
      } );
    }
  } );
}
 
template< typename FUNC >
void copyNeighborData( ObjectManagerBase const & srcManager,
                       string const & mapKey,
                       std::vector< integer > const & ranks,
                       ObjectManagerBase & dstManager,
                       FUNC && copyFunc )
{
  arrayView1d< localIndex const > const map = srcManager.getReference< array1d< localIndex > >( mapKey );
  for( integer const rank : ranks )
  {
    NeighborData const & srcData = srcManager.getNeighborData( rank );
    NeighborData & dstData = dstManager.getNeighborData( rank );
    copyFunc( srcData.ghostsToSend(), map, dstData.ghostsToSend() );
    copyFunc( srcData.ghostsToReceive(), map, dstData.ghostsToReceive() );
    copyFunc( srcData.adjacencyList(), map, dstData.adjacencyList() );
    copyFunc( srcData.matchedPartitionBoundary(), map, dstData.matchedPartitionBoundary() );
  }
}
 
void copySets( ObjectManagerBase const & srcManager,
               string const & mapKey,
               ObjectManagerBase & dstManager );
 
namespace internal
{
 
// We don't have versions of IS_VALID_EXPRESSION for more args (and I don't want to add them).
// So second argument (count) is hardcoded to std::ptrdiff_t below.
 
IS_VALID_EXPRESSION( isCallableWithoutArgs, T, std::declval< T >()() );
IS_VALID_EXPRESSION_2( isCallableWithArg, T, U, std::declval< T >()( std::declval< U >() ) );
IS_VALID_EXPRESSION_2( isCallableWithArgAndCount, T, U, std::declval< T >()( std::declval< U >(), std::ptrdiff_t{} ) );
 
template< typename T, typename FUNC >
std::enable_if_t< isCallableWithoutArgs< FUNC > >
callWithArgs( T const & val, std::ptrdiff_t const count, FUNC && func )
{
  GEOS_UNUSED_VAR( val, count );
  func();
}
 
template< typename T, typename FUNC >
std::enable_if_t< isCallableWithArg< FUNC, T > >
callWithArgs( T const & val, std::ptrdiff_t const count, FUNC && func )
{
  GEOS_UNUSED_VAR( count );
  func( val );
}
 
template< typename T, typename FUNC >
std::enable_if_t< isCallableWithArgAndCount< FUNC, T > >
callWithArgs( T const & val, std::ptrdiff_t const count, FUNC && func )
{
  func( val, count );
}
 
} // namespace internal
 
template< typename ITER, typename FUNC >
void forUniqueValues( ITER first, ITER const last, FUNC && func )
{
  if( first == last )
    return;
  LvArray::sortedArrayManipulation::makeSorted( first, last );
  using T = typename std::iterator_traits< ITER >::value_type;
  while( first != last )
  {
    T const & curr = *first;
    ITER const it = std::find_if( first, last, [&curr]( T const & v ) { return v != curr; } );
    internal::callWithArgs( curr, std::distance( first, it ), std::forward< FUNC >( func ) );
    first = it;
  }
}
 
template< integer MAX_NEIGHBORS, typename L2C_MAP, typename C2L_MAP, typename PRED, typename FUNC >
void forUniqueNeighbors( localIndex const locIdx,
                         L2C_MAP const & locToConn,
                         C2L_MAP const & connToLoc,
                         PRED && connPred,
                         FUNC && func )
{
  localIndex neighbors[MAX_NEIGHBORS];
  integer numNeighbors = 0;
  for( localIndex const connIdx : locToConn[locIdx] )
  {
    if( connIdx >= 0 && connPred( connIdx ) )
    {
      for( localIndex const nbrIdx : connToLoc[connIdx] )
      {
        GEOS_ERROR_IF_GE_MSG( numNeighbors, MAX_NEIGHBORS, "Too many neighbors, need to increase stack limit" );
        neighbors[numNeighbors++] = nbrIdx;
      }
    }
  }
  forUniqueValues( neighbors, neighbors + numNeighbors, std::forward< FUNC >( func ) );
}
 
template< integer MAX_NEIGHBORS, typename L2C_MAP, typename C2L_MAP, typename FUNC >
void forUniqueNeighbors( localIndex const locIdx,
                         L2C_MAP const & locToConn,
                         C2L_MAP const & connToLoc,
                         FUNC && func )
{
  forUniqueNeighbors< MAX_NEIGHBORS >( locIdx,
                                       locToConn,
                                       connToLoc,
                                       []( auto ){ return true; },
                                       std::forward< FUNC >( func ) );
}
 
template< integer MAX_NEIGHBORS, typename NBR_MAP, typename VAL_FUNC, typename VAL_PRED, typename FUNC >
void forUniqueNeighborValues( localIndex const locIdx,
                              NBR_MAP const & neighbors,
                              VAL_FUNC && valueFunc,
                              VAL_PRED && pred,
                              FUNC && func )
{
  using T = std::remove_cv_t< std::remove_reference_t< decltype( valueFunc( localIndex {} ) ) >>;
  T nbrValues[MAX_NEIGHBORS];
  integer numValues = 0;
  for( localIndex const nbrIdx : neighbors[locIdx] )
  {
    GEOS_ERROR_IF_GE_MSG( numValues, MAX_NEIGHBORS, "Too many neighbors, need to increase stack limit" );
    T const value = valueFunc( nbrIdx );
    if( pred( value ) )
    {
      nbrValues[numValues++] = value;
    }
  }
  forUniqueValues( nbrValues, nbrValues + numValues, std::forward< FUNC >( func ) );
}
 
template< integer MAX_NEIGHBORS, typename NBR_MAP, typename VAL_FUNC, typename FUNC >
void forUniqueNeighborValues( localIndex const locIdx,
                              NBR_MAP const & neighbors,
                              VAL_FUNC && valueFunc,
                              FUNC && func )
{
  forUniqueNeighborValues< MAX_NEIGHBORS >( locIdx,
                                            neighbors,
                                            std::forward< VAL_FUNC >( valueFunc ),
                                            []( auto ){ return true; },
                                            std::forward< FUNC >( func ) );
}
 
template< typename INDEX_TYPE >
ArrayOfSets< INDEX_TYPE >
buildFineObjectToSubdomainMap( MeshObjectManager const & fineObjectManager,
                               arrayView1d< INDEX_TYPE const > const & subdomains,
                               string_array const & boundaryObjectSets )
{
  GEOS_MARK_FUNCTION;
 
  MeshObjectManager::MapViewConst const dualMap = fineObjectManager.toDualRelation().toViewConst();
 
  // count the row lengths
  array1d< localIndex > rowCounts( fineObjectManager.size() );
  forAll< parallelHostPolicy >( fineObjectManager.size(), [=, rowCounts = rowCounts.toView()]( localIndex const objIdx )
  {
    localIndex count = 0;
    forUniqueNeighborValues< 256 >( objIdx, dualMap, subdomains, [&count]
    {
      ++count;
    } );
    rowCounts[objIdx] = count;
  } );
  for( string const & setName : boundaryObjectSets )
  {
    SortedArrayView< localIndex const > const set = fineObjectManager.getSet( setName ).toViewConst();
    forAll< parallelHostPolicy >( set.size(), [=, rowCounts = rowCounts.toView()]( localIndex const i )
    {
      ++rowCounts[set[i]];
    } );
  }
 
  // Resize from row lengths
  ArrayOfSets< INDEX_TYPE > objectToSubdomain;
  objectToSubdomain.template resizeFromCapacities< parallelHostPolicy >( rowCounts.size(), rowCounts.data() );
 
  // Fill the map
  INDEX_TYPE numBoundaries = 0;
  for( string const & setName : boundaryObjectSets )
  {
    ++numBoundaries;
    SortedArrayView< localIndex const > const set = fineObjectManager.getSet( setName ).toViewConst();
    forAll< parallelHostPolicy >( set.size(), [=, objToSub = objectToSubdomain.toView()]( localIndex const i )
    {
      objToSub.insertIntoSet( set[i], -numBoundaries );
    } );
  }
  forAll< parallelHostPolicy >( fineObjectManager.size(), [=, objToSub = objectToSubdomain.toView()]( localIndex const objIdx )
  {
    for( localIndex const dualIdx : dualMap[objIdx] )
    {
      objToSub.insertIntoSet( objIdx, subdomains[dualIdx] );
    }
  } );
 
  return objectToSubdomain;
}
 
template< typename INDEX_TYPE >
ArrayOfSets< INDEX_TYPE >
addBoundarySubdomains( MeshObjectManager const & manager,
                       ArrayOfSetsView< INDEX_TYPE const > const & inputMap,
                       string_array const & boundaryObjectSets )
{
  GEOS_MARK_FUNCTION;
 
  array1d< localIndex > rowCounts( manager.size() );
  forAll< parallelHostPolicy >( manager.size(), [=, rowCounts = rowCounts.toView()]( localIndex const objIdx )
  {
    rowCounts[objIdx] = inputMap.sizeOfSet( objIdx );
  } );
  for( string const & setName: boundaryObjectSets )
  {
    SortedArrayView< localIndex const > const set = manager.getSet( setName ).toViewConst();
    forAll< parallelHostPolicy >( set.size(), [=, rowCounts = rowCounts.toView()]( localIndex const i )
    {
      ++rowCounts[set[i]];
    } );
  }
 
  // Resize from row lengths
  ArrayOfSets< INDEX_TYPE > outputMap;
  outputMap.template resizeFromCapacities< parallelHostPolicy >( rowCounts.size(), rowCounts.data() );
 
  // Fill the map
  INDEX_TYPE numBoundaries = 0;
  for( string const & setName: boundaryObjectSets )
  {
    ++numBoundaries;
    SortedArrayView< localIndex const > const set = manager.getSet( setName ).toViewConst();
    forAll< parallelHostPolicy >( set.size(), [=, outputMap = outputMap.toView()]( localIndex const i )
    {
      outputMap.insertIntoSet( set[i], -numBoundaries );
    } );
  }
  forAll< parallelHostPolicy >( manager.size(), [=, outputMap = outputMap.toView()]( localIndex const objIdx )
  {
    outputMap.insertIntoSet( objIdx, inputMap[objIdx].begin(), inputMap[objIdx].end() );
  } );
 
  return outputMap;
}
 
array1d< localIndex >
findCoarseNodesByDualPartition( MeshObjectManager::MapViewConst const & nodeToDual,
                                MeshObjectManager::MapViewConst const & dualToNode,
                                ArrayOfSetsView< globalIndex const > const & nodeToSubdomain,
                                integer const minSubdomains,
                                bool allowMultiNodes = true );
 
} // namespace meshUtils
} // namespace multiscale
} // namespace geos
 
#endif //GEOS_LINEARALGEBRA_MULTISCALE_MESHUTILS_HPP