ParticleManager.hpp

Go to the documentation of this file.

/*
 * ------------------------------------------------------------------------------------------------------------
 * SPDX-License-Identifier: LGPL-2.1-only
 *
 * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
 * Copyright (c) 2018-2024 TotalEnergies
 * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
 * Copyright (c) 2023-2024 Chevron
 * Copyright (c) 2019-     GEOS/GEOSX Contributors
 * All rights reserved
 *
 * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
 * ------------------------------------------------------------------------------------------------------------
 */
 
#ifndef GEOS_MESH_PARTICLEREGIONMANAGER_HPP
#define GEOS_MESH_PARTICLEREGIONMANAGER_HPP
 
#include "generators/ParticleBlock.hpp"
#include "generators/ParticleBlockManager.hpp"
#include "constitutive/ConstitutiveManager.hpp"
#include "ParticleRegion.hpp"
#include "ParticleSubRegion.hpp"
#include "mesh/ObjectManagerBase.hpp"
#include "dataRepository/ReferenceWrapper.hpp"
 
namespace geos
{
 
class MeshManager;
class DomainPartition;
 
class ParticleManager : public ObjectManagerBase
{
public:
 
  struct groupKeyStruct : public ObjectManagerBase::groupKeyStruct
  {
    static constexpr auto particleRegionsGroup() { return "particleRegionsGroup"; }
  };
 
  template< typename VIEWTYPE >
  using ParticleViewAccessor = array1d< array1d< VIEWTYPE > >;
 
  template< typename VIEWTYPE >
  using ParticleView = typename ParticleViewAccessor< VIEWTYPE >::NestedViewType;
 
  template< typename VIEWTYPE >
  using ParticleViewConst = typename ParticleViewAccessor< VIEWTYPE >::NestedViewTypeConst;
 
  template< typename VIEWTYPE >
  using ParticleReferenceAccessor = array1d< array1d< ReferenceWrapper< VIEWTYPE > > >;
 
  template< typename VIEWTYPE >
  using MaterialViewAccessor = array1d< array1d< array1d< VIEWTYPE > > >;
 
  template< typename CONSTITUTIVE_TYPE >
  using ConstitutiveRelationAccessor = array1d< array1d< array1d< CONSTITUTIVE_TYPE * > > >;
 
  static string catalogName()
  { return "ParticleManager"; }
 
  virtual string getCatalogName() const override final
  { return ParticleManager::catalogName(); }
 
  ParticleManager( string const & name, Group * const parent );
 
  virtual ~ParticleManager() override;
 
  template< typename T = ParticleSubRegionBase >
  localIndex getNumberOfParticles() const
  {
    localIndex numParticle = 0;
    this->forParticleSubRegions< T >( [&]( ParticleSubRegionBase const & particleSubRegion )
    {
      numParticle += particleSubRegion.size();
    } );
    return numParticle;
  }
 
//  void Initialize(  ){}
 
  void generateMesh( ParticleBlockManagerABC & particleBlockManager );
 
  virtual Group * createChild( string const & childKey, string const & childName ) override;
//  virtual void ReadXMLsub( xmlWrapper::xmlNode const & targetNode ) override;
 
  virtual void expandObjectCatalogs() override;
 
  virtual void setSchemaDeviations( xmlWrapper::xmlNode schemaRoot,
                                    xmlWrapper::xmlNode schemaParent,
                                    integer documentationType ) override;
 
  using Group::resize;
 
  void resize( integer_array const & numParticles,
               string_array const & regionNames,
               string_array const & particleTypes );
 
  virtual void setMaxGlobalIndex() override final;
 
  subGroupMap const & getRegions() const
  {
    return this->getGroup( groupKeyStruct::particleRegionsGroup() ).getSubGroups();
  }
 
  subGroupMap & getRegions()
  {
    return this->getGroup( groupKeyStruct::particleRegionsGroup() ).getSubGroups();
  }
 
  template< typename T=ParticleRegionBase, typename KEY_TYPE=void >
  T const & getRegion( KEY_TYPE const & key ) const
  {
    return this->getGroup( groupKeyStruct::particleRegionsGroup() ).getGroup< T >( key );
  }
 
  template< typename T=ParticleRegionBase, typename KEY_TYPE=void >
  T & getRegion( KEY_TYPE const & key )
  {
    return this->getGroup( groupKeyStruct::particleRegionsGroup() ).getGroup< T >( key );
  }
 
  template< typename T=ParticleRegionBase >
  bool hasRegion( string const & name ) const
  {
    return this->getGroup( groupKeyStruct::particleRegionsGroup() ).hasGroup< T >( name );
  }
 
  localIndex numRegions() const
  {
    return this->getRegions().size();
  }
 
  localIndex numParticleBlocks() const;
 
  template< typename REGIONTYPE = ParticleRegionBase, typename ... REGIONTYPES, typename LAMBDA >
  void forParticleRegions( LAMBDA && lambda )
  {
    this->getGroup( groupKeyStruct::particleRegionsGroup() ).forSubGroups< REGIONTYPE, REGIONTYPES... >( std::forward< LAMBDA >( lambda ) );
  }
 
  template< typename REGIONTYPE = ParticleRegionBase, typename ... REGIONTYPES, typename LAMBDA >
  void forParticleRegions( LAMBDA && lambda ) const
  {
    this->getGroup( groupKeyStruct::particleRegionsGroup() ).forSubGroups< REGIONTYPE, REGIONTYPES... >( std::forward< LAMBDA >( lambda ) );
  }
 
