InternalMeshGenerator.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 Total, S.A
 * 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_GENERATORS_INTERNALMESHGENERATOR_HPP
#define GEOS_MESH_GENERATORS_INTERNALMESHGENERATOR_HPP
 
#include "common/format/EnumStrings.hpp"
#include "mesh/generators/MeshGeneratorBase.hpp"
#include "mesh/generators/CellBlockManager.hpp"
#include "mesh/mpiCommunications/SpatialPartition.hpp"
 
namespace geos
{
 
 
class InternalMeshGenerator : public MeshGeneratorBase
{
public:
 
  InternalMeshGenerator( const string & name, Group * const parent );
 
  virtual ~InternalMeshGenerator() override = default;
 
  static string catalogName() { return "InternalMesh"; }
 
 
  void importFieldOnArray( Block block,
                           string const & blockName,
                           string const & meshFieldName,
                           bool isMaterialField,
                           dataRepository::WrapperBase & wrapper ) const override
  {
    GEOS_UNUSED_VAR( block );
    GEOS_UNUSED_VAR( blockName );
    GEOS_UNUSED_VAR( meshFieldName );
    GEOS_UNUSED_VAR( isMaterialField );
    GEOS_UNUSED_VAR( wrapper );
  }
 
  virtual inline bool isCartesian() const
  {
    return true;
  }
 
  virtual void reduceNumNodesForPeriodicBoundary( SpatialPartition & partition,
                                                  integer (& numNodes) [3] )
  {
    GEOS_UNUSED_VAR( partition, numNodes );
  };
 
  virtual void
  setNodeGlobalIndicesOnPeriodicBoundary( SpatialPartition & partition,
                                          int (& index)[3] )
  {
    GEOS_UNUSED_VAR( partition, index );
  }
 
  virtual void setConnectivityForPeriodicBoundaries( int ( & globalIJK )[3],
                                                     integer const ( &numNodesInDir )[3],
                                                     int const ( &firstElemIndexInPartition )[3],
                                                     localIndex ( & nodeOfBox )[8] )
  {
    GEOS_UNUSED_VAR( globalIJK, numNodesInDir, firstElemIndexInPartition, nodeOfBox );
  }
 
  void setConnectivityForPeriodicBoundary( int const component,
                                           int const (&globalIJK)[3],
                                           integer const (&numNodesInDir)[3],
                                           int const (&firstElemIndexInPartition)[3],
                                           localIndex ( &nodeOfBox )[8] );
 
  virtual void coordinateTransformation( arrayView2d< real64, nodes::REFERENCE_POSITION_USD > X, std::map< string, SortedArray< localIndex > > & nodeSets )
  {
    GEOS_UNUSED_VAR( X );
    GEOS_UNUSED_VAR( nodeSets );
  }
 
 
protected:
 
  struct viewKeyStruct
  {
    constexpr static char const * xCoordsString() { return "xCoords"; }
    constexpr static char const * yCoordsString() { return "yCoords"; }
    constexpr static char const * zCoordsString() { return "zCoords"; }
    constexpr static char const * xElemsString() { return "nx"; }
    constexpr static char const * yElemsString() { return "ny"; }
    constexpr static char const * zElemsString() { return "nz"; }
    constexpr static char const * xBiasString() { return "xBias"; }
    constexpr static char const * yBiasString() { return "yBias"; }
    constexpr static char const * zBiasString() { return "zBias"; }
    constexpr static char const * cellBlockNamesString() { return "cellBlockNames"; }
    constexpr static char const * elementTypesString() { return "elementTypes"; }
    constexpr static char const * trianglePatternString() { return "trianglePattern"; }
    constexpr static char const * meshTypeString() { return "meshType"; }
    constexpr static char const * positionToleranceString() { return "positionTolerance"; }
  };
 
  void postInputInitialization() override;
 
  int m_dim;
 
  real64 m_min[3];
 
  real64 m_max[3];
 
  real64 m_coordinatePrecision;
 
  array1d< real64 > m_vertices[3];
 
  array1d< integer > m_nElems[3];
 
  array1d< real64 > m_nElemBias[3];
 
  array1d< real64 > m_setCoords[3];
 
private:
 
  array1d< string > m_regionNames;
 
  array1d< integer > m_firstElemIndexForBlock[3];
 
