CellBlock.hpp

Go to the documentation of this file.

/*
 * ------------------------------------------------------------------------------------------------------------
 * SPDX-License-Identifier: LGPL-2.1-only
 *
 * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC
 * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University
 * Copyright (c) 2018-2020 Total, S.A
 * Copyright (c) 2019-     GEOSX Contributors
 * All rights reserved
 *
 * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
 * ------------------------------------------------------------------------------------------------------------
 */

#ifndef GEOSX_MESH_CELLBLOCK_HPP_
#define GEOSX_MESH_CELLBLOCK_HPP_

#include "ElementSubRegionBase.hpp"
#include "FaceManager.hpp"
#include "meshUtilities/ComputationalGeometry.hpp"
#include "rajaInterface/GEOS_RAJA_Interface.hpp"


class StableTimeStep;

namespace geosx
{

class CellBlock : public ElementSubRegionBase
{
public:

  using NodeMapType = InterObjectRelation< array2d< localIndex, cells::NODE_MAP_PERMUTATION > >;
  using EdgeMapType = FixedOneToManyRelation;
  using FaceMapType = FixedOneToManyRelation;


  static const string CatalogName()
  { return "CellBlock"; }

  virtual const string getCatalogName() const override final
  { return CellBlock::CatalogName(); }



  CellBlock() = delete;

  CellBlock( string const & name, Group * const parent );

  CellBlock( const CellBlock & init ) = delete;

  virtual ~CellBlock() override;



  virtual void CalculateElementGeometricQuantities( NodeManager const & nodeManager,
                                                    FaceManager const & facemanager ) override;

  virtual void setupRelatedObjectsInRelations( MeshLevel const * const mesh ) override;

  void calculateElementCenters( arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & X ) const
  {
    arrayView2d< real64 > const & elementCenters = m_elementCenter;
    localIndex nNodes = numNodesPerElement();

    forAll< parallelHostPolicy >( size(), [=]( localIndex const k )
    {
      LvArray::tensorOps::copy< 3 >( elementCenters[ k ], X[ m_toNodesRelation( k, 0 ) ] );
      for( localIndex a = 1; a < nNodes; ++a )
      {
        LvArray::tensorOps::add< 3 >( elementCenters[ k ], X[ m_toNodesRelation( k, a ) ] );
      }

      LvArray::tensorOps::scale< 3 >( elementCenters[ k ], 1.0 / nNodes );
    } );
  }

  inline void CalculateCellVolumesKernel( localIndex const k,
                                          arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & X ) const
  {
    LvArray::tensorOps::fill< 3 >( m_elementCenter[ k ], 0 );

    real64 Xlocal[10][3];

    for( localIndex a = 0; a < m_numNodesPerElement; ++a )
    {
      LvArray::tensorOps::copy< 3 >( Xlocal[ a ], X[ m_toNodesRelation( k, a ) ] );
      LvArray::tensorOps::add< 3 >( m_elementCenter[ k ], X[ m_toNodesRelation( k, a ) ] );
    }
    LvArray::tensorOps::scale< 3 >( m_elementCenter[ k ], 1.0 / m_numNodesPerElement );

    if( m_numNodesPerElement == 8 )
    {
      m_elementVolume[k] = computationalGeometry::HexVolume( Xlocal );
    }
    else if( m_numNodesPerElement == 4 )
    {
      m_elementVolume[k] = computationalGeometry::TetVolume( Xlocal );
    }
    else if( m_numNodesPerElement == 6 )
    {
      m_elementVolume[k] = computationalGeometry::WedgeVolume( Xlocal );
    }
    else if( m_numNodesPerElement == 5 )
    {
      m_elementVolume[k] = computationalGeometry::PyramidVolume( Xlocal );
    }
    else
    {
      GEOSX_ERROR( "GEOX does not support cells with " << m_numNodesPerElement << " nodes" );
    }
  }



  virtual void SetElementType( string const & elementType ) override;

  localIndex GetNumFaceNodes( localIndex const elementIndex,
                              localIndex const localFaceIndex ) const;

  localIndex GetFaceNodes( localIndex const elementIndex,
                           localIndex const localFaceIndex,
                           localIndex * const nodeIndices ) const;


  void GetFaceNodes( localIndex const elementIndex,
                     localIndex const localFaceIndex,
                     localIndex_array & nodeIndices ) const;


  NodeMapType & nodeList() { return m_toNodesRelation; }

  NodeMapType const & nodeList() const { return m_toNodesRelation; }

  localIndex & nodeList( localIndex const k, localIndex a ) { return m_toNodesRelation( k, a ); }

  localIndex const & nodeList( localIndex const k, localIndex a ) const { return m_toNodesRelation( k, a ); }

  FixedOneToManyRelation & edgeList() { return m_toEdgesRelation; }

  FixedOneToManyRelation const & edgeList() const { return m_toEdgesRelation; }

  FixedOneToManyRelation & faceList() { return m_toFacesRelation; }

  FixedOneToManyRelation const & faceList() const { return m_toFacesRelation; }



  template< typename T >
  T & AddProperty( string const & propertyName )
  {
    m_externalPropertyNames.emplace_back( propertyName );
    return this->registerWrapper< T >( propertyName )->reference();
  }

  template< typename LAMBDA >
  void forExternalProperties( LAMBDA && lambda )
  {
    for( auto & externalPropertyName : m_externalPropertyNames )
    {
      dataRepository::WrapperBase * const wrapper = this->getWrapperBase( externalPropertyName );
      lambda( wrapper );
    }
  }


protected:

  NodeMapType m_toNodesRelation;

  EdgeMapType m_toEdgesRelation;

  FaceMapType m_toFacesRelation;

private:
  string_array m_externalPropertyNames;

};

}

#endif /* GEOSX_MESH_CELLBLOCK_HPP_ */