  template< typename REGIONTYPE = ParticleRegionBase, typename ... REGIONTYPES, typename LOOKUP_CONTAINER, typename LAMBDA >
  void forParticleRegions( LOOKUP_CONTAINER const & targetRegions, LAMBDA && lambda )
  {
    this->getGroup( groupKeyStruct::particleRegionsGroup() ).forSubGroups< REGIONTYPE, REGIONTYPES... >( targetRegions, std::forward< LAMBDA >( lambda ) );
  }
 
  template< typename REGIONTYPE = ParticleRegionBase, typename ... REGIONTYPES, typename LOOKUP_CONTAINER, typename LAMBDA >
  void forParticleRegions( LOOKUP_CONTAINER const & targetRegions, LAMBDA && lambda ) const
  {
    this->getGroup( groupKeyStruct::particleRegionsGroup() ).forSubGroups< REGIONTYPE, REGIONTYPES... >( targetRegions, std::forward< LAMBDA >( lambda ) );
  }
 
  template< typename LAMBDA >
  void forParticleRegionsComplete( LAMBDA lambda ) const
  {
    forParticleRegionsComplete< ParticleRegion >( std::forward< LAMBDA >( lambda ) );
  }
 
  template< typename LAMBDA >
  void forParticleRegionsComplete( LAMBDA lambda )
  {
    forParticleRegionsComplete< ParticleRegion >( std::forward< LAMBDA >( lambda ) );
  }
 
  template< typename REGIONTYPE, typename ... REGIONTYPES, typename LAMBDA >
  void forParticleRegionsComplete( LAMBDA lambda )
  {
    for( localIndex er=0; er<this->numRegions(); ++er )
    {
      ParticleRegionBase & particleRegion = this->getRegion( er );
 
      Group::applyLambdaToContainer< REGIONTYPE, REGIONTYPES... >( particleRegion, [&]( auto & castedRegion )
      {
        lambda( er, castedRegion );
      } );
    }
  }
 
  template< typename REGIONTYPE, typename ... REGIONTYPES, typename LAMBDA >
  void forParticleRegionsComplete( LAMBDA lambda ) const
  {
    for( localIndex er=0; er<this->numRegions(); ++er )
    {
      ParticleRegionBase const & particleRegion = this->getRegion( er );
 
      Group::applyLambdaToContainer< REGIONTYPE, REGIONTYPES... >( particleRegion, [&]( auto const & castedRegion )
      {
        lambda( er, castedRegion );
      } );
    }
  }
 
  template< typename LOOKUP_CONTAINER, typename LAMBDA >
  void forParticleRegionsComplete( LOOKUP_CONTAINER const & targetRegions, LAMBDA lambda ) const
  {
    forParticleRegionsComplete< ParticleRegion >( targetRegions, std::forward< LAMBDA >( lambda ) );
  }
 
  template< typename LOOKUP_CONTAINER, typename LAMBDA >
  void forParticleRegionsComplete( LOOKUP_CONTAINER const & targetRegions, LAMBDA lambda )
  {
    forParticleRegionsComplete< ParticleRegion >( targetRegions, std::forward< LAMBDA >( lambda ) );
  }
 
  template< typename REGIONTYPE, typename ... REGIONTYPES, typename LOOKUP_CONTAINER, typename LAMBDA >
  void forParticleRegionsComplete( LOOKUP_CONTAINER const & targetRegions, LAMBDA lambda )
  {
    forParticleRegions< REGIONTYPE, REGIONTYPES... >( targetRegions, [&] ( localIndex const targetIndex,
                                                                           auto & particleRegion )
    {
      lambda( targetIndex, particleRegion.getIndexInParent(), particleRegion );
    } );
  }
 
  template< typename REGIONTYPE, typename ... REGIONTYPES, typename LOOKUP_CONTAINER, typename LAMBDA >
  void forParticleRegionsComplete( LOOKUP_CONTAINER const & targetRegions, LAMBDA lambda ) const
  {
    forParticleRegions< REGIONTYPE, REGIONTYPES... >( targetRegions, [&] ( localIndex const targetIndex,
                                                                           auto const & particleRegion )
    {
      lambda( targetIndex, particleRegion.getIndexInParent(), particleRegion );
    } );
  }
 
  template< typename LAMBDA >
  void forParticleSubRegions( LAMBDA && lambda )
  {
    forParticleSubRegions< ParticleSubRegion >( std::forward< LAMBDA >( lambda ) );
  }
 
  template< typename LAMBDA >
  void forParticleSubRegions( LAMBDA && lambda ) const
  {
    forParticleSubRegions< ParticleSubRegion >( std::forward< LAMBDA >( lambda ) );
  }
 
  template< typename LOOKUP_CONTAINER, typename LAMBDA >
  void forParticleSubRegions( LOOKUP_CONTAINER const & targetRegions, LAMBDA && lambda )
  {
    forParticleSubRegions< ParticleSubRegion >( targetRegions, std::forward< LAMBDA >( lambda ) );
  }
 
  template< typename LOOKUP_CONTAINER, typename LAMBDA >
  void forParticleSubRegions( LOOKUP_CONTAINER const & targetRegions, LAMBDA && lambda ) const
  {
    forParticleSubRegions< ParticleSubRegion >( targetRegions, std::forward< LAMBDA >( lambda ) );
  }
 