  array1d< integer > m_lastElemIndexForBlock[3];
 
  globalIndex m_numElemsTotal[3];
 
  array1d< string > m_elementType;
 
  array1d< integer > m_numElePerBox;
 
  int m_trianglePattern;
 
  real64 m_fPerturb = 0.0;
 
  int m_randSeed = 0;
 
  real64 m_skewAngle = 0;
 
  real64 m_skewCenter[3] = { 0, 0, 0 };
 
 
 
  virtual void fillCellBlockManager( CellBlockManager & cellBlockManager, SpatialPartition & partition ) override;
 
  inline globalIndex nodeGlobalIndex( int const index[3] )
  {
    return index[0] + index[1]*(m_numElemsTotal[0]+1) + index[2]*(m_numElemsTotal[0]+1)*(m_numElemsTotal[1]+1);
  }
 
  inline globalIndex elemGlobalIndex( int const index[3] )
  {
    return index[0] + index[1]*m_numElemsTotal[0] + index[2]*m_numElemsTotal[0]*m_numElemsTotal[1];
  }
 
  template< typename OUT_VECTOR >
  inline void getNodePosition( int const (&a)[3], int trianglePattern, OUT_VECTOR && X )
  {
    real64 xInterval( 0 );
 
    int xPosIndex = 0;
    if( trianglePattern == 1 )
    {
      int startingIndex = 0;
      int endingIndex = 0;
      int block = 0;
      for( block=0; block<m_nElems[0].size(); ++block )
      {
        startingIndex = endingIndex;
        endingIndex = startingIndex + m_nElems[0][block];
      }
      xPosIndex = endingIndex;
    }
 
    for( int i=0; i<3; ++i )
    {
      if( m_setCoords[i].size()>0 )
      {
        X[i] = m_setCoords[i][a[i]];
      }
      else
      {
 
        int startingIndex = 0;
        int endingIndex = 0;
        int block = 0;
        for( block=0; block<m_nElems[i].size(); ++block )
        {
          startingIndex = endingIndex;
          endingIndex = startingIndex + m_nElems[i][block];
          if( a[i]>=startingIndex && a[i]<=endingIndex )
          {
            break;
          }
        }
        real64 min = m_vertices[i][block];
        real64 max = m_vertices[i][block+1];
 
 
        X[i] = min + (max-min) * ( double( a[i] - startingIndex ) / m_nElems[i][block] );
 
        // First check if m_nElemBias contains values
        // Otherwise the next test will cause a segfault when looking for "block"
        if( m_nElemBias[i].size()>0 )
        {
          // Verify that the bias is non-zero and applied to more than one block:
          if( ( !isZero( m_nElemBias[i][block] ) ) && (m_nElems[i][block]>1))
          {
            GEOS_ERROR_IF( fabs( m_nElemBias[i][block] ) >= 1,
                           getWrapperDataContext( i == 0 ? viewKeyStruct::xBiasString() :
                                                  i == 1 ? viewKeyStruct::yBiasString() :
                                                  viewKeyStruct::zBiasString() ) <<
                           ", block index = " << block << " : Mesh bias must between -1 and 1!" );
 
            real64 len = max -  min;
            real64 xmean = len / m_nElems[i][block];
            real64 x0 = xmean * double( a[i] - startingIndex );
            real64 chi = m_nElemBias[i][block]/(xmean/len - 1.0);
            real64 dx = -x0*chi + x0*x0*chi/len;
            X[i] += dx;
          }
        }
 
        // This is for creating regular triangle pattern
        if( i==0 ) xInterval = (max-min) / m_nElems[i][block];
        if( trianglePattern == 1 && i == 1 && a[1] % 2 == 1 && a[0] != 0 && a[0] != xPosIndex )
          X[0] -= xInterval * 0.5;
      }
    }
  }
 
  template< typename OUT_VECTOR >
  inline void getElemCenterPosition( int const k[3], OUT_VECTOR && X )
  {
    for( int i=0; i<3; ++i )
    {
      X[i] = m_min[i] + (m_max[i]-m_min[i]) * ( ( k[i] + 0.5 ) / m_numElemsTotal[i] );
    }
  }
 
public:
 
 
};
 
} /* namespace geos */
 
#endif /* GEOS_MESH_GENERATORS_INTERNALMESHGENERATOR_HPP */