Loading [MathJax]/extensions/tex2jax.js
Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ReadUtil.cpp
Go to the documentation of this file.
1 /** 2  * MOAB, a Mesh-Oriented datABase, is a software component for creating, 3  * storing and accessing finite element mesh data. 4  * 5  * Copyright 2004 Sandia Corporation. Under the terms of Contract 6  * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government 7  * retains certain rights in this software. 8  * 9  * This library is free software; you can redistribute it and/or 10  * modify it under the terms of the GNU Lesser General Public 11  * License as published by the Free Software Foundation; either 12  * version 2.1 of the License, or (at your option) any later version. 13  * 14  */ 15  16 #ifdef WIN32 17 #ifdef _DEBUG 18 // turn off warnings that say they debugging identifier has been truncated 19 // this warning comes up when using some STL containers 20 #pragma warning( disable : 4786 ) 21 #endif 22 #endif 23  24 #include "ReadUtil.hpp" 25 #include "moab/Core.hpp" 26 #include "AEntityFactory.hpp" 27 #include "moab/Error.hpp" 28 #include "SequenceManager.hpp" 29 #include "VertexSequence.hpp" 30 #include "ElementSequence.hpp" 31  32 namespace moab 33 { 34  35 #define RR \ 36  if( MB_SUCCESS != result ) return result 37  38 ReadUtil::ReadUtil( Core* mdb, Error* /*error_handler*/ ) : ReadUtilIface(), mMB( mdb ) {} 39  40 ErrorCode ReadUtil::get_node_coords( const int /*num_arrays*/, 41  const int num_nodes, 42  const int preferred_start_id, 43  EntityHandle& actual_start_handle, 44  std::vector< double* >& arrays, 45  int sequence_size ) 46 { 47  ErrorCode error; 48  EntitySequence* seq = 0; 49  50  if( num_nodes < 1 ) 51  { 52  actual_start_handle = 0; 53  arrays.clear(); 54  return MB_INDEX_OUT_OF_RANGE; 55  } 56  57  // Create an entity sequence for these nodes 58  error = mMB->sequence_manager()->create_entity_sequence( MBVERTEX, num_nodes, 0, preferred_start_id, 59  actual_start_handle, seq, sequence_size ); 60  61  if( error != MB_SUCCESS ) return error; 62  63  if( seq->start_handle() > actual_start_handle || seq->end_handle() < actual_start_handle || 64  seq->end_handle() - actual_start_handle + 1 < (unsigned)num_nodes ) 65  return MB_FAILURE; 66  67  arrays.resize( 3 ); 68  69  error = static_cast< VertexSequence* >( seq )->get_coordinate_arrays( arrays[0], arrays[1], arrays[2] ); 70  for( unsigned i = 0; i < arrays.size(); ++i ) 71  if( arrays[i] ) arrays[i] += ( actual_start_handle - seq->start_handle() ); 72  73  return error; 74 } 75  76 ErrorCode ReadUtil::get_element_connect( const int num_elements, 77  const int verts_per_element, 78  const EntityType mdb_type, 79  const int preferred_start_id, 80  EntityHandle& actual_start_handle, 81  EntityHandle*& array, 82  int sequence_size ) 83 { 84  ErrorCode error; 85  EntitySequence* seq; 86  87  if( num_elements < 1 ) 88  { 89  actual_start_handle = 0; 90  array = 0; 91  return MB_INDEX_OUT_OF_RANGE; 92  } 93  94  // if (mdb_type <= MBVERTEX || mdb_type >= MBPOLYHEDRON || mdb_type == MBPOLYGON) 95  // return MB_TYPE_OUT_OF_RANGE; 96  97  // Make an entity sequence to hold these elements 98  error = 99  mMB->sequence_manager()->create_entity_sequence( mdb_type, num_elements, verts_per_element, preferred_start_id, 100  actual_start_handle, seq, sequence_size ); 101  if( MB_SUCCESS != error ) return error; 102  103  if( seq->start_handle() > actual_start_handle || seq->end_handle() < actual_start_handle || 104  seq->end_handle() - actual_start_handle + 1 < (unsigned)num_elements ) 105  return MB_FAILURE; 106  107  // Get an array for the connectivity 108  array = static_cast< ElementSequence* >( seq )->get_connectivity_array(); 109  if( !array ) return MB_FAILURE; 110  array += 111  ( actual_start_handle - seq->start_handle() ) * static_cast< ElementSequence* >( seq )->nodes_per_element(); 112  113  return error; 114 } 115  116 ErrorCode ReadUtil::create_entity_sets( EntityID num_sets, 117  const unsigned* flags, 118  EntityID start_id, 119  EntityHandle& start_handle ) 120 { 121  if( num_sets < 1 ) 122  { 123  start_handle = 0; 124  return MB_INDEX_OUT_OF_RANGE; 125  } 126  127  ErrorCode error; 128  EntitySequence* seq; 129  error = mMB->sequence_manager()->create_meshset_sequence( num_sets, start_id, flags, start_handle, seq ); 130  if( MB_SUCCESS != error ) return error; 131  132  if( seq->start_handle() > start_handle || seq->end_handle() < start_handle || 133  seq->end_handle() - start_handle + 1 < (EntityHandle)num_sets ) 134  return MB_FAILURE; 135  136  return MB_SUCCESS; 137 } 138  139 ErrorCode ReadUtil::update_adjacencies( const EntityHandle start_handle, 140  const int number_elements, 141  const int number_vertices_per_element, 142  const EntityHandle* conn_array ) 143 { 144  EntityHandle tmp_hndl = start_handle; 145  AEntityFactory* adj_fact = mMB->a_entity_factory(); 146  147  // Iterate over the elements and update adjacency information 148  if( adj_fact != NULL && adj_fact->vert_elem_adjacencies() ) 149  { 150  int j = 0; 151  for( int i = 0; i < number_elements; i++ ) 152  { 153  adj_fact->notify_create_entity( tmp_hndl, ( conn_array + j ), number_vertices_per_element ); 154  tmp_hndl++; 155  j += number_vertices_per_element; 156  } 157  } 158  159  return MB_SUCCESS; 160 } 161  162 ErrorCode ReadUtil::gather_related_ents( Range& partition, Range& related_ents, EntityHandle* file_set ) 163 { 164  // Loop over any sets, getting contained ents 165  std::pair< Range::const_iterator, Range::const_iterator > pair_it = partition.equal_range( MBENTITYSET ); 166  167  ErrorCode result = MB_SUCCESS; 168  for( Range::const_iterator rit = pair_it.first; rit != pair_it.second; ++rit ) 169  { 170  ErrorCode tmp_result = mMB->get_entities_by_handle( *rit, related_ents, Interface::UNION ); 171  if( MB_SUCCESS != tmp_result ) result = tmp_result; 172  } 173  RR; 174  175  // Gather adjacent ents of other dimensions 176  Range tmp_ents; 177  for( int dim = 3; dim >= 0; dim-- ) 178  { 179  tmp_ents.clear(); 180  ErrorCode tmp_result = mMB->get_adjacencies( related_ents, dim, false, tmp_ents, Interface::UNION ); 181  if( MB_SUCCESS != tmp_result ) 182  result = tmp_result; 183  else 184  related_ents.merge( tmp_ents ); 185  } 186  RR; 187  188  // Related ents includes the partition itself 189  related_ents.merge( partition ); 190  191  // Get contains-related sets 192  Range tmp_ents3, last_related; 193  if( file_set ) 194  result = mMB->get_entities_by_type( *file_set, MBENTITYSET, tmp_ents3 ); 195  else 196  result = mMB->get_entities_by_type( 0, MBENTITYSET, tmp_ents3 );RR; 197  198  while( related_ents.size() != last_related.size() ) 199  { 200  last_related = related_ents; 201  for( Range::iterator rit = tmp_ents3.begin(); rit != tmp_ents3.end(); ++rit ) 202  { 203  if( related_ents.find( *rit ) != related_ents.end() ) continue; 204  205  tmp_ents.clear(); 206  result = mMB->get_entities_by_handle( *rit, tmp_ents, true );RR; 207  Range tmp_ents2 = intersect( tmp_ents, related_ents ); 208  209  // If the intersection is not empty, set is related 210  if( !tmp_ents2.empty() ) related_ents.insert( *rit ); 211  } 212  } 213  214  // Get child-related sets 215  last_related.clear(); 216  while( related_ents.size() != last_related.size() ) 217  { 218  last_related = related_ents; 219  std::pair< Range::const_iterator, Range::const_iterator > it_pair = last_related.equal_range( MBENTITYSET ); 220  221  for( Range::const_iterator rit = it_pair.first; rit != it_pair.second; ++rit ) 222  { 223  // Get all children and add to related ents 224  tmp_ents.clear(); 225  result = mMB->get_child_meshsets( *rit, tmp_ents, 0 );RR; 226  related_ents.merge( tmp_ents ); 227  } 228  } 229  230  // Get parent-related sets 231  last_related.clear(); 232  while( related_ents.size() != last_related.size() ) 233  { 234  last_related = related_ents; 235  std::pair< Range::const_iterator, Range::const_iterator > it_pair = last_related.equal_range( MBENTITYSET ); 236  237  for( Range::const_iterator rit = it_pair.first; rit != it_pair.second; ++rit ) 238  { 239  // Get all parents and add to related ents 240  tmp_ents.clear(); 241  result = mMB->get_parent_meshsets( *rit, tmp_ents, 0 );RR; 242  related_ents.merge( tmp_ents ); 243  } 244  } 245  246  return MB_SUCCESS; 247 } 248  249 ErrorCode ReadUtil::get_ordered_vertices( EntityHandle* bound_ents, 250  int* sense, 251  int bound_size, 252  int dim, 253  EntityHandle* bound_verts, 254  EntityType& etype ) 255 { 256  // Get dimension of bounding entities 257  int bound_dim = CN::Dimension( TYPE_FROM_HANDLE( bound_ents[0] ) ); 258  int indices[MAX_SUB_ENTITY_VERTICES]; 259  const EntityHandle* connect = NULL; 260  std::vector< EntityHandle > tmp_connect; 261  262  // Find the right entity type based on # bounding ents 263  int numv = 0, num_connect = 0; 264  ErrorCode result; 265  for( EntityType t = MBEDGE; t < MBENTITYSET; t++ ) 266  { 267  int nindex = CN::NumSubEntities( t, bound_dim ); 268  if( CN::Dimension( t ) != dim || nindex != bound_size ) continue; 269  270  // Fill in vertices from bounding entity vertices 271  int nverts = CN::VerticesPerEntity( t ); 272  std::fill( bound_verts, bound_verts + nverts, 0 ); 273  for( int index = 0; index < nindex; index++ ) 274  { 275  result = mMB->get_connectivity( bound_ents[index], connect, num_connect, false, &tmp_connect ); 276  if( MB_SUCCESS != result ) return result; 277  278  CN::SubEntityVertexIndices( t, bound_dim, index, indices ); 279  280  for( int c = 0; c < num_connect; c++ ) 281  { 282  if( !bound_verts[indices[c]] ) 283  { 284  bound_verts[indices[c]] = ( sense[index] > 0 ) ? connect[c] : connect[num_connect - c - 1]; 285  numv++; 286  } 287  } 288  if( numv == nverts ) 289  { 290  etype = t; 291  return MB_SUCCESS; 292  } 293  } 294  } 295  296  // If we get here, we didn't get full connectivity 297  etype = MBMAXTYPE; 298  return MB_FAILURE; 299 } 300  301 static ErrorCode check_int_tag( Interface* mb, Tag tag ) 302 { 303  int size; 304  DataType type; 305  ErrorCode rval = mb->tag_get_bytes( tag, size ); 306  if( MB_SUCCESS != rval ) return rval; 307  if( size != sizeof( int ) ) return MB_TYPE_OUT_OF_RANGE; 308  rval = mb->tag_get_data_type( tag, type ); 309  if( type != MB_TYPE_OPAQUE && type != MB_TYPE_INTEGER ) return MB_TYPE_OUT_OF_RANGE; 310  311  return MB_SUCCESS; 312 } 313  314 ErrorCode ReadUtil::assign_ids( Tag id_tag, const Range& ents, int start ) 315 { 316  ErrorCode rval = check_int_tag( mMB, id_tag ); 317  if( MB_SUCCESS != rval ) return rval; 318  319  Range tmp_range; 320  std::vector< int > data; 321  for( Range::const_pair_iterator i = ents.pair_begin(); i != ents.pair_end(); ++i ) 322  { 323  data.resize( i->second + 1 - i->first ); 324  for( std::vector< int >::iterator j = data.begin(); j != data.end(); ++j ) 325  *j = start++; 326  tmp_range.clear(); 327  tmp_range.insert( i->first, i->second ); 328  rval = mMB->tag_set_data( id_tag, tmp_range, &data[0] ); 329  if( MB_SUCCESS != rval ) return rval; 330  } 331  332  return MB_SUCCESS; 333 } 334  335 ErrorCode ReadUtil::assign_ids( Tag id_tag, const EntityHandle* ents, size_t num_ents, int start ) 336 { 337  ErrorCode rval = check_int_tag( mMB, id_tag ); 338  if( MB_SUCCESS != rval ) return rval; 339  340  std::vector< int > data; 341  const EntityHandle* const end = ents + num_ents; 342  const EntityHandle* i = ents; 343  while( i != end ) 344  { 345  const EntityHandle* next = std::find( i, end, 0u ); 346  size_t size = next - i; 347  if( !size ) 348  { 349  ++i; 350  continue; 351  } 352  353  int id = start + ( i - ents ); 354  data.resize( size ); 355  for( std::vector< int >::iterator j = data.begin(); j != data.end(); ++j ) 356  *j = id++; 357  358  rval = mMB->tag_set_data( id_tag, i, size, &data[0] ); 359  if( MB_SUCCESS != rval ) return rval; 360  } 361  362  return MB_SUCCESS; 363 } 364  365 ErrorCode ReadUtil::create_gather_set( EntityHandle& gather_set ) 366 { 367  ErrorCode rval = mMB->create_meshset( MESHSET_SET, gather_set ); 368  if( MB_SUCCESS != rval ) return rval; 369  370  Tag gather_set_tag; 371  rval = mMB->tag_get_handle( "GATHER_SET", 1, MB_TYPE_INTEGER, gather_set_tag, MB_TAG_CREAT | MB_TAG_SPARSE ); 372  if( MB_SUCCESS != rval ) return rval; 373  374  int gather_val = 1; 375  rval = mMB->tag_set_data( gather_set_tag, &gather_set, 1, &gather_val ); 376  if( MB_SUCCESS != rval ) return rval; 377  378  return MB_SUCCESS; 379 } 380  381 ErrorCode ReadUtil::get_gather_set( EntityHandle& gather_set ) 382 { 383  Tag gather_set_tag; 384  ErrorCode rval = mMB->tag_get_handle( "GATHER_SET", 1, MB_TYPE_INTEGER, gather_set_tag, MB_TAG_SPARSE ); 385  if( MB_SUCCESS != rval ) return rval; 386  387  int gather_val = 1; 388  void* vals[] = { &gather_val }; 389  Range gather_sets; 390  rval = mMB->get_entities_by_type_and_tag( 0, MBENTITYSET, &gather_set_tag, vals, 1, gather_sets ); 391  if( MB_SUCCESS != rval ) return rval; 392  393  if( gather_sets.empty() ) return MB_ENTITY_NOT_FOUND; 394  395  gather_set = gather_sets[0]; 396  397  return MB_SUCCESS; 398 } 399  400 } // namespace moab