FEMHelpers.hpp

/*------------------------------------------------------------------------*/
/*                 Copyright 2010, 2011 Sandia Corporation.                     */
/*  Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive   */
/*  license for use of this work by or on behalf of the U.S. Government.  */
/*  Export of this program may require a license from the                 */
/*  United States Government.                                             */
/*------------------------------------------------------------------------*/

#ifndef stk_mesh_FEMHelpers_hpp
#define stk_mesh_FEMHelpers_hpp

#include <stk_mesh/base/Types.hpp>

#include <stk_mesh/fem/FEMMetaData.hpp>
#include <stk_mesh/fem/CellTopology.hpp>
// This is needed for ElementNode class
#include <stk_mesh/fem/TopologyDimensions.hpp>

namespace stk_classic {
namespace mesh {

class Bucket;
class Entity;

namespace fem {

//----------------------------------------------------------------------
Entity & declare_element( BulkData & mesh ,
                          Part & part ,
                          const EntityId elem_id ,
                          const EntityId node_id[] );


Entity & declare_element_side( BulkData & mesh ,
                               const stk_classic::mesh::EntityId global_side_id ,
                               Entity & elem ,
                               const unsigned local_side_id ,
                               Part * part = NULL );

Entity & declare_element_edge( BulkData & mesh ,
                               const stk_classic::mesh::EntityId global_side_id ,
                               Entity & elem ,
                               const unsigned local_side_id ,
                               Part * part = NULL );

bool element_side_polarity( const Entity & elem ,
                            const Entity & side , int local_side_id = -1 );

Entity & declare_element_side( Entity & elem ,
                               Entity & side ,
                               const unsigned local_side_id ,
                               Part * part = NULL );



Entity & declare_element_edge( Entity & elem ,
                               Entity & edge ,
                               const unsigned local_edge_id ,
                               Part * part = NULL );



template< class Top >
Part &declare_part(FEMMetaData& meta_data, const std::string &name) {
  return meta_data.declare_part(name, shards::getCellTopologyData<Top>());
}

const CellTopologyData * get_subcell_nodes(
    const Entity     & entity ,
    EntityRank         subcell_rank ,
    unsigned           subcell_identifier ,
    EntityVector     & subcell_nodes
    );

int get_entity_subcell_id( const Entity            & entity ,
                           const EntityRank          subcell_rank,
                           const CellTopologyData  * side_topology,
                           const EntityVector      & side_nodes );

bool comm_mesh_counts( BulkData & ,
                       std::vector<size_t> & counts ,
                       bool = false );

typedef Field<double*,stk_classic::mesh::ElementNode> ElementNodePointerField ;

template< class NodeField >
inline
ElementNodePointerField &
declare_element_node_pointer_field(
  FEMMetaData & fmd , const std::string & s ,
  NodeField & node_field )
{
  const unsigned num_states = node_field.number_of_states();

  ElementNodePointerField & f =
    fmd.template declare_field< ElementNodePointerField >( s, num_states );

  for ( unsigned i = 0 ; i < num_states ; ++i ) {
    FieldState state = (FieldState) i;
    fmd.declare_field_relation(
      f.field_of_state( state ) ,
      fem::get_element_node_stencil(fmd.spatial_dimension()) ,
      node_field.field_of_state( state ) );
  }

  return f ;
}

template< class Traits >
void get_parts_with_topology(stk_classic::mesh::BulkData& mesh,
                             stk_classic::mesh::PartVector& parts,
                             bool skip_topology_root_parts=false)
{
  parts.clear();

  stk_classic::mesh::fem::FEMMetaData & fem_meta = stk_classic::mesh::fem::FEMMetaData::get(mesh);

  const stk_classic::mesh::PartVector& all_parts = fem_meta.get_parts();

  stk_classic::mesh::PartVector::const_iterator
    iter = all_parts.begin(),
    iter_end = all_parts.end();

  const CellTopologyData* topology = shards::getCellTopologyData<Traits>();

  for(; iter!=iter_end; ++iter) {
    stk_classic::mesh::Part* part =  *iter;
    if (fem_meta.get_cell_topology(*part).getCellTopologyData() == topology) {
      if (skip_topology_root_parts && stk_classic::mesh::fem::is_cell_topology_root_part(*part)) {
        continue;
      }
      parts.push_back(part);
    }
  }
}

inline
unsigned get_spatial_dimension(const Entity& entity)
{
  // expose some dot-chain to ensure everything inlined
  return entity.bucket().mesh().mesh_meta_data().get_spatial_dimension();
}

/* The Fmwk uses an enum to identify nodes, edges, faces and elements. The
   toolkit is similar, but the the toolkit rank depends on the spatial
   dimension. For 3D parts they are identical. With 2D, the element rank is
   2 (not 3). etc.
*/

#ifdef SIERRA_MIGRATION

inline
unsigned convert_fmwk_rank_to_stk(unsigned fmwk_rank, unsigned spatial_dim)
{
  ThrowAssert(spatial_dim > 0);
  ThrowAssert(spatial_dim < 4);
  ThrowAssert(fmwk_rank <= 4); // up to four basic entities types and constraints
  static int map2Stk[4][5] = { {-1,-1,-1,-1, 4},
                               { 0,-1,-1, 1, 4},
                               { 0, 1,-1, 2, 4},
                               { 0, 1, 2, 3, 4}};
  int stk_rank = map2Stk[spatial_dim][fmwk_rank];

  return static_cast<unsigned>(stk_rank);
}

inline
unsigned convert_stk_rank_to_fmwk(unsigned stk_rank,  unsigned spatial_dim)
{
  ThrowAssert(spatial_dim > 0);
  ThrowAssert(spatial_dim < 4);
  ThrowAssert(stk_rank <= 4); // up to four basic entities types and constraints
  static int map2Fmwk[4][5]={ {-1,-1,-1,-1, 4},
                              { 0, 3,-1,-1, 4},
                              { 0, 1, 3,-1, 4},
                              { 0, 1, 2, 3, 4}};
  int fmwk_rank=map2Fmwk[spatial_dim][stk_rank];

  return static_cast<unsigned>(fmwk_rank);
}

inline
unsigned get_derived_type(const Entity& entity)
{
  return convert_stk_rank_to_fmwk(entity.entity_rank(), get_spatial_dimension(entity));
}

#endif

} //namespace fem
} //namespace mesh
} //namespace stk_classic
#endif