  template< typename SUBREGIONTYPE, typename ... SUBREGIONTYPES, typename LAMBDA >
  void forParticleSubRegions( LAMBDA && lambda )
  {
    forParticleSubRegionsComplete< SUBREGIONTYPE, SUBREGIONTYPES... >(
      [lambda = std::forward< LAMBDA >( lambda )]( localIndex const,
                                                   localIndex const,
                                                   ParticleRegionBase &,
                                                   auto & subRegion )
    {
      lambda( subRegion );
    }
      );
  }
 
  template< typename SUBREGIONTYPE, typename ... SUBREGIONTYPES, typename LAMBDA >
  void forParticleSubRegions( LAMBDA && lambda ) const
  {
    forParticleSubRegionsComplete< SUBREGIONTYPE, SUBREGIONTYPES... >(
      [lambda = std::forward< LAMBDA >( lambda )]( localIndex const,
                                                   localIndex const,
                                                   ParticleRegionBase const &,
                                                   auto const & subRegion )
    {
      lambda( subRegion );
    } );
  }
 
  template< typename SUBREGIONTYPE, typename ... SUBREGIONTYPES, typename LOOKUP_CONTAINER, typename LAMBDA >
  void forParticleSubRegions( LOOKUP_CONTAINER const & targetRegions, LAMBDA && lambda )
  {
    forParticleSubRegionsComplete< SUBREGIONTYPE, SUBREGIONTYPES... >( targetRegions,
                                                                       [lambda = std::forward< LAMBDA >( lambda )]( localIndex const targetIndex,
                                                                                                                    localIndex const,
                                                                                                                    localIndex const,
                                                                                                                    ParticleRegionBase &,
                                                                                                                    auto & subRegion )
    {
      lambda( targetIndex, subRegion );
    } );
  }
 
  template< typename SUBREGIONTYPE, typename ... SUBREGIONTYPES, typename LOOKUP_CONTAINER, typename LAMBDA >
  void forParticleSubRegions( LOOKUP_CONTAINER const & targetRegions, LAMBDA && lambda ) const
  {
    forParticleSubRegionsComplete< SUBREGIONTYPE, SUBREGIONTYPES... >( targetRegions,
                                                                       [lambda = std::forward< LAMBDA >( lambda )]( localIndex const targetIndex,
                                                                                                                    localIndex const,
                                                                                                                    localIndex const,
                                                                                                                    ParticleRegionBase const &,
                                                                                                                    auto const & subRegion )
    {
      lambda( targetIndex, subRegion );
    } );
  }
 
  template< typename LAMBDA >
  void forParticleSubRegionsComplete( LAMBDA && lambda ) const
  {
    forParticleSubRegionsComplete< ParticleSubRegion >( std::forward< LAMBDA >( lambda ) );
  }
 
  template< typename LAMBDA >
  void forParticleSubRegionsComplete( LAMBDA && lambda )
  {
    forParticleSubRegionsComplete< ParticleSubRegion >( std::forward< LAMBDA >( lambda ) );
  }
 
  template< typename LOOKUP_CONTAINER, typename LAMBDA >
  void forParticleSubRegionsComplete( LOOKUP_CONTAINER const & targetRegions, LAMBDA && lambda )
  {
    forParticleSubRegionsComplete< ParticleSubRegion >( targetRegions, std::forward< LAMBDA >( lambda ) );
  }
 
  template< typename LOOKUP_CONTAINER, typename LAMBDA >
  void forParticleSubRegionsComplete( LOOKUP_CONTAINER const & targetRegions, LAMBDA && lambda ) const
  {
    forParticleSubRegionsComplete< ParticleSubRegion >( targetRegions, std::forward< LAMBDA >( lambda ) );
  }
 
  template< typename SUBREGIONTYPE, typename ... SUBREGIONTYPES, typename LAMBDA >
  void forParticleSubRegionsComplete( LAMBDA && lambda )
  {
    for( localIndex er=0; er<this->numRegions(); ++er )
    {
      ParticleRegionBase & particleRegion = this->getRegion( er );
 
      for( localIndex esr=0; esr<particleRegion.numSubRegions(); ++esr )
      {
        ParticleSubRegionBase & subRegion = particleRegion.getSubRegion( esr );
 
        Group::applyLambdaToContainer< SUBREGIONTYPE, SUBREGIONTYPES... >( subRegion, [&]( auto & castedSubRegion )
        {
          lambda( er, esr, particleRegion, castedSubRegion );
        } );
      }
    }
  }
 
  template< typename SUBREGIONTYPE, typename ... SUBREGIONTYPES, typename LAMBDA >
  void forParticleSubRegionsComplete( LAMBDA && lambda ) const
  {
    for( localIndex er=0; er<this->numRegions(); ++er )
    {
      ParticleRegionBase const & particleRegion = this->getRegion( er );
 
      for( localIndex esr=0; esr<particleRegion.numSubRegions(); ++esr )
      {
        ParticleSubRegionBase const & subRegion = particleRegion.getSubRegion( esr );
 
        Group::applyLambdaToContainer< SUBREGIONTYPE, SUBREGIONTYPES... >( subRegion, [&]( auto const & castedSubRegion )
        {
          lambda( er, esr, particleRegion, castedSubRegion );
        } );
      }
    }
  }
 
