GEOS
MeshUtils.hpp
Go to the documentation of this file.
1 /*
2  * ------------------------------------------------------------------------------------------------------------
3  * SPDX-License-Identifier: LGPL-2.1-only
4  *
5  * Copyright (c) 2018-2019 Lawrence Livermore National Security LLC
6  * Copyright (c) 2018-2019 The Board of Trustees of the Leland Stanford Junior University
7  * Copyright (c) 2018-2019 Total, S.A
8  * Copyright (c) 2019- GEOS/GEOSX Contributors
9  * All right reserved
10  *
11  * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
12  * ------------------------------------------------------------------------------------------------------------
13  */
14 
19 #ifndef GEOS_LINEARALGEBRA_MULTISCALE_MESHUTILS_HPP
20 #define GEOS_LINEARALGEBRA_MULTISCALE_MESHUTILS_HPP
21 
22 #include "common/DataTypes.hpp"
25 
26 namespace geos
27 {
28 
29 class DomainPartition;
30 
31 namespace multiscale
32 {
33 namespace meshUtils
34 {
35 
45 {
46 public:
47 
55  template< typename T >
56  ScopedDataRegistrar( ObjectManagerBase & manager, string key, T & data )
57  : m_manager( manager ),
58  m_key( std::move( key ) )
59  {
60  m_manager.registerWrapper( m_key, &data );
61  }
62 
63  ScopedDataRegistrar( ScopedDataRegistrar const & ) = delete;
65 
70  {
71  m_manager.deregisterWrapper( m_key );
72  }
73 
78  void sync( DomainPartition & domain ) const;
79 
80 private:
81 
82  ObjectManagerBase & m_manager;
83  string const m_key;
84 };
85 
97 template< typename T, typename U = T >
100  array1d< U > & dst )
101 {
102  dst.clear();
103  dst.reserve( src.size() );
104  for( T const & val : src )
105  {
106  U const newVal = map[val];
107  if( newVal >= 0 )
108  {
109  dst.emplace_back( newVal );
110  }
111  }
112 }
113 
126 template< typename T, typename U = T >
128  arrayView1d< U const > const & map,
129  array1d< U > & dst )
130 {
131  SortedArray< U > values;
132  values.reserve( src.size() );
133  for( T const & val : src )
134  {
135  U const newVal = map[val];
136  if( newVal >= 0 )
137  {
138  values.insert( newVal );
139  }
140  }
141  dst.clear();
142  dst.reserve( values.size() );
143  for( U const & val : values )
144  {
145  dst.emplace_back( val );
146  }
147 }
148 
157 template< typename T, typename U = T >
159  arrayView1d< U const > const & map,
160  SortedArray< U > & dst )
161 {
162  dst.reserve( src.size() );
163  for( T const & val : src )
164  {
165  U const newVal = map[val];
166  if( newVal >= 0 )
167  {
168  dst.insert( newVal );
169  }
170  }
171 }
172 
187 template< typename POLICY, typename T, int NDIM, int USD, typename INDEX, typename FUNC >
189  ArrayView< T, NDIM, USD > const & dst,
190  FUNC && src )
191 {
192  forAll< POLICY >( map.size(), [=]( INDEX const srcIdx )
193  {
194  INDEX const dstIdx = map[srcIdx];
195  if( dstIdx >= 0 )
196  {
197  LvArray::forValuesInSliceWithIndices( dst[dstIdx], [&]( T & dstValue, auto ... indices )
198  {
199  dstValue = src( srcIdx, indices ... );
200  } );
201  }
202  } );
203 }
204 
219 template< typename POLICY, typename T, int NDIM, int USD, typename INDEX, typename FUNC >
221  ArrayView< T, NDIM, USD > const & dst,
222  FUNC && src )
223 {
224  GEOS_ASSERT_EQ( dst.size(), map.size() );
225  forAll< POLICY >( map.size(), [=]( INDEX const dstIndex )
226  {
227  INDEX const srcIndex = map[dstIndex];
228  if( srcIndex >= 0 )
229  {
230  LvArray::forValuesInSliceWithIndices( dst[dstIndex], [&]( T & dstValue, auto ... indices )
231  {
232  dstValue = src( srcIndex, indices ... );
233  } );
234  }
235  } );
236 }
237 
247 template< typename FUNC >
248 void copyNeighborData( ObjectManagerBase const & srcManager,
249  string const & mapKey,
250  std::vector< integer > const & ranks,
251  ObjectManagerBase & dstManager,
252  FUNC && copyFunc )
253 {
255  for( integer const rank : ranks )
256  {
257  NeighborData const & srcData = srcManager.getNeighborData( rank );
258  NeighborData & dstData = dstManager.getNeighborData( rank );
259  copyFunc( srcData.ghostsToSend(), map, dstData.ghostsToSend() );
260  copyFunc( srcData.ghostsToReceive(), map, dstData.ghostsToReceive() );
261  copyFunc( srcData.adjacencyList(), map, dstData.adjacencyList() );
262  copyFunc( srcData.matchedPartitionBoundary(), map, dstData.matchedPartitionBoundary() );
263  }
264 }
265 
272 void copySets( ObjectManagerBase const & srcManager,
273  string const & mapKey,
274  ObjectManagerBase & dstManager );
275 
276 namespace internal
277 {
278 
279 // We don't have versions of IS_VALID_EXPRESSION for more args (and I don't want to add them).
280 // So second argument (count) is hardcoded to std::ptrdiff_t below.
281 
282 IS_VALID_EXPRESSION( isCallableWithoutArgs, T, std::declval< T >()() );
283 IS_VALID_EXPRESSION_2( isCallableWithArg, T, U, std::declval< T >()( std::declval< U >() ) );
284 IS_VALID_EXPRESSION_2( isCallableWithArgAndCount, T, U, std::declval< T >()( std::declval< U >(), std::ptrdiff_t{} ) );
285 
286 template< typename T, typename FUNC >
287 std::enable_if_t< isCallableWithoutArgs< FUNC > >
288 callWithArgs( T const & val, std::ptrdiff_t const count, FUNC && func )
289 {
290  GEOS_UNUSED_VAR( val, count );
291  func();
292 }
293 
294 template< typename T, typename FUNC >
295 std::enable_if_t< isCallableWithArg< FUNC, T > >
296 callWithArgs( T const & val, std::ptrdiff_t const count, FUNC && func )
297 {
298  GEOS_UNUSED_VAR( count );
299  func( val );
300 }
301 
302 template< typename T, typename FUNC >
303 std::enable_if_t< isCallableWithArgAndCount< FUNC, T > >
304 callWithArgs( T const & val, std::ptrdiff_t const count, FUNC && func )
305 {
306  func( val, count );
307 }
308 
309 } // namespace internal
310 
320 template< typename ITER, typename FUNC >
321 void forUniqueValues( ITER first, ITER const last, FUNC && func )
322 {
323  if( first == last )
324  return;
325  LvArray::sortedArrayManipulation::makeSorted( first, last );
326  using T = typename std::iterator_traits< ITER >::value_type;
327  while( first != last )
328  {
329  T const & curr = *first;
330  ITER const it = std::find_if( first, last, [&curr]( T const & v ) { return v != curr; } );
331  internal::callWithArgs( curr, std::distance( first, it ), std::forward< FUNC >( func ) );
332  first = it;
333  }
334 }
335 
349 template< integer MAX_NEIGHBORS, typename L2C_MAP, typename C2L_MAP, typename PRED, typename FUNC >
350 void forUniqueNeighbors( localIndex const locIdx,
351  L2C_MAP const & locToConn,
352  C2L_MAP const & connToLoc,
353  PRED && connPred,
354  FUNC && func )
355 {
356  localIndex neighbors[MAX_NEIGHBORS];
357  integer numNeighbors = 0;
358  for( localIndex const connIdx : locToConn[locIdx] )
359  {
360  if( connIdx >= 0 && connPred( connIdx ) )
361  {
362  for( localIndex const nbrIdx : connToLoc[connIdx] )
363  {
364  GEOS_ERROR_IF_GE_MSG( numNeighbors, MAX_NEIGHBORS, "Too many neighbors, need to increase stack limit" );
365  neighbors[numNeighbors++] = nbrIdx;
366  }
367  }
368  }
369  forUniqueValues( neighbors, neighbors + numNeighbors, std::forward< FUNC >( func ) );
370 }
371 
385 template< integer MAX_NEIGHBORS, typename L2C_MAP, typename C2L_MAP, typename FUNC >
386 void forUniqueNeighbors( localIndex const locIdx,
387  L2C_MAP const & locToConn,
388  C2L_MAP const & connToLoc,
389  FUNC && func )
390 {
391  forUniqueNeighbors< MAX_NEIGHBORS >( locIdx,
392  locToConn,
393  connToLoc,
394  []( auto ){ return true; },
395  std::forward< FUNC >( func ) );
396 }
397 
411 template< integer MAX_NEIGHBORS, typename NBR_MAP, typename VAL_FUNC, typename VAL_PRED, typename FUNC >
413  NBR_MAP const & neighbors,
414  VAL_FUNC && valueFunc,
415  VAL_PRED && pred,
416  FUNC && func )
417 {
418  using T = std::remove_cv_t< std::remove_reference_t< decltype( valueFunc( localIndex {} ) ) >>;
419  T nbrValues[MAX_NEIGHBORS];
420  integer numValues = 0;
421  for( localIndex const nbrIdx : neighbors[locIdx] )
422  {
423  GEOS_ERROR_IF_GE_MSG( numValues, MAX_NEIGHBORS, "Too many neighbors, need to increase stack limit" );
424  T const value = valueFunc( nbrIdx );
425  if( pred( value ) )
426  {
427  nbrValues[numValues++] = value;
428  }
429  }
430  forUniqueValues( nbrValues, nbrValues + numValues, std::forward< FUNC >( func ) );
431 }
432 
446 template< integer MAX_NEIGHBORS, typename NBR_MAP, typename VAL_FUNC, typename FUNC >
448  NBR_MAP const & neighbors,
449  VAL_FUNC && valueFunc,
450  FUNC && func )
451 {
452  forUniqueNeighborValues< MAX_NEIGHBORS >( locIdx,
453  neighbors,
454  std::forward< VAL_FUNC >( valueFunc ),
455  []( auto ){ return true; },
456  std::forward< FUNC >( func ) );
457 }
458 
471 template< typename INDEX_TYPE >
474  arrayView1d< INDEX_TYPE const > const & subdomains,
475  string_array const & boundaryObjectSets )
476 {
478 
479  MeshObjectManager::MapViewConst const dualMap = fineObjectManager.toDualRelation().toViewConst();
480 
481  // count the row lengths
482  array1d< localIndex > rowCounts( fineObjectManager.size() );
483  forAll< parallelHostPolicy >( fineObjectManager.size(), [=, rowCounts = rowCounts.toView()]( localIndex const objIdx )
484  {
485  localIndex count = 0;
486  forUniqueNeighborValues< 256 >( objIdx, dualMap, subdomains, [&count]
487  {
488  ++count;
489  } );
490  rowCounts[objIdx] = count;
491  } );
492  for( string const & setName : boundaryObjectSets )
493  {
494  SortedArrayView< localIndex const > const set = fineObjectManager.getSet( setName ).toViewConst();
495  forAll< parallelHostPolicy >( set.size(), [=, rowCounts = rowCounts.toView()]( localIndex const i )
496  {
497  ++rowCounts[set[i]];
498  } );
499  }
500 
501  // Resize from row lengths
502  ArrayOfSets< INDEX_TYPE > objectToSubdomain;
503  objectToSubdomain.template resizeFromCapacities< parallelHostPolicy >( rowCounts.size(), rowCounts.data() );
504 
505  // Fill the map
506  INDEX_TYPE numBoundaries = 0;
507  for( string const & setName : boundaryObjectSets )
508  {
509  ++numBoundaries;
510  SortedArrayView< localIndex const > const set = fineObjectManager.getSet( setName ).toViewConst();
511  forAll< parallelHostPolicy >( set.size(), [=, objToSub = objectToSubdomain.toView()]( localIndex const i )
512  {
513  objToSub.insertIntoSet( set[i], -numBoundaries );
514  } );
515  }
516  forAll< parallelHostPolicy >( fineObjectManager.size(), [=, objToSub = objectToSubdomain.toView()]( localIndex const objIdx )
517  {
518  for( localIndex const dualIdx : dualMap[objIdx] )
519  {
520  objToSub.insertIntoSet( objIdx, subdomains[dualIdx] );
521  }
522  } );
523 
524  return objectToSubdomain;
525 }
526 
539 template< typename INDEX_TYPE >
540 ArrayOfSets< INDEX_TYPE >
542  ArrayOfSetsView< INDEX_TYPE const > const & inputMap,
543  string_array const & boundaryObjectSets )
544 {
546 
547  array1d< localIndex > rowCounts( manager.size() );
548  forAll< parallelHostPolicy >( manager.size(), [=, rowCounts = rowCounts.toView()]( localIndex const objIdx )
549  {
550  rowCounts[objIdx] = inputMap.sizeOfSet( objIdx );
551  } );
552  for( string const & setName: boundaryObjectSets )
553  {
554  SortedArrayView< localIndex const > const set = manager.getSet( setName ).toViewConst();
555  forAll< parallelHostPolicy >( set.size(), [=, rowCounts = rowCounts.toView()]( localIndex const i )
556  {
557  ++rowCounts[set[i]];
558  } );
559  }
560 
561  // Resize from row lengths
562  ArrayOfSets< INDEX_TYPE > outputMap;
563  outputMap.template resizeFromCapacities< parallelHostPolicy >( rowCounts.size(), rowCounts.data() );
564 
565  // Fill the map
566  INDEX_TYPE numBoundaries = 0;
567  for( string const & setName: boundaryObjectSets )
568  {
569  ++numBoundaries;
570  SortedArrayView< localIndex const > const set = manager.getSet( setName ).toViewConst();
571  forAll< parallelHostPolicy >( set.size(), [=, outputMap = outputMap.toView()]( localIndex const i )
572  {
573  outputMap.insertIntoSet( set[i], -numBoundaries );
574  } );
575  }
576  forAll< parallelHostPolicy >( manager.size(), [=, outputMap = outputMap.toView()]( localIndex const objIdx )
577  {
578  outputMap.insertIntoSet( objIdx, inputMap[objIdx].begin(), inputMap[objIdx].end() );
579  } );
580 
581  return outputMap;
582 }
583 
607  MeshObjectManager::MapViewConst const & dualToNode,
608  ArrayOfSetsView< globalIndex const > const & nodeToSubdomain,
609  integer const minSubdomains,
610  bool allowMultiNodes = true );
611 
612 } // namespace meshUtils
613 } // namespace multiscale
614 } // namespace geos
615 
616 #endif //GEOS_LINEARALGEBRA_MULTISCALE_MESHUTILS_HPP
#define GEOS_UNUSED_VAR(...)
Mark an unused variable and silence compiler warnings.
Definition: GeosxMacros.hpp:84
#define GEOS_ASSERT_EQ(lhs, rhs)
Assert that two values compare equal in debug builds.
Definition: Logger.hpp:410
#define GEOS_ERROR_IF_GE_MSG(lhs, rhs, msg)
Raise a hard error if one value compares greater than or equal to the other.
Definition: Logger.hpp:307
void mapIndexSet(SortedArrayView< T const > const &src, arrayView1d< U const > const &map, SortedArray< U > &dst)
Transform a set of mesh indices using a map.
Definition: MeshUtils.hpp:158
void mapIndexArrayUnique(arrayView1d< T const > const &src, arrayView1d< U const > const &map, array1d< U > &dst)
Transform an array of mesh indices using a map, removing duplicates.
Definition: MeshUtils.hpp:127
void fillArrayBySrcIndex(arrayView1d< INDEX const > const &map, ArrayView< T, NDIM, USD > const &dst, FUNC &&src)
Populate an array using a function source and an index map.
Definition: MeshUtils.hpp:220
void fillArrayByDstIndex(arrayView1d< INDEX const > const &map, ArrayView< T, NDIM, USD > const &dst, FUNC &&src)
Populate an array using a function source and an index map.
Definition: MeshUtils.hpp:188
void forUniqueNeighbors(localIndex const locIdx, L2C_MAP const &locToConn, C2L_MAP const &connToLoc, PRED &&connPred, FUNC &&func)
Call a function on unique indices of topological neighbors visited through location-connector adjacen...
Definition: MeshUtils.hpp:350
void forUniqueNeighborValues(localIndex const locIdx, NBR_MAP const &neighbors, VAL_FUNC &&valueFunc, VAL_PRED &&pred, FUNC &&func)
Call a function on unique values of topological neighbors visited through location-connector adjacenc...
Definition: MeshUtils.hpp:412
ArrayOfSets< INDEX_TYPE > addBoundarySubdomains(MeshObjectManager const &manager, ArrayOfSetsView< INDEX_TYPE const > const &inputMap, string_array const &boundaryObjectSets)
Insert virtual boundary subdomains into an adjacency map.
Definition: MeshUtils.hpp:541
array1d< localIndex > findCoarseNodesByDualPartition(MeshObjectManager::MapViewConst const &nodeToDual, MeshObjectManager::MapViewConst const &dualToNode, ArrayOfSetsView< globalIndex const > const &nodeToSubdomain, integer const minSubdomains, bool allowMultiNodes=true)
Find "coarse" nodes in a mesh in which dual objects have been partitioned into subdomains.
void mapIndexArray(arrayView1d< T const > const &src, arrayView1d< U const > const &map, array1d< U > &dst)
Transform an array of mesh indices using a map.
Definition: MeshUtils.hpp:98
void copyNeighborData(ObjectManagerBase const &srcManager, string const &mapKey, std::vector< integer > const &ranks, ObjectManagerBase &dstManager, FUNC &&copyFunc)
Sets up NeighborData of a new object manager mapped from existing one.
Definition: MeshUtils.hpp:248
void forUniqueValues(ITER first, ITER const last, FUNC &&func)
Call a function on unique values from a previously collected range.
Definition: MeshUtils.hpp:321
ArrayOfSets< INDEX_TYPE > buildFineObjectToSubdomainMap(MeshObjectManager const &fineObjectManager, arrayView1d< INDEX_TYPE const > const &subdomains, string_array const &boundaryObjectSets)
Build a map from mesh objects (nodes/cells) to subdomains defined by a partitioning of dual objects (...
Definition: MeshUtils.hpp:473
void copySets(ObjectManagerBase const &srcManager, string const &mapKey, ObjectManagerBase &dstManager)
Copy and update sets from an existing to a new object manager.
#define GEOS_MARK_FUNCTION
Mark function with both Caliper and NVTX if enabled.
Partition of the decomposed physical domain. It also manages the connexion information to its neighbo...
The ObjectManagerBase is the base object of all object managers in the mesh data hierachy.
NeighborData & getNeighborData(int const rank)
Get neighbor data for given rank.
SortedArray< localIndex > & getSet(string const &setName)
Get a set by name.
Wrapper< TBASE > & registerWrapper(string const &name, wrapperMap::KeyIndex::index_type *const rkey=nullptr)
Create and register a Wrapper around a new object.
GEOS_DECLTYPE_AUTO_RETURN getReference(LOOKUP_TYPE const &lookup) const
Look up a wrapper and get reference to wrapped object.
Definition: Group.hpp:1275
void deregisterWrapper(string const &name)
Removes a Wrapper from this group.
localIndex size() const
Get the "size" of the group, which determines the number of elements in resizable wrappers.
Definition: Group.hpp:1317
Base template for ordered and unordered maps.
Mesh object manager used in multiscale preconditioners to keep a simplified (node/cell only) represen...
decltype(std::declval< MapType >().base().toViewConst()) MapViewConst
Alias for relation map const view type.
Utility class to manage registration of temporary scope-bound data in the data repository.
Definition: MeshUtils.hpp:45
ScopedDataRegistrar(ObjectManagerBase &manager, string key, T &data)
Constructor.
Definition: MeshUtils.hpp:56
void sync(DomainPartition &domain) const
Synchronize data across ranks.
ArrayView< T, 1 > arrayView1d
Alias for 1D array view.
Definition: DataTypes.hpp:179
stdVector< string > string_array
A 1-dimensional array of geos::string types.
Definition: DataTypes.hpp:361
LvArray::ArrayOfSetsView< T, INDEX_TYPE const, LvArray::ChaiBuffer > ArrayOfSetsView
View of array of variable-sized sets. See LvArray::ArrayOfSetsView for details.
Definition: DataTypes.hpp:293
std::set< T > set
A set of local indices.
Definition: DataTypes.hpp:262
GEOS_LOCALINDEX_TYPE localIndex
Local index type (for indexing objects within an MPI partition).
Definition: DataTypes.hpp:84
LvArray::SortedArray< T, localIndex, LvArray::ChaiBuffer > SortedArray
A sorted array of local indices.
Definition: DataTypes.hpp:266
LvArray::ArrayOfSets< T, INDEX_TYPE, LvArray::ChaiBuffer > ArrayOfSets
Array of variable-sized sets. See LvArray::ArrayOfSets for details.
Definition: DataTypes.hpp:289
LvArray::SortedArrayView< T, localIndex, LvArray::ChaiBuffer > SortedArrayView
A sorted array view of local indices.
Definition: DataTypes.hpp:270
int integer
Signed integer type.
Definition: DataTypes.hpp:81
Array< T, 1 > array1d
Alias for 1D array.
Definition: DataTypes.hpp:175
LvArray::ArrayView< T, NDIM, USD, localIndex, LvArray::ChaiBuffer > ArrayView
Multidimensional array view type. See LvArray:ArrayView for details.
Definition: DataTypes.hpp:147