  template< typename SUBREGIONTYPE, typename ... SUBREGIONTYPES, typename LOOKUP_CONTAINER, typename LAMBDA >
  void forParticleSubRegionsComplete( LOOKUP_CONTAINER const & targetRegions, LAMBDA && lambda )
  {
    forParticleRegions( targetRegions, [&] ( localIndex const targetIndex, ParticleRegionBase & particleRegion )
    {
      localIndex const er = particleRegion.getIndexInParent();
 
      if( er>-1 )
      {
        for( localIndex esr=0; esr<particleRegion.numSubRegions(); ++esr )
        {
          ParticleSubRegionBase & subRegion = particleRegion.getSubRegion( esr );
 
          Group::applyLambdaToContainer< SUBREGIONTYPE, SUBREGIONTYPES... >( subRegion, [&]( auto & castedSubRegion )
          {
            lambda( targetIndex, er, esr, particleRegion, castedSubRegion );
          } );
        }
      }
    } );
  }
 
  template< typename SUBREGIONTYPE, typename ... SUBREGIONTYPES, typename LOOKUP_CONTAINER, typename LAMBDA >
  void forParticleSubRegionsComplete( LOOKUP_CONTAINER const & targetRegions, LAMBDA && lambda ) const
  {
    forParticleRegions( targetRegions, [&] ( localIndex const targetIndex, ParticleRegionBase const & particleRegion )
    {
      localIndex const er = particleRegion.getIndexInParent();
 
      if( er>-1 )
      {
        for( localIndex esr=0; esr<particleRegion.numSubRegions(); ++esr )
        {
          ParticleSubRegionBase const & subRegion = particleRegion.getSubRegion( esr );
 
          Group::applyLambdaToContainer< SUBREGIONTYPE, SUBREGIONTYPES... >( subRegion, [&]( auto const & castedSubRegion )
          {
            lambda( targetIndex, er, esr, particleRegion, castedSubRegion );
          } );
        }
      }
    } );
  }
 
 
  template< typename FIELD_TRAIT >
  ParticleViewAccessor< traits::ViewTypeConst< typename FIELD_TRAIT::type > >
  constructFieldAccessor( string const & neighborName = string() ) const;
 
  template< typename VIEWTYPE, typename LHS=VIEWTYPE >
  ParticleViewAccessor< LHS >
  constructViewAccessor( string const & name, string const & neighborName = string() ) const;
 
  template< typename VIEWTYPE, typename LHS=VIEWTYPE >
  ParticleViewAccessor< LHS >
  constructViewAccessor( string const & name, string const & neighborName = string() );
 
  template< typename T, int NDIM, typename PERM = defaultLayout< NDIM > >
  ParticleViewAccessor< ArrayView< T const, NDIM, getUSD< PERM > > >
  constructArrayViewAccessor( string const & name, string const & neighborName = string() ) const;
 
  template< typename VIEWTYPE >
  ParticleViewAccessor< ReferenceWrapper< VIEWTYPE > >
  constructReferenceAccessor( string const & viewName, string const & neighborName = string() ) const;
 
  template< typename VIEWTYPE >
  ParticleViewAccessor< ReferenceWrapper< VIEWTYPE > >
  constructReferenceAccessor( string const & viewName, string const & neighborName = string() );
 
  template< typename VIEWTYPE, typename LHS=VIEWTYPE >
  MaterialViewAccessor< LHS >
  constructFullMaterialViewAccessor( string const & viewName,
                                     constitutive::ConstitutiveManager const & cm ) const;
 
  template< typename VIEWTYPE, typename LHS=VIEWTYPE >
  MaterialViewAccessor< LHS >
  constructFullMaterialViewAccessor( string const & viewName,
                                     constitutive::ConstitutiveManager const & cm );
 
  template< typename FIELD_TRAIT >
  ParticleViewAccessor< traits::ViewTypeConst< typename FIELD_TRAIT::type > >
  constructMaterialFieldAccessor( string_array const & regionNames,
                                  string_array const & materialNames,
                                  bool const allowMissingViews = false ) const;
 
  template< typename MATERIAL_TYPE, typename FIELD_TRAIT >
  ParticleViewAccessor< traits::ViewTypeConst< typename FIELD_TRAIT::type > >
  constructMaterialFieldAccessor( bool const allowMissingViews = false ) const;
 
  template< typename VIEWTYPE, typename LHS=VIEWTYPE >
  ParticleViewAccessor< LHS >
  constructMaterialViewAccessor( string const & viewName,
                                 string_array const & regionNames,
                                 string const & materialKeyName,
                                 bool const allowMissingViews = false ) const;
 
  template< typename VIEWTYPE, typename LHS=VIEWTYPE >
  ParticleViewAccessor< LHS >
  constructMaterialViewAccessor( string const & viewName,
                                 string_array const & regionNames,
                                 string const & materialKeyName,
                                 bool const allowMissingViews = false );
 
  template< typename T, int NDIM, typename PERM = defaultLayout< NDIM > >
  ParticleViewAccessor< ArrayView< T const, NDIM, getUSD< PERM > > >
  constructMaterialArrayViewAccessor( string const & viewName,
                                      string_array const & regionNames,
                                      string const & materialKeyName,
                                      bool const allowMissingViews = false ) const;
 
  template< typename MATERIALTYPE, typename VIEWTYPE, typename LHS=VIEWTYPE >
  ParticleViewAccessor< LHS >
  constructMaterialViewAccessor( string const & viewName ) const;
 
  template< typename MATERIALTYPE, typename T, int NDIM, typename PERM = defaultLayout< NDIM > >
  ParticleViewAccessor< ArrayView< T const, NDIM, getUSD< PERM > > >
  constructMaterialArrayViewAccessor( string const & viewName ) const;
 
  template< typename CONSTITUTIVE_TYPE >
  ConstitutiveRelationAccessor< CONSTITUTIVE_TYPE >
  constructFullConstitutiveAccessor( constitutive::ConstitutiveManager const & cm ) const;
 
 
  template< typename CONSTITUTIVE_TYPE >
  ConstitutiveRelationAccessor< CONSTITUTIVE_TYPE >
  constructFullConstitutiveAccessor( constitutive::ConstitutiveManager const & cm );
 
  using Group::packSize;
  using Group::pack;
  using ObjectManagerBase::packGlobalMapsSize;
  using ObjectManagerBase::packGlobalMaps;
  using ObjectManagerBase::unpackGlobalMaps;
 
  int PackSize( string_array const & wrapperNames,
                ParticleViewAccessor< arrayView1d< localIndex > > const & packList ) const;
 
  int Pack( buffer_unit_type * & buffer,
            string_array const & wrapperNames,
            ParticleViewAccessor< arrayView1d< localIndex > > const & packList ) const;
 
  using ObjectManagerBase::unpack;
 
  int Unpack( buffer_unit_type const * & buffer,
              ParticleViewAccessor< arrayView1d< localIndex > > & packList );
 
  int Unpack( buffer_unit_type const * & buffer,
              ParticleReferenceAccessor< array1d< localIndex > > & packList );
 
  int PackGlobalMapsSize( ParticleViewAccessor< arrayView1d< localIndex > > const & packList ) const;
 
  int PackGlobalMaps( buffer_unit_type * & buffer,
                      ParticleViewAccessor< arrayView1d< localIndex > > const & packList ) const;
 
  int UnpackGlobalMaps( buffer_unit_type const * & buffer,
                        ParticleViewAccessor< ReferenceWrapper< localIndex_array > > & packList );
 
  void updateMaps();
 
 
private:
 
  template< bool DOPACK >
  int PackPrivate( buffer_unit_type * & buffer,
                   string_array const & wrapperNames,
                   ParticleViewAccessor< arrayView1d< localIndex > > const & viewAccessor ) const;
 
  template< bool DOPACK >
  int PackGlobalMapsPrivate( buffer_unit_type * & buffer,
                             ParticleViewAccessor< arrayView1d< localIndex > > const & viewAccessor ) const;
 
  template< typename T >
  int unpackPrivate( buffer_unit_type const * & buffer,
                     T & packList );
 
  ParticleManager( const ParticleManager & );
 
  ParticleManager & operator=( const ParticleManager & );
 
};
 
 
template< typename VIEWTYPE, typename LHS >
ParticleManager::ParticleViewAccessor< LHS >
ParticleManager::constructViewAccessor( string const & viewName, string const & neighborName ) const
{
  ParticleViewAccessor< LHS > viewAccessor;
  viewAccessor.resize( numRegions() );
  for( typename dataRepository::indexType kReg=0; kReg<numRegions(); ++kReg )
  {
    ParticleRegionBase const & particleRegion = getRegion( kReg );
    viewAccessor[kReg].resize( particleRegion.numSubRegions() );
 
    for( typename dataRepository::indexType kSubReg = 0; kSubReg < particleRegion.numSubRegions(); ++kSubReg )
    {
      Group const * group = &particleRegion.getSubRegion( kSubReg );
 
      if( !neighborName.empty() )
      {
        group = &group->getGroup( ObjectManagerBase::groupKeyStruct::neighborDataString() ).getGroup( neighborName );
      }
 
      dataRepository::Wrapper< VIEWTYPE > const * const wrapper = group->getWrapperPointer< VIEWTYPE >( viewName );
      if( wrapper )
      {
        viewAccessor[kReg][kSubReg] = wrapper->reference();
      }
    }
  }
  return viewAccessor;
}
 
 
template< typename VIEWTYPE, typename LHS >
ParticleManager::ParticleViewAccessor< LHS >
ParticleManager::constructViewAccessor( string const & viewName, string const & neighborName )
{
  ParticleViewAccessor< LHS > viewAccessor;
  viewAccessor.resize( numRegions() );
  for( typename dataRepository::indexType kReg=0; kReg<numRegions(); ++kReg )
  {
    ParticleRegionBase & particleRegion = getRegion( kReg );
    viewAccessor[kReg].resize( particleRegion.numSubRegions() );
 
    for( typename dataRepository::indexType kSubReg = 0; kSubReg < particleRegion.numSubRegions(); ++kSubReg )
    {
      Group * group = &particleRegion.getSubRegion( kSubReg );
 
      if( !neighborName.empty() )
      {
        group = &group->getGroup( ObjectManagerBase::groupKeyStruct::neighborDataString() ).getGroup( neighborName );
      }
 
      dataRepository::Wrapper< VIEWTYPE > * const wrapper = group->getWrapperPointer< VIEWTYPE >( viewName );
      if( wrapper )
      {
        viewAccessor[kReg][kSubReg] = wrapper->reference();
      }
    }
  }
  return viewAccessor;
}
 
template< typename FIELD_TRAIT >
ParticleManager::ParticleViewAccessor< traits::ViewTypeConst< typename FIELD_TRAIT::type > >
ParticleManager::constructFieldAccessor( string const & neighborName ) const
{
  return constructViewAccessor< typename FIELD_TRAIT::type,
                                traits::ViewTypeConst< typename FIELD_TRAIT::type > >( FIELD_TRAIT::key(), neighborName );
}
 
template< typename T, int NDIM, typename PERM >
ParticleManager::ParticleViewAccessor< ArrayView< T const, NDIM, getUSD< PERM > > >
ParticleManager::constructArrayViewAccessor( string const & name, string const & neighborName ) const
{
  return constructViewAccessor< Array< T, NDIM, PERM >,
                                ArrayView< T const, NDIM, getUSD< PERM > >
                                >( name, neighborName );
}
 
template< typename VIEWTYPE >
ParticleManager::ParticleViewAccessor< ReferenceWrapper< VIEWTYPE > >
ParticleManager::constructReferenceAccessor( string const & viewName, string const & neighborName ) const
{
  ParticleViewAccessor< ReferenceWrapper< VIEWTYPE > > viewAccessor;
  viewAccessor.resize( numRegions() );
  for( typename dataRepository::indexType kReg=0; kReg<numRegions(); ++kReg )
  {
    ParticleRegionBase const & particleRegion = getRegion( kReg );
    viewAccessor[kReg].resize( particleRegion.numSubRegions() );
 
    for( typename dataRepository::indexType kSubReg=0; kSubReg<particleRegion.numSubRegions(); ++kSubReg )
    {
      Group const * group = &particleRegion.getSubRegion( kSubReg );
 
      if( !neighborName.empty() )
      {
        group = &group->getGroup( ObjectManagerBase::groupKeyStruct::neighborDataString() ).getGroup( neighborName );
      }
 
      if( group->hasWrapper( viewName ) )
      {
        viewAccessor[kReg][kSubReg].set( group->getReference< VIEWTYPE >( viewName ) );
      }
    }
  }
  return viewAccessor;
}
 
template< typename VIEWTYPE >
ParticleManager::ParticleViewAccessor< ReferenceWrapper< VIEWTYPE > >
ParticleManager::constructReferenceAccessor( string const & viewName, string const & neighborName )
{
  ParticleViewAccessor< ReferenceWrapper< VIEWTYPE > > viewAccessor;
  viewAccessor.resize( numRegions() );
  for( typename dataRepository::indexType kReg=0; kReg<numRegions(); ++kReg )
  {
    ParticleRegionBase & particleRegion = getRegion( kReg );
    viewAccessor[kReg].resize( particleRegion.numSubRegions() );
 
    for( typename dataRepository::indexType kSubReg=0; kSubReg<particleRegion.numSubRegions(); ++kSubReg )
    {
      Group * group = &particleRegion.getSubRegion( kSubReg );
 
      if( !neighborName.empty() )
      {
        group = &group->getGroup( ObjectManagerBase::groupKeyStruct::neighborDataString() ).getGroup( neighborName );
      }
 
      if( group->hasWrapper( viewName ) )
      {
        viewAccessor[kReg][kSubReg].set( group->getReference< VIEWTYPE >( viewName ) );
      }
    }
  }
  return viewAccessor;
}
 
template< typename VIEWTYPE, typename LHS >
ParticleManager::MaterialViewAccessor< LHS >
ParticleManager::constructFullMaterialViewAccessor( string const & viewName,
                                                    constitutive::ConstitutiveManager const & cm ) const
{
  MaterialViewAccessor< LHS > accessor;
  accessor.resize( numRegions() );
  for( localIndex kReg=0; kReg<numRegions(); ++kReg )
  {
    ParticleRegionBase const & particleRegion = getRegion( kReg );
    accessor[kReg].resize( particleRegion.numSubRegions() );
 
    for( localIndex kSubReg=0; kSubReg<particleRegion.numSubRegions(); ++kSubReg )
    {
      ParticleSubRegionBase const & subRegion = particleRegion.getSubRegion( kSubReg );
      dataRepository::Group const & constitutiveGroup = subRegion.getConstitutiveModels();
 
      accessor[kReg][kSubReg].resize( cm.numSubGroups() );
 
      for( localIndex matIndex=0; matIndex<cm.numSubGroups(); ++matIndex )
      {
        string const & constitutiveName = cm.getGroup( matIndex ).getName();
        dataRepository::Group const * const constitutiveRelation = constitutiveGroup.getGroupPointer( constitutiveName );
        if( constitutiveRelation != nullptr )
        {
          dataRepository::Wrapper< VIEWTYPE > const * const wrapper = constitutiveRelation->getWrapperPointer< VIEWTYPE >( viewName );
          if( wrapper )
          {
            accessor[kReg][kSubReg][matIndex] = wrapper->reference();
          }
        }
      }
    }
  }
  return accessor;
}
 
template< typename VIEWTYPE, typename LHS >
ParticleManager::MaterialViewAccessor< LHS >
ParticleManager::constructFullMaterialViewAccessor( string const & viewName,
                                                    constitutive::ConstitutiveManager const & cm )
{
  MaterialViewAccessor< LHS > accessor;
  accessor.resize( numRegions() );
  for( localIndex kReg=0; kReg<numRegions(); ++kReg )
  {
    ParticleRegionBase & particleRegion = getRegion( kReg );
    accessor[kReg].resize( particleRegion.numSubRegions() );
 
    for( localIndex kSubReg=0; kSubReg<particleRegion.numSubRegions(); ++kSubReg )
    {
      ParticleSubRegionBase & subRegion = particleRegion.getSubRegion( kSubReg );
      dataRepository::Group & constitutiveGroup = subRegion.getConstitutiveModels();
 
      accessor[kReg][kSubReg].resize( cm.numSubGroups() );
 
      for( localIndex matIndex=0; matIndex<cm.numSubGroups(); ++matIndex )
      {
        string const & constitutiveName = cm.getGroup( matIndex ).getName();
        dataRepository::Group * const constitutiveRelation = constitutiveGroup.getGroupPointer( constitutiveName );
        if( constitutiveRelation != nullptr )
        {
          dataRepository::Wrapper< VIEWTYPE > * const wrapper = constitutiveRelation->getWrapperPointer< VIEWTYPE >( viewName );
          if( wrapper )
          {
            accessor[kReg][kSubReg][matIndex] = wrapper->reference();
          }
        }
      }
    }
  }
  return accessor;
}
 
template< typename VIEWTYPE, typename LHS >
ParticleManager::ParticleViewAccessor< LHS >
ParticleManager::constructMaterialViewAccessor( string const & viewName,
                                                string_array const & regionNames,
                                                string const & materialKeyName,
                                                bool const allowMissingViews ) const
{
  ParticleViewAccessor< LHS > accessor;
 
  // Resize the accessor to all regions and subregions
  accessor.resize( numRegions() );
  for( localIndex kReg = 0; kReg < numRegions(); ++kReg )
  {
    accessor[kReg].resize( getRegion( kReg ).numSubRegions() );
  }
 
  subGroupMap const & regionMap = getRegions();
 
  // Loop only over regions named and populate according to given material names
  for( size_t k = 0; k < regionNames.size(); ++k )
  {
    localIndex const er = regionMap.getIndex( regionNames[k] );
    if( er >=0 )
    {
      GEOS_ERROR_IF_EQ_MSG( er, subGroupMap::KeyIndex::invalid_index, "Region not found: " << regionNames[k] );
      ParticleRegionBase const & region = getRegion( er );
 
      region.forParticleSubRegionsIndex( [&]( localIndex const esr,
                                              ParticleSubRegionBase const & subRegion )
      {
        string const & materialName = subRegion.getReference< string >( materialKeyName );
        dataRepository::Group const & constitutiveRelation = subRegion.getConstitutiveModel( materialName );
 
        dataRepository::Wrapper< VIEWTYPE > const * const wrapper = constitutiveRelation.getWrapperPointer< VIEWTYPE >( viewName );
        if( wrapper )
        {
          accessor[er][esr] = wrapper->reference();
        }
        else
        {
          GEOS_ERROR_IF( !allowMissingViews, "Material " << materialKeyName[k] << " does not contain " << viewName );
        }
      } );
    }
  }
  return accessor;
}
 
template< typename VIEWTYPE, typename LHS >
ParticleManager::ParticleViewAccessor< LHS >
ParticleManager::constructMaterialViewAccessor( string const & viewName,
                                                string_array const & regionNames,
                                                string const & materialKeyName,
                                                bool const allowMissingViews )
{
  ParticleViewAccessor< LHS > accessor;
 
  // Resize the accessor to all regions and subregions
  accessor.resize( numRegions() );
  for( localIndex kReg = 0; kReg < numRegions(); ++kReg )
  {
    accessor[kReg].resize( getRegion( kReg ).numSubRegions() );
  }
 
  subGroupMap const & regionMap = getRegions();
 
  // Loop only over regions named and populate according to given material names
  for( size_t k = 0; k < regionNames.size(); ++k )
  {
    localIndex const er = regionMap.getIndex( regionNames[k] );
    if( er >=0 )
    {
      GEOS_ERROR_IF_EQ_MSG( er, subGroupMap::KeyIndex::invalid_index, "Region not found: " << regionNames[k] );
      ParticleRegionBase & region = getRegion( er );
 
      region.forParticleSubRegionsIndex( [&]( localIndex const esr, ParticleSubRegionBase & subRegion )
      {
        string const & materialName = subRegion.getReference< string >( materialKeyName );
        dataRepository::Group const & constitutiveRelation = subRegion.getConstitutiveModel( materialName );
 
        dataRepository::Wrapper< VIEWTYPE > * const wrapper = constitutiveRelation.getWrapperPointer< VIEWTYPE >( viewName );
        if( wrapper )
        {
          accessor[er][esr] = wrapper->reference();
        }
        else
        {
          GEOS_ERROR_IF( !allowMissingViews, "Material " << materialName << " does not contain " << viewName );
        }
      } );
    }
  }
  return accessor;
}
 
template< typename FIELD_TRAIT >
ParticleManager::ParticleViewAccessor< traits::ViewTypeConst< typename FIELD_TRAIT::type > >
ParticleManager::constructMaterialFieldAccessor( string_array const & regionNames,
                                                 string_array const & materialNames,
                                                 bool const allowMissingViews ) const
{
  return constructMaterialViewAccessor< typename FIELD_TRAIT::type,
                                        traits::ViewTypeConst< typename FIELD_TRAIT::type > >( FIELD_TRAIT::key(),
                                                                                               regionNames,
                                                                                               materialNames,
                                                                                               allowMissingViews );
}
 
template< typename MATERIAL_TYPE, typename FIELD_TRAIT >
ParticleManager::ParticleViewAccessor< traits::ViewTypeConst< typename FIELD_TRAIT::type > >
ParticleManager::constructMaterialFieldAccessor( bool const allowMissingViews ) const
{
  GEOS_UNUSED_VAR( allowMissingViews );
  return constructMaterialViewAccessor< MATERIAL_TYPE, typename FIELD_TRAIT::type,
                                        traits::ViewTypeConst< typename FIELD_TRAIT::type > >( FIELD_TRAIT::key() );
}
 
template< typename T, int NDIM, typename PERM >
ParticleManager::ParticleViewAccessor< ArrayView< T const, NDIM, getUSD< PERM > > >
ParticleManager::constructMaterialArrayViewAccessor( string const & viewName,
                                                     string_array const & regionNames,
                                                     string const & materialKeyName,
                                                     bool const allowMissingViews ) const
{
  return constructMaterialViewAccessor< Array< T, NDIM, PERM >, ArrayView< T const, NDIM, getUSD< PERM > > >( viewName,
                                                                                                              regionNames,
                                                                                                              materialKeyName,
                                                                                                              allowMissingViews );
}
 
template< typename MATERIALTYPE, typename VIEWTYPE, typename LHS >
ParticleManager::ParticleViewAccessor< LHS >
ParticleManager::constructMaterialViewAccessor( string const & viewName ) const
{
  ParticleViewAccessor< LHS > accessor( numRegions() );
 
  // Resize the accessor to all regions and subregions
  for( localIndex er = 0; er < numRegions(); ++er )
  {
    accessor[er].resize( getRegion( er ).numSubRegions() );
  }
 
  // Loop only over regions named and populate according to given material names
  for( localIndex er = 0; er < numRegions(); ++er )
  {
    ParticleRegionBase const & region = getRegion( er );
 
    region.forParticleSubRegionsIndex( [&]( localIndex const esr,
                                            ParticleSubRegionBase const & subRegion )
    {
      dataRepository::Group const & constitutiveGroup = subRegion.getConstitutiveModels();
 
      string materialName;
      constitutiveGroup.forSubGroups< MATERIALTYPE >( [&]( MATERIALTYPE const & constitutiveRelation )
      {
        materialName = constitutiveRelation.getName();
        if( constitutiveRelation.template hasWrapper<>( viewName ) )  //NOTE (matteo): I have added this check to allow for the view to be
                                                                      // missing. I am not sure this is the default behaviour we want
                                                                      // though.
        {
          accessor[er][esr] = constitutiveRelation.template getReference< VIEWTYPE >( viewName );
        }
      } );
    } );
  }
  return accessor;
}
 
template< typename MATERIALTYPE, typename T, int NDIM, typename PERM >
ParticleManager::ParticleViewAccessor< ArrayView< T const, NDIM, getUSD< PERM > > >
ParticleManager::constructMaterialArrayViewAccessor( string const & viewName ) const
{
  return constructMaterialViewAccessor< MATERIALTYPE, Array< T, NDIM, PERM >, ArrayView< T const, NDIM, getUSD< PERM > > >( viewName );
}
 
template< typename CONSTITUTIVE_TYPE >
ParticleManager::ConstitutiveRelationAccessor< CONSTITUTIVE_TYPE >
ParticleManager::constructFullConstitutiveAccessor( constitutive::ConstitutiveManager const & cm ) const
{
  ConstitutiveRelationAccessor< CONSTITUTIVE_TYPE > accessor;
  accessor.resize( numRegions() );
  for( localIndex kReg=0; kReg<numRegions(); ++kReg )
  {
    ParticleRegionBase const & particleRegion = getRegion( kReg );
    accessor[kReg].resize( particleRegion.numSubRegions() );
 
    for( localIndex kSubReg=0; kSubReg<particleRegion.numSubRegions(); ++kSubReg )
    {
      ParticleSubRegionBase const & subRegion = particleRegion.getSubRegion( kSubReg );
      dataRepository::Group const & constitutiveGroup = subRegion.getConstitutiveModels();
      accessor[kReg][kSubReg].resize( cm.numSubGroups() );
 
      for( localIndex matIndex=0; matIndex<cm.numSubGroups(); ++matIndex )
      {
        string const & constitutiveName = cm.getGroup( matIndex ).getName();
 
        CONSTITUTIVE_TYPE * const
        constitutiveRelation = constitutiveGroup.getGroupPointer< CONSTITUTIVE_TYPE >( constitutiveName );
        if( constitutiveRelation != nullptr )
        {
          accessor[kReg][kSubReg][matIndex] = constitutiveRelation;
        }
      }
    }
  }
  return accessor;
}
 
template< typename CONSTITUTIVE_TYPE >
ParticleManager::ConstitutiveRelationAccessor< CONSTITUTIVE_TYPE >
ParticleManager::constructFullConstitutiveAccessor( constitutive::ConstitutiveManager const & cm )
{
  ConstitutiveRelationAccessor< CONSTITUTIVE_TYPE > accessor;
  accessor.resize( numRegions() );
  for( localIndex kReg=0; kReg<numRegions(); ++kReg )
  {
    ParticleRegionBase & particleRegion = getRegion( kReg );
    accessor[kReg].resize( particleRegion.numSubRegions() );
 
    for( localIndex kSubReg=0; kSubReg<particleRegion.numSubRegions(); ++kSubReg )
    {
      ParticleSubRegionBase & subRegion = particleRegion.getSubRegion( kSubReg );
      dataRepository::Group & constitutiveGroup = subRegion.getConstitutiveModels();
      accessor[kReg][kSubReg].resize( cm.numSubGroups() );
 
      for( localIndex matIndex=0; matIndex<cm.numSubGroups(); ++matIndex )
      {
        string const & constitutiveName = cm.getGroup( matIndex ).getName();
 
        CONSTITUTIVE_TYPE * const
        constitutiveRelation = constitutiveGroup.getGroupPointer< CONSTITUTIVE_TYPE >( constitutiveName );
        if( constitutiveRelation != nullptr )
        {
          accessor[kReg][kSubReg][matIndex] = constitutiveRelation;
        }
      }
    }
  }
  return accessor;
}
 
}
#endif /* GEOS_MESH_PARTICLEREGIONMANAGER_HPP */