Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
Core.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 // turn off warnings that say they debugging identifier has been truncated
18 // this warning comes up when using some STL containers
19 #pragma warning( disable : 4786 )
20 #endif
21 
22 #define MOAB_MPE_LOG "moab.mpe"
23 
24 #include <iostream>
25 #include <sstream>
26 #include <vector>
27 #include <string>
28 #include <algorithm>
29 #include "moab/Core.hpp"
30 #include "MeshSetSequence.hpp"
31 #include "ElementSequence.hpp"
32 #include "VertexSequence.hpp"
33 #include <cassert>
34 #include "AEntityFactory.hpp"
35 #include "ReadUtil.hpp"
36 #include "WriteUtil.hpp"
37 #include "moab/CN.hpp"
39 #include "SequenceManager.hpp"
40 #include "moab/Error.hpp"
41 #include "moab/ReaderWriterSet.hpp"
42 #include "moab/ReaderIface.hpp"
43 #include "moab/WriterIface.hpp"
44 #include "moab/ScdInterface.hpp"
45 #include "moab/SetIterator.hpp"
46 
47 #include "BitTag.hpp"
48 #include "DenseTag.hpp"
49 #include "MeshTag.hpp"
50 #include "SparseTag.hpp"
51 #include "VarLenDenseTag.hpp"
52 #include "VarLenSparseTag.hpp"
53 
54 #include <sys/stat.h>
55 #include <cerrno>
56 #include <cstring>
57 
58 #ifdef MOAB_HAVE_AHF
59 #include "moab/HalfFacetRep.hpp"
60 #endif
61 
62 #ifdef MOAB_HAVE_MPI
63 /* Leave ParallelComm.hpp before mpi.h or MPICH2 will fail
64  * because its C++ headers do not like SEEK_* macros.
65  */
66 #include "moab/ParallelComm.hpp"
67 #include "moab_mpi.h"
68 #include "ReadParallel.hpp"
69 #endif
70 
71 #ifdef MOAB_HAVE_HDF5
72 #ifdef MOAB_HAVE_HDF5_PARALLEL
73 #include "WriteHDF5Parallel.hpp"
75 #define DefaultWriterName "WriteHDF5Parallel"
76 #else
77 #include "WriteHDF5.hpp"
79 #define DefaultWriterName "WriteHDF5"
80 #endif
81 #elif defined( MOAB_HAVE_NETCDF )
82 #include "WriteNCDF.hpp"
84 #define DefaultWriterName "WriteNCDF"
85 #else
86 #include "WriteVtk.hpp"
88 #define DefaultWriterName "WriteVtk"
89 #endif
90 #include "MBTagConventions.hpp"
91 #include "ExoIIUtil.hpp"
92 #include "EntitySequence.hpp"
93 #include "moab/FileOptions.hpp"
94 #ifdef LINUX
95 #include <dlfcn.h>
96 #include <dirent.h>
97 #endif
98 
99 #ifdef MOAB_HAVE_MPI
100 #include "moab_mpe.h"
101 #endif
102 
103 // MOAB used to use a NULL handle list and a zero handle count
104 // to indicate that tag functions are to operate on the global/mesh
105 // value of the tag. For the 3.0 release MOAB also accepted a handle
106 // with a value of 0 to indicate the same. Now we want to drop the
107 // old NULL list method because a) it is one less special case that
108 // must be handled in the tag get/set paths and b) it aviods unexpected
109 // segfaults when applications requested tag data with an empty list
110 // of handles.
111 //
112 // Define this constant to revert to the old behavior, but also print
113 // a warning.
114 #define ALLOW_NULL_FOR_MESH_TAG
115 //
116 // Define this to print an error and abort() when a NULL handle list
117 // is passed. This is intended as an interim solution to help catch
118 // spots in code that haven't been updated for the change.
119 #undef DISALLOW_EMPTY_HANDLE_LIST_FOR_TAGS
120 //
121 // The eventual goal is to define neither of the above, eliminating the
122 // check and allowing applications to safely request tag data for no
123 // entities.
124 
126 {
127  std::cerr << "WARNING: Accepting empty array to indicate mesh tag" << std::endl;
128 }
129 
130 #ifdef ALLOW_NULL_FOR_MESH_TAG
131 #define CHECK_MESH_NULL \
132  EntityHandle root = 0; \
133  if( NULL == entity_handles && 0 == num_entities ) \
134  { \
135  entity_handles = &root; \
136  num_entities = 1; \
137  warn_null_array_mesh_tag(); \
138  }
139 #elif defined( DISALLOW_EMPTY_HANDLE_LIST_FOR_TAGS )
140 #define CHECK_MESH_NULL \
141  if( NULL == entity_handles ) \
142  { \
143  std::cerr << "ERROR: Deprecated NULL handle list at " __FILE__ ":" << __LINE__ << std::endl; \
144  abort(); \
145  }
146 #else
147 #define CHECK_MESH_NULL
148 #endif
149 
150 namespace moab
151 {
152 
153 using namespace std;
154 
155 static inline const MeshSet* get_mesh_set( const SequenceManager* sm, EntityHandle h )
156 {
157  const EntitySequence* seq;
158  if( MBENTITYSET != TYPE_FROM_HANDLE( h ) || MB_SUCCESS != sm->find( h, seq ) ) return 0;
159  return reinterpret_cast< const MeshSetSequence* >( seq )->get_set( h );
160 }
161 
163 {
164  EntitySequence* seq;
165  if( MBENTITYSET != TYPE_FROM_HANDLE( h ) || MB_SUCCESS != sm->find( h, seq ) ) return 0;
166  return reinterpret_cast< MeshSetSequence* >( seq )->get_set( h );
167 }
168 
169 //! Constructor
171 {
172  if( initialize() != MB_SUCCESS )
173  {
174  printf( "Error initializing moab::Core\n" );
175  exit( 1 );
176  }
177 }
178 
179 //! destructor
181 {
182  if( mMBWriteUtil ) delete mMBWriteUtil;
183  if( mMBReadUtil ) delete mMBReadUtil;
184  if( scdInterface ) delete scdInterface;
185 
186  mMBWriteUtil = NULL;
187  mMBReadUtil = NULL;
188  scdInterface = NULL;
189 
190  deinitialize();
191 }
192 
194 {
195 #ifdef MOAB_HAVE_MPI
196  int flag;
197  if( MPI_SUCCESS == MPI_Initialized( &flag ) )
198  {
199  if( flag )
200  {
201  writeMPELog = !MPE_Initialized_logging();
202  if( writeMPELog ) (void)MPE_Init_log();
203  }
204  }
205 #endif
206 
207  initErrorHandlerInCore = false;
209  {
211  initErrorHandlerInCore = true;
212  }
213 
214  geometricDimension = 3;
215  materialTag = 0;
216  neumannBCTag = 0;
217  dirichletBCTag = 0;
218  geomDimensionTag = 0;
219  globalIdTag = 0;
220 
221  sequenceManager = new( std::nothrow ) SequenceManager;
222  if( !sequenceManager ) return MB_MEMORY_ALLOCATION_FAILED;
223 
224  aEntityFactory = new( std::nothrow ) AEntityFactory( this );
225  if( !aEntityFactory ) return MB_MEMORY_ALLOCATION_FAILED;
226 
227  mError = new( std::nothrow ) Error;
228  if( !mError ) return MB_MEMORY_ALLOCATION_FAILED;
229 
230  mMBWriteUtil = NULL;
231  mMBReadUtil = NULL;
232  scdInterface = NULL;
233 
234  // Readers and writers try to get pointers to above utils.
235  // Do this after pointers are initialized. (Pointers should
236  // really be initialized in constructor to avoid this kind
237  // of thing -- j.kraftcheck.)
238  readerWriterSet = new( std::nothrow ) ReaderWriterSet( this );
239  if( !readerWriterSet ) return MB_MEMORY_ALLOCATION_FAILED;
240 
241  material_tag();
242  neumannBC_tag();
243  dirichletBC_tag();
244  geom_dimension_tag();
245  globalId_tag();
246 
247 #ifdef MOAB_HAVE_AHF
248  ahfRep = new HalfFacetRep( this );
249  if( !ahfRep ) return MB_MEMORY_ALLOCATION_FAILED;
250  mesh_modified = false;
251 #endif
252 
253  return MB_SUCCESS;
254 }
255 
257 {
258  return 0;
259 }
260 
262 {
263 
264 #ifdef MOAB_HAVE_MPI
265  std::vector< ParallelComm* > pc_list;
266  ParallelComm::get_all_pcomm( this, pc_list );
267  for( std::vector< ParallelComm* >::iterator vit = pc_list.begin(); vit != pc_list.end(); ++vit )
268  delete *vit;
269 #endif
270 
271 #ifdef MOAB_HAVE_AHF
272  delete ahfRep;
273  ahfRep = 0;
274 #endif
275 
276  if( aEntityFactory ) delete aEntityFactory;
277 
278  aEntityFactory = 0;
279 
280  while( !tagList.empty() )
281  tag_delete( tagList.front() );
282 
283  if( sequenceManager ) delete sequenceManager;
284 
285  sequenceManager = 0;
286 
287  delete readerWriterSet;
288  readerWriterSet = 0;
289 
290  if( mError ) delete mError;
291  mError = 0;
292 
293 #ifdef MOAB_HAVE_MPI
294  if( writeMPELog )
295  {
296  const char* default_log = MOAB_MPE_LOG;
297  const char* logfile = getenv( "MPE_LOG_FILE" );
298  if( !logfile ) logfile = default_log;
299  MPE_Finish_log( logfile );
300  }
301 #endif
302 
303  if( initErrorHandlerInCore ) MBErrorHandler_Finalize();
304 }
305 
306 ErrorCode Core::query_interface_type( const std::type_info& interface_type, void*& ptr )
307 {
308  if( interface_type == typeid( ReadUtilIface ) )
309  {
310  if( !mMBReadUtil ) mMBReadUtil = new ReadUtil( this, mError );
311  ptr = static_cast< ReadUtilIface* >( mMBReadUtil );
312  }
313  else if( interface_type == typeid( WriteUtilIface ) )
314  {
315  if( !mMBWriteUtil ) mMBWriteUtil = new WriteUtil( this );
316  ptr = static_cast< WriteUtilIface* >( mMBWriteUtil );
317  }
318  else if( interface_type == typeid( ReaderWriterSet ) )
319  {
320  ptr = reader_writer_set();
321  }
322  else if( interface_type == typeid( Error ) )
323  {
324  ptr = mError;
325  }
326  else if( interface_type == typeid( ExoIIInterface ) )
327  {
328  ptr = static_cast< ExoIIInterface* >( new ExoIIUtil( this ) );
329  }
330  else if( interface_type == typeid( ScdInterface ) )
331  {
332  if( !scdInterface ) scdInterface = new ScdInterface( this );
333  ptr = scdInterface;
334  }
335  else
336  {
337  ptr = 0;
338  return MB_FAILURE;
339  }
340  return MB_SUCCESS;
341 }
342 
343 ErrorCode Core::release_interface_type( const std::type_info& interface_type, void* iface )
344 {
345  if( interface_type == typeid( ExoIIInterface ) )
346  delete static_cast< ExoIIInterface* >( iface );
347  else if( interface_type != typeid( ReadUtilIface ) && interface_type != typeid( WriteUtilIface ) &&
348  interface_type != typeid( ReaderWriterSet ) && interface_type != typeid( Error ) &&
349  interface_type != typeid( ScdInterface ) )
350  return MB_FAILURE;
351 
352  return MB_SUCCESS;
353 }
354 
356 {
357  *iface = 0;
358  if( uuid == IDD_MBUnknown ) *iface = this;
359  if( uuid == IDD_MBCore )
360  *iface = this;
361  else
362  return 0;
363  return 1;
364 }
365 
366 float Core::impl_version( std::string* version_string )
367 {
368  if( version_string ) *version_string = MOAB_VERSION_STRING;
369 
370  return MOAB_VERSION_MAJOR + MOAB_VERSION_MINOR / 100.0f;
371 }
372 
373 //! get the type from a handle, returns type
374 EntityType Core::type_from_handle( const EntityHandle handle ) const
375 {
376  if( !handle ) // root set
377  return MBENTITYSET;
378  else
379  return TYPE_FROM_HANDLE( handle );
380 }
381 
382 //! get the id from a handle, returns id
384 {
385  return ID_FROM_HANDLE( handle );
386 }
387 
388 //! get a handle from an id and type
389 ErrorCode Core::handle_from_id( const EntityType entity_type, const EntityID id, EntityHandle& handle ) const
390 {
391  int err;
392  handle = CREATE_HANDLE( entity_type, id, err );
393 
394  // check to see if handle exists
395  const EntitySequence* dummy_seq = 0;
396  ErrorCode error_code = sequence_manager()->find( handle, dummy_seq );
397  return error_code;
398 }
399 
400 int Core::dimension_from_handle( const EntityHandle handle ) const
401 {
402  if( !handle ) // root set
403  return 4;
404  else
405  return CN::Dimension( TYPE_FROM_HANDLE( handle ) );
406 }
407 
408 //! load mesh from data in file
409 //! NOTE: if there is mesh already present, the new mesh will be added
410 ErrorCode Core::load_mesh( const char* file_name, const int* block_id_list, const int num_blocks )
411 {
412  const char* name = block_id_list ? MATERIAL_SET_TAG_NAME : 0;
413  return load_file( file_name, 0, 0, name, block_id_list, num_blocks );
414 }
415 
416 ErrorCode Core::load_file( const char* file_name,
417  const EntityHandle* file_set,
418  const char* setoptions,
419  const char* set_tag_name,
420  const int* set_tag_vals,
421  int num_set_tag_vals )
422 {
423  FileOptions opts( setoptions );
424  ErrorCode rval;
425  ReaderIface::IDTag t = { set_tag_name, set_tag_vals, num_set_tag_vals };
426  ReaderIface::SubsetList sl = { &t, 1, 0, 0 };
427 
428  assert( !file_set || ( *file_set && is_valid( *file_set ) ) );
429  if( file_set && !*file_set )
430  {
431  MB_SET_GLB_ERR( MB_FAILURE, "Non-NULL file set pointer should point to non-NULL set" );
432  }
433 
434  // if reading in parallel, call a different reader
435  std::string parallel_opt;
436  rval = opts.get_option( "PARALLEL", parallel_opt );
437  if( MB_SUCCESS == rval )
438  {
439 #ifdef MOAB_HAVE_MPI
440  ParallelComm* pcomm = 0;
441  int pcomm_id;
442  rval = opts.get_int_option( "PARALLEL_COMM", pcomm_id );
443  if( MB_ENTITY_NOT_FOUND == rval ) rval = opts.get_int_option( "PCOMM", pcomm_id );
444  if( rval == MB_SUCCESS )
445  {
446  pcomm = ParallelComm::get_pcomm( this, pcomm_id );
447  if( !pcomm ) return MB_ENTITY_NOT_FOUND;
448  }
449  else if( rval != MB_ENTITY_NOT_FOUND )
450  return rval;
451  if( set_tag_name && num_set_tag_vals )
452  {
453  rval = ReadParallel( this, pcomm ).load_file( file_name, file_set, opts, &sl );MB_CHK_ERR( rval );
454  }
455  else
456  {
457  rval = ReadParallel( this, pcomm ).load_file( file_name, file_set, opts );MB_CHK_ERR( rval );
458  }
459 #else
460  MB_SET_GLB_ERR( MB_FAILURE, "PARALLEL option not valid, this instance compiled for serial execution" );
461 #endif
462  }
463  else
464  {
465  if( set_tag_name && num_set_tag_vals )
466  {
467  rval = serial_load_file( file_name, file_set, opts, &sl );MB_CHK_ERR( rval );
468  }
469  else
470  {
471  rval = serial_load_file( file_name, file_set, opts );MB_CHK_ERR( rval );
472  }
473  }
474 
475  if( MB_SUCCESS == rval && !opts.all_seen() )
476  {
477  std::string bad_opt;
478  if( MB_SUCCESS == opts.get_unseen_option( bad_opt ) )
479  {
480  MB_SET_ERR( MB_UNHANDLED_OPTION, "Unrecognized option: \"" << bad_opt << "\"" );
481  }
482  else
483  {
484  MB_SET_ERR( MB_UNHANDLED_OPTION, "Unrecognized option" );
485  }
486  }
487 
488  return MB_SUCCESS;
489 }
490 
491 void Core::clean_up_failed_read( const Range& initial_ents, std::vector< Tag > initial_tags )
492 {
493  Range new_ents;
494  get_entities_by_handle( 0, new_ents );
495  new_ents = subtract( new_ents, initial_ents );
496  delete_entities( new_ents );
497 
498  std::vector< Tag > all_tags, new_tags;
499  tag_get_tags( all_tags );
500  std::sort( initial_tags.begin(), initial_tags.end() );
501  std::sort( all_tags.begin(), all_tags.end() );
502  std::set_difference( all_tags.begin(), all_tags.end(), initial_tags.begin(), initial_tags.end(),
503  std::back_inserter( new_tags ) );
504  while( !new_tags.empty() )
505  {
506  tag_delete( new_tags.back() );
507  new_tags.pop_back();
508  }
509 }
510 
511 ErrorCode Core::serial_load_file( const char* file_name,
512  const EntityHandle* file_set,
513  const FileOptions& opts,
514  const ReaderIface::SubsetList* subsets,
515  const Tag* id_tag )
516 {
517  int status;
518 #if defined( WIN32 ) || defined( WIN64 ) || defined( MSC_VER )
519  struct _stat stat_data;
520  status = _stat( file_name, &stat_data );
521 #else
522  struct stat stat_data;
523  status = stat( file_name, &stat_data );
524 #endif
525  if( status )
526  {
527  MB_SET_GLB_ERR( MB_FILE_DOES_NOT_EXIST, file_name << ": " << strerror( errno ) );
528  }
529 #if defined( WIN32 ) || defined( WIN64 ) || defined( MSC_VER )
530  else if( stat_data.st_mode & _S_IFDIR )
531  {
532 #else
533  else if( S_ISDIR( stat_data.st_mode ) )
534  {
535 #endif
536  MB_SET_GLB_ERR( MB_FILE_DOES_NOT_EXIST, file_name << ": Cannot read directory/folder" );
537  }
538 
539  const ReaderWriterSet* set = reader_writer_set();
540 
541  Range initial_ents;
542  ErrorCode rval = get_entities_by_handle( 0, initial_ents );MB_CHK_ERR( rval );
543 
544  std::vector< Tag > initial_tags;
545  rval = tag_get_tags( initial_tags );MB_CHK_ERR( rval );
546 
547  // otherwise try using the file extension to select a reader
548  std::string ext = set->extension_from_filename( file_name );
549 
550  // Try all the readers
552  rval = MB_FAILURE;
553  bool tried_one = false;
554  for( iter = set->begin(); iter != set->end(); ++iter )
555  {
556  if( !iter->reads_extension( ext.c_str() ) ) continue;
557 
558  ReaderIface* reader = iter->make_reader( this );
559  if( NULL != reader )
560  {
561  tried_one = true;
562  rval = reader->load_file( file_name, file_set, opts, subsets, id_tag );
563  delete reader;
564  if( MB_SUCCESS == rval ) break;
565  clean_up_failed_read( initial_ents, initial_tags );
566  }
567  }
568 
569  if( MB_SUCCESS != rval && !tried_one )
570  {
571  // didn't recognize the extension; try all of them now
572  for( iter = set->begin(); iter != set->end(); ++iter )
573  {
574  ReaderIface* reader = iter->make_reader( this );
575  if( !reader ) continue;
576  rval = reader->load_file( file_name, file_set, opts, subsets, id_tag );
577  delete reader;
578  if( MB_SUCCESS == rval )
579  break;
580  else
581  clean_up_failed_read( initial_ents, initial_tags );
582  }
583  }
584 
585  if( MB_SUCCESS != rval )
586  {
587  clean_up_failed_read( initial_ents, initial_tags );
588  MB_SET_ERR( rval, "Failed to load file after trying all possible readers" );
589  }
590  else if( file_set )
591  {
592  Range new_ents;
593  get_entities_by_handle( 0, new_ents );
594  new_ents = subtract( new_ents, initial_ents );
595 
596  // Check if gather set exists
597  EntityHandle gather_set;
598  rval = mMBReadUtil->get_gather_set( gather_set );
599  if( MB_SUCCESS == rval )
600  {
601  // Exclude gather set itself
602  new_ents.erase( gather_set );
603 
604  // Exclude gather set entities
605  Range gather_ents;
606  rval = get_entities_by_handle( gather_set, gather_ents );
607  if( MB_SUCCESS == rval ) new_ents = subtract( new_ents, gather_ents );
608  }
609 
610  rval = add_entities( *file_set, new_ents );
611  }
612 
613  return rval;
614 } // namespace moab
615 
616 ErrorCode Core::serial_read_tag( const char* file_name,
617  const char* tag_name,
618  const FileOptions& opts,
619  std::vector< int >& vals,
620  const ReaderIface::SubsetList* subsets )
621 {
622  ErrorCode rval = MB_FAILURE;
623  const ReaderWriterSet* set = reader_writer_set();
624 
625  // otherwise try using the file extension to select a reader
626  ReaderIface* reader = set->get_file_extension_reader( file_name );
627  if( reader )
628  {
629  rval = reader->read_tag_values( file_name, tag_name, opts, vals, subsets );
630  delete reader;
631  }
632  else
633  {
634  // Try all the readers
636  for( iter = set->begin(); iter != set->end(); ++iter )
637  {
638  reader = iter->make_reader( this );
639  if( NULL != reader )
640  {
641  rval = reader->read_tag_values( file_name, tag_name, opts, vals, subsets );
642  delete reader;
643  if( MB_SUCCESS == rval ) break;
644  }
645  }
646  }
647 
648  return rval;
649 }
650 
651 ErrorCode Core::write_mesh( const char* file_name, const EntityHandle* output_list, const int num_sets )
652 {
653  return write_file( file_name, 0, 0, output_list, num_sets );
654 }
655 
656 ErrorCode Core::write_file( const char* file_name,
657  const char* file_type,
658  const char* options_string,
659  const EntityHandle* output_sets,
660  int num_output_sets,
661  const Tag* tag_list,
662  int num_tags )
663 {
664  Range range;
665  std::copy( output_sets, output_sets + num_output_sets, range_inserter( range ) );
666  return write_file( file_name, file_type, options_string, range, tag_list, num_tags );
667 }
668 
669 ErrorCode Core::write_file( const char* file_name,
670  const char* file_type,
671  const char* options_string,
672  const Range& output_sets,
673  const Tag* tag_list,
674  int num_tags )
675 {
676  // convert range to vector
677  std::vector< EntityHandle > list( output_sets.size() );
678  std::copy( output_sets.begin(), output_sets.end(), list.begin() );
679 
680  // parse some options
681  FileOptions opts( options_string );
682  ErrorCode rval;
683 
684  rval = opts.get_null_option( "CREATE" );
685  if( rval == MB_TYPE_OUT_OF_RANGE )
686  {
687  MB_SET_GLB_ERR( MB_FAILURE, "Unexpected value for CREATE option" );
688  }
689  bool overwrite = ( rval == MB_ENTITY_NOT_FOUND );
690 
691  // Get the file writer
692  std::string ext = ReaderWriterSet::extension_from_filename( file_name );
693  std::vector< std::string > qa_records;
694  const EntityHandle* list_ptr = list.empty() ? (EntityHandle*)0 : &list[0];
695 
696  rval = MB_TYPE_OUT_OF_RANGE;
697 
698  // Try all possible writers
699  for( ReaderWriterSet::iterator i = reader_writer_set()->begin(); i != reader_writer_set()->end(); ++i )
700  {
701 
702  if( ( file_type && !i->name().compare( file_type ) ) || i->writes_extension( ext.c_str() ) )
703  {
704 
705  WriterIface* writer = i->make_writer( this );
706 
707  // write the file
708  rval =
709  writer->write_file( file_name, overwrite, opts, list_ptr, list.size(), qa_records, tag_list, num_tags );
710  delete writer;
711  if( MB_SUCCESS == rval ) break;
712  printf( "Writer with name %s for file %s using extension %s (file type \"%s\") was "
713  "unsuccessful\n",
714  i->name().c_str(), file_name, ext.c_str(), file_type );
715  }
716  }
717 
718  if( file_type && rval == MB_TYPE_OUT_OF_RANGE )
719  {
720  MB_SET_ERR( rval, "Unrecognized file type \"" << file_type << "\"" );
721  }
722  // Should we use default writer (e.g. HDF5)?
723  else if( MB_SUCCESS != rval )
724  {
725  DefaultWriter writer( this );
726  printf( "Using default writer %s for file %s \n", DefaultWriterName, file_name );
727  rval = writer.write_file( file_name, overwrite, opts, list_ptr, list.size(), qa_records, tag_list, num_tags );
728  }
729 
730  if( MB_SUCCESS == rval && !opts.all_seen() )
731  {
732  std::string bad_opt;
733  if( MB_SUCCESS == opts.get_unseen_option( bad_opt ) )
734  {
735  MB_SET_ERR( MB_UNHANDLED_OPTION, "Unrecognized option: \"" << bad_opt << "\"" );
736  }
737  else
738  {
739  MB_SET_ERR( MB_UNHANDLED_OPTION, "Unrecognized option" );
740  }
741  }
742 
743  return MB_SUCCESS;
744 }
745 
746 //! deletes all mesh entities from this datastore
748 {
749 
750  ErrorCode result = MB_SUCCESS;
751 
752  // perform all deinitialization procedures to clean up
753  if( aEntityFactory ) delete aEntityFactory;
754  aEntityFactory = new AEntityFactory( this );
755 
756  for( std::list< TagInfo* >::iterator i = tagList.begin(); i != tagList.end(); ++i )
757  {
758  result = ( *i )->release_all_data( sequenceManager, mError, false );MB_CHK_ERR( result );
759  }
760 
761  sequenceManager->clear();
762 
763  return MB_SUCCESS;
764 }
765 
766 //! get overall geometric dimension
768 {
769  dim = geometricDimension;
770  return MB_SUCCESS;
771 }
772 
773 //! set overall geometric dimension
774 /** Returns error if setting to 3 dimensions, mesh has been created, and
775  * there are only 2 dimensions on that mesh
776  */
778 {
779  // check to see if current dimension is smaller
780  if( geometricDimension < dim )
781  {
782  // need to check the number of entities
783  int num;
784  /*ErrorCode result = */ get_number_entities_by_dimension( 0, geometricDimension, num );
785 
786  // test written to be more readable but possibly less efficient
787  // if (MB_SUCCESS != result) return MB_FAILURE;
788  // else if (0 != num && dim == 2 && ycoordTag == 0) return MB_FAILURE;
789  // else if (0 != num && dim == 3 && (ycoordTag == 0 || zcoordTag == 0)) return MB_FAILURE;
790  // TODO -- replace this with not using xcoordTag, etc...
791  }
792 
793  // if we got here, it's ok to set dimension
794  geometricDimension = dim;
795  return MB_SUCCESS;
796 }
797 
798 //! get blocked vertex coordinates for all vertices
799 /** Blocked = all x, then all y, etc.
800  */
801 ErrorCode Core::get_vertex_coordinates( std::vector< double >& coords ) const
802 {
803  // INEFFICIENT implementation for now, until we get blocked tag access
804  Range vertices;
805  ErrorCode result = get_entities_by_type( 0, MBVERTEX, vertices );MB_CHK_ERR( result );
806 
807  // the least we can do is resize the vector and only go through the
808  // vertex list once
809  int num_verts = vertices.size();
810  int vec_pos = 0;
811  double xyz[3];
812  coords.resize( geometricDimension * num_verts );
813  for( Range::iterator it = vertices.begin(); it != vertices.end(); ++it )
814  {
815  result = get_coords( &( *it ), 1, xyz );MB_CHK_ERR( result );
816 
817  coords[vec_pos] = xyz[0];
818  coords[num_verts + vec_pos] = xyz[1];
819  coords[2 * num_verts + vec_pos] = xyz[2];
820 
821  vec_pos++;
822  }
823 
824  return MB_SUCCESS;
825 }
826 
829  double*& xcoords_ptr,
830  double*& ycoords_ptr,
831  double*& zcoords_ptr,
832  int& count )
833 {
834  EntitySequence* seq;
835  ErrorCode rval = sequence_manager()->find( *iter, seq );
836  if( MB_SUCCESS != rval )
837  {
838  xcoords_ptr = ycoords_ptr = zcoords_ptr = NULL;
839  MB_SET_ERR( rval, "Couldn't find sequence for start handle" );
840  }
841  VertexSequence* vseq = dynamic_cast< VertexSequence* >( seq );
842  if( !vseq )
843  {
844  MB_SET_ERR( MB_ENTITY_NOT_FOUND, "Couldn't find sequence for start handle" );
845  }
846 
847  unsigned int offset = *iter - vseq->data()->start_handle();
848  xcoords_ptr = reinterpret_cast< double* >( vseq->data()->get_sequence_data( 0 ) ) + offset;
849  ycoords_ptr = reinterpret_cast< double* >( vseq->data()->get_sequence_data( 1 ) ) + offset;
850  zcoords_ptr = reinterpret_cast< double* >( vseq->data()->get_sequence_data( 2 ) ) + offset;
851 
852  EntityHandle real_end = std::min( seq->end_handle(), *( iter.end_of_block() ) );
853  if( *end ) real_end = std::min( real_end, *end );
854  count = real_end - *iter + 1;
855 
856  return MB_SUCCESS;
857 }
858 
859 ErrorCode Core::get_coords( const Range& entities, double* coords ) const
860 {
861  const TypeSequenceManager& vert_data = sequence_manager()->entity_map( MBVERTEX );
863 
864  Range::const_pair_iterator i = entities.const_pair_begin();
865  EntityHandle first = i->first;
866  while( i != entities.const_pair_end() && TYPE_FROM_HANDLE( i->first ) == MBVERTEX )
867  {
868 
869  seq_iter = vert_data.lower_bound( first );
870  if( seq_iter == vert_data.end() || first < ( *seq_iter )->start_handle() ) return MB_ENTITY_NOT_FOUND;
871  const VertexSequence* vseq = reinterpret_cast< const VertexSequence* >( *seq_iter );
872 
873  EntityID offset = first - vseq->start_handle();
874  EntityID count;
875  if( i->second <= vseq->end_handle() )
876  {
877  count = i->second - first + 1;
878  ++i;
879  if( i != entities.const_pair_end() ) first = i->first;
880  }
881  else
882  {
883  count = vseq->end_handle() - first + 1;
884  first = vseq->end_handle() + 1;
885  }
886 
887  double const *x, *y, *z;
888  ErrorCode rval = vseq->get_coordinate_arrays( x, y, z );MB_CHK_ERR( rval );
889  x += offset;
890  y += offset;
891  z += offset;
892  for( EntityID j = 0; j < count; ++j )
893  {
894  coords[3 * j] = x[j];
895  coords[3 * j + 1] = y[j];
896  coords[3 * j + 2] = z[j];
897  }
898  coords = &coords[3 * count];
899  }
900 
901  // for non-vertices...
902  ErrorCode rval = MB_SUCCESS;
903  for( Range::const_iterator rit( &( *i ), i->first ); rit != entities.end(); ++rit )
904  {
905  rval = get_coords( &( *rit ), 1, coords );MB_CHK_ERR( rval );
906  coords += 3;
907  }
908 
909  return rval;
910 }
911 
912 /**\author Jason Kraftcheck <[email protected]> - 2007-5-15 */
913 ErrorCode Core::get_coords( const Range& entities, double* x_coords, double* y_coords, double* z_coords ) const
914 {
915  const TypeSequenceManager& vert_data = sequence_manager()->entity_map( MBVERTEX );
917 
918  Range::const_pair_iterator i = entities.const_pair_begin();
919  EntityHandle first = i->first;
920  while( i != entities.const_pair_end() && TYPE_FROM_HANDLE( i->first ) == MBVERTEX )
921  {
922 
923  seq_iter = vert_data.lower_bound( first );
924  if( seq_iter == vert_data.end() || first < ( *seq_iter )->start_handle() ) return MB_ENTITY_NOT_FOUND;
925  const VertexSequence* vseq = reinterpret_cast< const VertexSequence* >( *seq_iter );
926 
927  EntityID offset = first - vseq->start_handle();
928  EntityID count;
929  if( i->second <= vseq->end_handle() )
930  {
931  count = i->second - first + 1;
932  ++i;
933  if( i != entities.const_pair_end() ) first = i->first;
934  }
935  else
936  {
937  count = vseq->end_handle() - first + 1;
938  first = vseq->end_handle() + 1;
939  }
940 
941  double const *x, *y, *z;
942  ErrorCode rval = vseq->get_coordinate_arrays( x, y, z );MB_CHK_ERR( rval );
943  if( x_coords )
944  {
945  memcpy( x_coords, x + offset, count * sizeof( double ) );
946  x_coords += count;
947  }
948  if( y_coords )
949  {
950  memcpy( y_coords, y + offset, count * sizeof( double ) );
951  y_coords += count;
952  }
953  if( z_coords )
954  {
955  memcpy( z_coords, z + offset, count * sizeof( double ) );
956  z_coords += count;
957  }
958  }
959 
960  // for non-vertices...
961  ErrorCode rval = MB_SUCCESS;
962  double xyz[3];
963  for( Range::const_iterator rit( &( *i ), i->first ); rit != entities.end(); ++rit )
964  {
965  rval = get_coords( &( *rit ), 1, xyz );MB_CHK_ERR( rval );
966  *x_coords++ = xyz[0];
967  *y_coords++ = xyz[1];
968  *z_coords++ = xyz[2];
969  }
970 
971  return rval;
972 }
973 
974 ErrorCode Core::get_coords( const EntityHandle* entities, const int num_entities, double* coords ) const
975 {
976  const EntitySequence* seq = NULL;
977  const VertexSequence* vseq = NULL;
978  const EntityHandle* const end = entities + num_entities;
979  const EntityHandle* iter = entities;
980  ErrorCode status = MB_SUCCESS;
981 
982  while( iter != end )
983  {
984  if( TYPE_FROM_HANDLE( *iter ) == MBVERTEX )
985  {
986  if( !seq )
987  {
988  seq = sequence_manager()->get_last_accessed_sequence( MBVERTEX );
989  vseq = static_cast< const VertexSequence* >( seq );
990  }
991  if( !vseq )
992  return MB_ENTITY_NOT_FOUND;
993  else if( vseq->start_handle() > *iter || vseq->end_handle() < *iter )
994  {
995  if( MB_SUCCESS != sequence_manager()->find( *iter, seq ) ) return MB_ENTITY_NOT_FOUND;
996  vseq = static_cast< const VertexSequence* >( seq );
997  }
998  vseq->get_coordinates( *iter, coords );
999  }
1000  else
1001  {
1002  static std::vector< EntityHandle > dum_conn( CN::MAX_NODES_PER_ELEMENT );
1003  static std::vector< double > dum_pos( 3 * CN::MAX_NODES_PER_ELEMENT );
1004  static const EntityHandle* conn;
1005  static int num_conn;
1006  status = get_connectivity( *iter, conn, num_conn, false, &dum_conn );MB_CHK_ERR( status );
1007  status = get_coords( conn, num_conn, &dum_pos[0] );MB_CHK_ERR( status );
1008  coords[0] = coords[1] = coords[2] = 0.0;
1009  for( int i = 0; i < num_conn; i++ )
1010  {
1011  coords[0] += dum_pos[3 * i];
1012  coords[1] += dum_pos[3 * i + 1];
1013  coords[2] += dum_pos[3 * i + 2];
1014  }
1015  coords[0] /= num_conn;
1016  coords[1] /= num_conn;
1017  coords[2] /= num_conn;
1018  }
1019  coords += 3;
1020  ++iter;
1021  }
1022 
1023  return status;
1024 }
1025 
1027  const double*& x,
1028  const double*& y,
1029  const double*& z ) const
1030 {
1032 
1033  if( TYPE_FROM_HANDLE( entity_handle ) == MBVERTEX )
1034  {
1035  const EntitySequence* seq = 0;
1036  status = sequence_manager()->find( entity_handle, seq );
1037 
1038  if( seq == 0 || status != MB_SUCCESS ) return MB_ENTITY_NOT_FOUND;
1039 
1040  status = static_cast< const VertexSequence* >( seq )->get_coordinates_ref( entity_handle, x, y, z );
1041  }
1042 
1043  return status;
1044 }
1045 
1046 //! set the coordinate information for this handle if it is of type Vertex
1047 //! otherwise, return an error
1048 ErrorCode Core::set_coords( const EntityHandle* entity_handles, const int num_entities, const double* coords )
1049 {
1050 
1051  ErrorCode status = MB_SUCCESS;
1052 
1053  int i, j = 0;
1054 
1055  for( i = 0; i < num_entities; i++ )
1056  {
1057  if( TYPE_FROM_HANDLE( entity_handles[i] ) == MBVERTEX )
1058  {
1059  EntitySequence* seq = 0;
1060  status = sequence_manager()->find( entity_handles[i], seq );
1061 
1062  if( seq != 0 && status == MB_SUCCESS )
1063  {
1064  status = static_cast< VertexSequence* >( seq )->set_coordinates( entity_handles[i], coords[j],
1065  coords[j + 1], coords[j + 2] );
1066  j += 3;
1067  }
1068  }
1069  else if( status == MB_SUCCESS )
1070  status = MB_TYPE_OUT_OF_RANGE;
1071  }
1072 
1073  return status;
1074 }
1075 
1076 //! set the coordinate information for this handle if it is of type Vertex
1077 //! otherwise, return an error
1078 ErrorCode Core::set_coords( Range entity_handles, const double* coords )
1079 {
1080 
1081  ErrorCode status = MB_SUCCESS;
1082 
1083  int j = 0;
1084 
1085  for( Range::iterator rit = entity_handles.begin(); rit != entity_handles.end(); ++rit )
1086  {
1087  if( TYPE_FROM_HANDLE( *rit ) == MBVERTEX )
1088  {
1089  EntitySequence* seq = 0;
1090  status = sequence_manager()->find( *rit, seq );
1091 
1092  if( seq != 0 && status == MB_SUCCESS )
1093  {
1094  status = static_cast< VertexSequence* >( seq )->set_coordinates( *rit, coords[j], coords[j + 1],
1095  coords[j + 2] );
1096  j += 3;
1097  }
1098  }
1099  else if( status == MB_SUCCESS )
1100  status = MB_TYPE_OUT_OF_RANGE;
1101  }
1102 
1103  return status;
1104 }
1105 
1107 {
1108  return sequenceManager->get_sequence_multiplier();
1109 }
1110 
1111 void Core::set_sequence_multiplier( double factor )
1112 {
1113  assert( factor >= 1.0 );
1114  sequenceManager->set_sequence_multiplier( factor );
1115 }
1116 
1117 //! get global connectivity array for specified entity type
1118 /** Assumes just vertices, no higher order nodes
1119  */
1120 ErrorCode Core::get_connectivity_by_type( const EntityType entity_type, std::vector< EntityHandle >& connect ) const
1121 {
1122  // inefficient implementation until we get blocked tag access
1123 
1124  // get the range of entities of this type
1125  Range this_range;
1126  ErrorCode result = get_entities_by_type( 0, entity_type, this_range );
1127 
1128  int num_ents = this_range.size();
1129  connect.reserve( num_ents * CN::VerticesPerEntity( entity_type ) );
1130 
1131  // now loop over these entities, getting connectivity for each
1132  for( Range::iterator this_it = this_range.begin(); this_it != this_range.end(); ++this_it )
1133  {
1134  const EntityHandle* connect_vec = NULL;
1135  result = get_connectivity( *this_it, connect_vec, num_ents, true );MB_CHK_ERR( result );
1136  connect.insert( connect.end(), &connect_vec[0], &connect_vec[num_ents] );
1137  }
1138 
1139  return MB_SUCCESS;
1140 }
1141 
1142 //! get the connectivity for element /handles. For non-element handles, return an error
1144  const int num_handles,
1145  Range& connectivity,
1146  bool corners_only ) const
1147 {
1148  std::vector< EntityHandle > tmp_connect;
1149  ErrorCode result = get_connectivity( entity_handles, num_handles, tmp_connect, corners_only );MB_CHK_ERR( result );
1150 
1151  std::sort( tmp_connect.begin(), tmp_connect.end() );
1152  std::copy( tmp_connect.rbegin(), tmp_connect.rend(), range_inserter( connectivity ) );
1153  return result;
1154 }
1155 
1156 //! get the connectivity for element /handles. For non-element handles, return an error
1158  const int num_handles,
1159  std::vector< EntityHandle >& connectivity,
1160  bool corners_only,
1161  std::vector< int >* offsets ) const
1162 {
1163  connectivity.clear(); // this seems wrong as compared to other API functions,
1164  // but changing it breaks lost of code, so I'm leaving
1165  // it in. - j.kraftcheck 2009-11-06
1166 
1167  ErrorCode rval;
1168  std::vector< EntityHandle > tmp_storage; // used only for structured mesh
1169  const EntityHandle* conn;
1170  int len;
1171  if( offsets ) offsets->push_back( 0 );
1172  for( int i = 0; i < num_handles; ++i )
1173  {
1174  rval = get_connectivity( entity_handles[i], conn, len, corners_only, &tmp_storage );MB_CHK_ERR( rval );
1175  connectivity.insert( connectivity.end(), conn, conn + len );
1176  if( offsets ) offsets->push_back( connectivity.size() );
1177  }
1178  return MB_SUCCESS;
1179 }
1180 
1181 //! get the connectivity for element handles. For non-element handles, return an error
1183  const EntityHandle*& connectivity,
1184  int& number_nodes,
1185  bool corners_only,
1186  std::vector< EntityHandle >* storage ) const
1187 {
1188  ErrorCode status;
1189 
1190  // Make sure the entity should have a connectivity.
1191  EntityType entity_type = TYPE_FROM_HANDLE( entity_handle );
1192 
1193  // WARNING: This is very dependent on the ordering of the EntityType enum
1194  if( entity_type < MBVERTEX || entity_type >= MBENTITYSET )
1195  return MB_TYPE_OUT_OF_RANGE;
1196 
1197  else if( entity_type == MBVERTEX )
1198  {
1199  return MB_FAILURE;
1200  }
1201 
1202  const EntitySequence* seq = 0;
1203 
1204  // We know that connectivity is stored in an EntitySequence so jump straight
1205  // to the entity sequence
1206  status = sequence_manager()->find( entity_handle, seq );
1207  if( seq == 0 || status != MB_SUCCESS ) return MB_ENTITY_NOT_FOUND;
1208 
1209  return static_cast< const ElementSequence* >( seq )->get_connectivity( entity_handle, connectivity, number_nodes,
1210  corners_only, storage );
1211 }
1212 
1213 //! set the connectivity for element handles. For non-element handles, return an error
1214 ErrorCode Core::set_connectivity( const EntityHandle entity_handle, EntityHandle* connect, const int num_connect )
1215 {
1216  ErrorCode status = MB_FAILURE;
1217 
1218  // Make sure the entity should have a connectivity.
1219  // WARNING: This is very dependent on the ordering of the EntityType enum
1220  EntityType entity_type = TYPE_FROM_HANDLE( entity_handle );
1221 
1222  EntitySequence* seq = 0;
1223 
1224  if( entity_type < MBVERTEX || entity_type > MBENTITYSET ) return MB_TYPE_OUT_OF_RANGE;
1225 
1226  status = sequence_manager()->find( entity_handle, seq );
1227  if( seq == 0 || status != MB_SUCCESS ) return ( status != MB_SUCCESS ? status : MB_ENTITY_NOT_FOUND );
1228 
1229  const EntityHandle* old_conn;
1230  int len;
1231  status = static_cast< ElementSequence* >( seq )->get_connectivity( entity_handle, old_conn, len );MB_CHK_ERR( status );
1232 
1233  aEntityFactory->notify_change_connectivity( entity_handle, old_conn, connect, num_connect );
1234 
1235  status = static_cast< ElementSequence* >( seq )->set_connectivity( entity_handle, connect, num_connect );
1236  if( status != MB_SUCCESS )
1237  aEntityFactory->notify_change_connectivity( entity_handle, connect, old_conn, num_connect );
1238 
1239  return status;
1240 }
1241 
1242 template < typename ITER >
1244  ITER begin,
1245  ITER end,
1246  int to_dimension,
1247  bool create_if_missing,
1248  Range& adj_entities )
1249 {
1250  const size_t DEFAULT_MAX_BLOCKS_SIZE = 4000;
1251  const size_t MAX_OUTER_ITERATIONS = 100;
1252 
1253  std::vector< EntityHandle > temp_vec, storage;
1254  std::vector< EntityHandle >::const_iterator ti;
1255  ErrorCode result = MB_SUCCESS, tmp_result;
1256  ITER i = begin;
1257  Range::iterator ins;
1258  const EntityHandle* conn;
1259  int conn_len;
1260 
1261  // Just copy any vertices from the input range into the output
1262  size_t remaining = end - begin;
1263  assert( begin + remaining == end );
1264 
1265  // How many entities to work with at once? 2000 or so shouldn't require
1266  // too much memory, but don't iterate in outer loop more than a
1267  // 1000 times (make it bigger if many input entiites.)
1268  const size_t block_size = std::max( DEFAULT_MAX_BLOCKS_SIZE, remaining / MAX_OUTER_ITERATIONS );
1269  while( remaining > 0 )
1270  {
1271  const size_t count = remaining > block_size ? block_size : remaining;
1272  remaining -= count;
1273  temp_vec.clear();
1274  for( size_t j = 0; j < count; ++i, ++j )
1275  {
1276  if( CN::Dimension( TYPE_FROM_HANDLE( *i ) ) == to_dimension )
1277  {
1278  temp_vec.push_back( *i );
1279  }
1280  else if( to_dimension == 0 && TYPE_FROM_HANDLE( *i ) != MBPOLYHEDRON )
1281  {
1282  tmp_result = gMB->get_connectivity( *i, conn, conn_len, false, &storage );
1283  if( MB_SUCCESS != tmp_result )
1284  {
1285  result = tmp_result;
1286  continue;
1287  }
1288  temp_vec.insert( temp_vec.end(), conn, conn + conn_len );
1289  }
1290  else
1291  {
1292  tmp_result = gMB->a_entity_factory()->get_adjacencies( *i, to_dimension, create_if_missing, temp_vec );
1293  if( MB_SUCCESS != tmp_result )
1294  {
1295  result = tmp_result;
1296  continue;
1297  }
1298  }
1299  }
1300 
1301  std::sort( temp_vec.begin(), temp_vec.end() );
1302  ins = adj_entities.begin();
1303  ti = temp_vec.begin();
1304  while( ti != temp_vec.end() )
1305  {
1306  EntityHandle first = *ti;
1307  EntityHandle second = *ti;
1308  for( ++ti; ti != temp_vec.end() && ( *ti - second <= 1 ); ++ti )
1309  second = *ti;
1310  ins = adj_entities.insert( ins, first, second );
1311  }
1312  }
1313  return result;
1314 }
1315 
1316 template < typename ITER >
1318  ITER begin,
1319  ITER end,
1320  const int to_dimension,
1321  const bool create_if_missing,
1322  std::vector< EntityHandle >& adj_entities )
1323 {
1324  const size_t SORT_THRESHOLD = 200;
1325  std::vector< EntityHandle > temp_vec;
1326  std::vector< EntityHandle >::iterator adj_it, w_it;
1327  ErrorCode result = MB_SUCCESS;
1328 
1329  if( begin == end )
1330  {
1331  adj_entities.clear(); // intersection
1332  return MB_SUCCESS;
1333  }
1334 
1335  // First iteration is a special case if input list is empty.
1336  // Rather than returning nothing (intersecting with empty
1337  // input list), we begin with the adjacencies for the first entity.
1338  if( adj_entities.empty() )
1339  {
1340  EntityType entity_type = TYPE_FROM_HANDLE( *begin );
1341  if( to_dimension == CN::Dimension( entity_type ) )
1342  adj_entities.push_back( *begin );
1343  else if( to_dimension == 0 && entity_type != MBPOLYHEDRON )
1344  {
1345  result = mb->get_connectivity( &( *begin ), 1, adj_entities );MB_CHK_ERR( result );
1346  }
1347  else
1348  {
1349  result = mb->a_entity_factory()->get_adjacencies( *begin, to_dimension, create_if_missing, adj_entities );MB_CHK_ERR( result );
1350  }
1351  ++begin;
1352  }
1353 
1354  for( ITER from_it = begin; from_it != end; ++from_it )
1355  {
1356  // running results kept in adj_entities; clear temp_vec, which is working space
1357  temp_vec.clear();
1358 
1359  // get the next set of adjacencies
1360  EntityType entity_type = TYPE_FROM_HANDLE( *from_it );
1361  if( to_dimension == CN::Dimension( entity_type ) )
1362  temp_vec.push_back( *from_it );
1363  else if( to_dimension == 0 && entity_type != MBPOLYHEDRON )
1364  {
1365  result = mb->get_connectivity( &( *from_it ), 1, temp_vec );MB_CHK_ERR( result );
1366  }
1367  else
1368  {
1369  result = mb->a_entity_factory()->get_adjacencies( *from_it, to_dimension, create_if_missing, temp_vec );MB_CHK_ERR( result );
1370  }
1371 
1372  // otherwise intersect with the current set of results
1373  w_it = adj_it = adj_entities.begin();
1374  if( temp_vec.size() * adj_entities.size() < SORT_THRESHOLD )
1375  {
1376  for( ; adj_it != adj_entities.end(); ++adj_it )
1377  if( std::find( temp_vec.begin(), temp_vec.end(), *adj_it ) != temp_vec.end() )
1378  {
1379  *w_it = *adj_it;
1380  ++w_it;
1381  }
1382  }
1383  else
1384  {
1385  std::sort( temp_vec.begin(), temp_vec.end() );
1386  for( ; adj_it != adj_entities.end(); ++adj_it )
1387  if( std::binary_search( temp_vec.begin(), temp_vec.end(), *adj_it ) )
1388  {
1389  *w_it = *adj_it;
1390  ++w_it;
1391  }
1392  }
1393  adj_entities.erase( w_it, adj_entities.end() );
1394 
1395  // we're intersecting, so if there are no more results, we're done
1396  if( adj_entities.empty() ) break;
1397  }
1398 
1399  return MB_SUCCESS;
1400 }
1401 
1402 template < typename ITER >
1404  ITER begin,
1405  ITER end,
1406  const int to_dimension,
1407  const bool create_if_missing,
1408  Range& adj_entities )
1409 {
1410  std::vector< EntityHandle > results;
1411  ErrorCode rval = moab::get_adjacencies_intersection( mb, begin, end, to_dimension, create_if_missing, results );MB_CHK_ERR( rval );
1412 
1413  if( adj_entities.empty() )
1414  {
1415  std::copy( results.begin(), results.end(), range_inserter( adj_entities ) );
1416  return MB_SUCCESS;
1417  }
1418 
1419  Range::iterator it = adj_entities.begin();
1420  while( it != adj_entities.end() )
1421  {
1422  if( std::find( results.begin(), results.end(), *it ) == results.end() )
1423  it = adj_entities.erase( it );
1424  else
1425  ++it;
1426  }
1427  return MB_SUCCESS;
1428 }
1429 
1430 ///////////////////////////////////////////////////////////////////
1431 //////////////////////////////////////////
1432 #ifdef MOAB_HAVE_AHF
1433 
1434 template < typename ITER >
1435 static inline ErrorCode get_adjacencies_intersection_ahf( Core* mb,
1436  ITER begin,
1437  ITER end,
1438  const int to_dimension,
1439  std::vector< EntityHandle >& adj_entities )
1440 {
1441  const size_t SORT_THRESHOLD = 200;
1442  std::vector< EntityHandle > temp_vec;
1443  std::vector< EntityHandle >::iterator adj_it, w_it;
1444  ErrorCode result = MB_SUCCESS;
1445 
1446  if( begin == end )
1447  {
1448  adj_entities.clear(); // intersection
1449  return MB_SUCCESS;
1450  }
1451 
1452  // First iteration is a special case if input list is empty.
1453  // Rather than returning nothing (intersecting with empty
1454  // input list), we begin with the adjacencies for the first entity.
1455  if( adj_entities.empty() )
1456  {
1457  EntityType entity_type = TYPE_FROM_HANDLE( *begin );
1458 
1459  if( to_dimension == 0 && entity_type != MBPOLYHEDRON )
1460  result = mb->get_connectivity( &( *begin ), 1, adj_entities );
1461  else
1462  result = mb->a_half_facet_rep()->get_adjacencies( *begin, to_dimension, adj_entities );
1463  if( MB_SUCCESS != result ) return result;
1464  ++begin;
1465  }
1466 
1467  for( ITER from_it = begin; from_it != end; ++from_it )
1468  {
1469  // running results kept in adj_entities; clear temp_vec, which is working space
1470  temp_vec.clear();
1471 
1472  // get the next set of adjacencies
1473  EntityType entity_type = TYPE_FROM_HANDLE( *from_it );
1474  if( to_dimension == 0 && entity_type != MBPOLYHEDRON )
1475  result = mb->get_connectivity( &( *from_it ), 1, temp_vec );
1476  else
1477  result = mb->a_half_facet_rep()->get_adjacencies( *from_it, to_dimension, temp_vec );
1478  if( MB_SUCCESS != result ) return result;
1479 
1480  // otherwise intersect with the current set of results
1481  w_it = adj_it = adj_entities.begin();
1482  if( temp_vec.size() * adj_entities.size() < SORT_THRESHOLD )
1483  {
1484  for( ; adj_it != adj_entities.end(); ++adj_it )
1485  if( std::find( temp_vec.begin(), temp_vec.end(), *adj_it ) != temp_vec.end() )
1486  {
1487  *w_it = *adj_it;
1488  ++w_it;
1489  }
1490  }
1491  else
1492  {
1493  std::sort( temp_vec.begin(), temp_vec.end() );
1494  for( ; adj_it != adj_entities.end(); ++adj_it )
1495  if( std::binary_search( temp_vec.begin(), temp_vec.end(), *adj_it ) )
1496  {
1497  *w_it = *adj_it;
1498  ++w_it;
1499  }
1500  }
1501  adj_entities.erase( w_it, adj_entities.end() );
1502 
1503  // we're intersecting, so if there are no more results, we're done
1504  if( adj_entities.empty() ) break;
1505  }
1506 
1507  return MB_SUCCESS;
1508 }
1509 #endif
1510 
1511 ///////////////////////////////////////////
1512 
1514  const int num_entities,
1515  const int to_dimension,
1516  const bool create_if_missing,
1517  std::vector< EntityHandle >& adj_entities,
1518  const int operation_type )
1519 {
1520 
1521 #ifdef MOAB_HAVE_AHF
1522  bool can_handle = true;
1523 
1524  if( to_dimension == 4 )
1525  can_handle = false; // NOT SUPPORTED: meshsets
1526  else if( create_if_missing )
1527  can_handle = false; // NOT SUPPORTED: create_if_missing
1528 
1529  bool mixed = ahfRep->check_mixed_entity_type(); // NOT SUPPORTED: mixed entity types or
1530  // polygonal/hedrals types
1531  if( mixed ) can_handle = false;
1532 
1533  if( mesh_modified ) // NOT SUPPORTED: modified mesh
1534  can_handle = false;
1535 
1536  if( can_handle )
1537  {
1538  ErrorCode result;
1539  if( operation_type == Interface::INTERSECT )
1540  return get_adjacencies_intersection_ahf( this, from_entities, from_entities + num_entities, to_dimension,
1541  adj_entities );
1542  else if( operation_type != Interface::UNION )
1543  return MB_FAILURE;
1544 
1545  // do union
1546 
1547  std::vector< EntityHandle > tmp_storage;
1548  const EntityHandle* conn;
1549  int len;
1550  for( int i = 0; i < num_entities; ++i )
1551  {
1552  if( to_dimension == 0 && TYPE_FROM_HANDLE( from_entities[0] ) != MBPOLYHEDRON )
1553  {
1554  result = get_connectivity( from_entities[i], conn, len, false, &tmp_storage );
1555  adj_entities.insert( adj_entities.end(), conn, conn + len );
1556  if( MB_SUCCESS != result ) return result;
1557  }
1558  else
1559  {
1560  result = ahfRep->get_adjacencies( from_entities[i], to_dimension, adj_entities );
1561  if( MB_SUCCESS != result ) return result;
1562  }
1563  }
1564  std::sort( adj_entities.begin(), adj_entities.end() );
1565  adj_entities.erase( std::unique( adj_entities.begin(), adj_entities.end() ), adj_entities.end() );
1566  }
1567  else
1568  {
1569 
1570 #endif
1571 
1572  if( operation_type == Interface::INTERSECT )
1573  return get_adjacencies_intersection( this, from_entities, from_entities + num_entities, to_dimension,
1574  create_if_missing, adj_entities );
1575  else if( operation_type != Interface::UNION )
1576  return MB_FAILURE;
1577 
1578  // do union
1579  ErrorCode result;
1580  std::vector< EntityHandle > tmp_storage;
1581  const EntityHandle* conn;
1582  int len;
1583  for( int i = 0; i < num_entities; ++i )
1584  {
1585  if( to_dimension == 0 && TYPE_FROM_HANDLE( from_entities[0] ) != MBPOLYHEDRON )
1586  {
1587  result = get_connectivity( from_entities[i], conn, len, false, &tmp_storage );MB_CHK_ERR( result );
1588  adj_entities.insert( adj_entities.end(), conn, conn + len );
1589  }
1590  else
1591  {
1592  result =
1593  aEntityFactory->get_adjacencies( from_entities[i], to_dimension, create_if_missing, adj_entities );MB_CHK_ERR( result );
1594  }
1595  }
1596  std::sort( adj_entities.begin(), adj_entities.end() );
1597  adj_entities.erase( std::unique( adj_entities.begin(), adj_entities.end() ), adj_entities.end() );
1598 
1599  // return MB_SUCCESS;
1600 
1601 #ifdef MOAB_HAVE_AHF
1602  }
1603 #endif
1604 
1605  return MB_SUCCESS;
1606 }
1607 
1609  const int num_entities,
1610  const int to_dimension,
1611  const bool create_if_missing,
1612  Range& adj_entities,
1613  const int operation_type )
1614 {
1615  if( operation_type == Interface::INTERSECT )
1616  return get_adjacencies_intersection( this, from_entities, from_entities + num_entities, to_dimension,
1617  create_if_missing, adj_entities );
1618  else if( operation_type == Interface::UNION )
1619  return get_adjacencies_union( this, from_entities, from_entities + num_entities, to_dimension,
1620  create_if_missing, adj_entities );
1621  else
1622  return MB_FAILURE;
1623 }
1624 
1625 ErrorCode Core::get_connectivity( const Range& from_entities, Range& adj_entities, bool corners_only ) const
1626 {
1627  const size_t DEFAULT_MAX_BLOCKS_SIZE = 4000;
1628  const size_t MAX_OUTER_ITERATIONS = 100;
1629 
1630  std::vector< EntityHandle > temp_vec, storage;
1631  std::vector< EntityHandle >::const_iterator ti;
1632  ErrorCode result = MB_SUCCESS, tmp_result;
1633  Range::const_iterator i = from_entities.begin();
1634  Range::iterator ins;
1635  const EntityHandle* conn;
1636  int conn_len;
1637 
1638  // Just copy any vertices from the input range into the output
1639  size_t remaining = from_entities.size();
1640  for( ; i != from_entities.end() && TYPE_FROM_HANDLE( *i ) == MBVERTEX; ++i )
1641  --remaining;
1642  adj_entities.merge( from_entities.begin(), i );
1643 
1644  // How many entities to work with at once? 2000 or so shouldn't require
1645  // too much memory, but don't iterate in outer loop more than a
1646  // 1000 times (make it bigger if many input entiites.)
1647  const size_t block_size = std::max( DEFAULT_MAX_BLOCKS_SIZE, remaining / MAX_OUTER_ITERATIONS );
1648  while( remaining > 0 )
1649  {
1650  const size_t count = remaining > block_size ? block_size : remaining;
1651  remaining -= count;
1652  temp_vec.clear();
1653  for( size_t j = 0; j < count; ++i, ++j )
1654  {
1655  tmp_result = get_connectivity( *i, conn, conn_len, corners_only, &storage );
1656  if( MB_SUCCESS != tmp_result )
1657  {
1658  result = tmp_result;
1659  continue;
1660  }
1661 
1662  const size_t oldsize = temp_vec.size();
1663  temp_vec.resize( oldsize + conn_len );
1664  memcpy( &temp_vec[oldsize], conn, sizeof( EntityHandle ) * conn_len );
1665  }
1666 
1667  std::sort( temp_vec.begin(), temp_vec.end() );
1668  ins = adj_entities.begin();
1669  ti = temp_vec.begin();
1670  while( ti != temp_vec.end() )
1671  {
1672  EntityHandle first = *ti;
1673  EntityHandle second = *ti;
1674  for( ++ti; ti != temp_vec.end() && ( *ti - second <= 1 ); ++ti )
1675  second = *ti;
1676  ins = adj_entities.insert( ins, first, second );
1677  }
1678  }
1679  return result;
1680 }
1681 
1684  EntityHandle*& connect,
1685  int& verts_per_entity,
1686  int& count )
1687 {
1688  // Make sure the entity should have a connectivity.
1689  EntityType entity_type = TYPE_FROM_HANDLE( *iter );
1690 
1691  // WARNING: This is very dependent on the ordering of the EntityType enum
1692  if( entity_type <= MBVERTEX || entity_type >= MBENTITYSET ) return MB_TYPE_OUT_OF_RANGE;
1693 
1694  EntitySequence* seq = NULL;
1695 
1696  // We know that connectivity is stored in an EntitySequence so jump straight
1697  // to the entity sequence
1698  ErrorCode rval = sequence_manager()->find( *iter, seq );
1699  if( !seq || rval != MB_SUCCESS ) return MB_ENTITY_NOT_FOUND;
1700 
1701  ElementSequence* eseq = dynamic_cast< ElementSequence* >( seq );
1702  assert( eseq != NULL );
1703 
1704  connect = eseq->get_connectivity_array();
1705  if( !connect )
1706  {
1707  MB_SET_ERR( MB_FAILURE, "Couldn't find connectivity array for start handle" );
1708  }
1709 
1710  connect += eseq->nodes_per_element() * ( *iter - eseq->start_handle() );
1711 
1712  EntityHandle real_end = std::min( eseq->end_handle(), *( iter.end_of_block() ) );
1713  if( *end ) real_end = std::min( real_end, *end );
1714  count = real_end - *iter + 1;
1715 
1716  verts_per_entity = eseq->nodes_per_element();
1717 
1718  return MB_SUCCESS;
1719 }
1720 
1721 ErrorCode Core::get_vertices( const Range& from_entities, Range& vertices )
1722 {
1723  Range range;
1724  ErrorCode rval = get_connectivity( from_entities, range );MB_CHK_ERR( rval );
1725 
1726  // If input contained polyhedra, connectivity will contain faces.
1727  // Get vertices from faces.
1728  if( !range.all_of_dimension( 0 ) )
1729  {
1730  Range::iterator it = range.upper_bound( MBVERTEX );
1731  Range polygons;
1732  polygons.merge( it, range.end() );
1733  range.erase( it, range.end() );
1734  rval = get_connectivity( polygons, range );MB_CHK_ERR( rval );
1735  }
1736 
1737  if( vertices.empty() )
1738  vertices.swap( range );
1739  else
1740  vertices.merge( range );
1741  return MB_SUCCESS;
1742 }
1743 
1744 ErrorCode Core::get_adjacencies( const Range& from_entities,
1745  const int to_dimension,
1746  const bool create_if_missing,
1747  Range& adj_entities,
1748  const int operation_type )
1749 {
1750  if( operation_type == Interface::INTERSECT )
1751  return get_adjacencies_intersection( this, from_entities.begin(), from_entities.end(), to_dimension,
1752  create_if_missing, adj_entities );
1753  else if( operation_type != Interface::UNION )
1754  return MB_FAILURE;
1755  else if( to_dimension == 0 )
1756  return get_vertices( from_entities, adj_entities );
1757  else
1758  return get_adjacencies_union( this, from_entities.begin(), from_entities.end(), to_dimension, create_if_missing,
1759  adj_entities );
1760 }
1761 
1762 ErrorCode Core::add_adjacencies( const EntityHandle entity_handle,
1763  const EntityHandle* adjacencies,
1764  const int num_handles,
1765  bool both_ways )
1766 {
1767  ErrorCode result = MB_SUCCESS;
1768 
1769  for( const EntityHandle* it = adjacencies; it != adjacencies + num_handles; it++ )
1770  {
1771  result = aEntityFactory->add_adjacency( entity_handle, *it, both_ways );MB_CHK_ERR( result );
1772  }
1773 
1774  return MB_SUCCESS;
1775 }
1776 
1777 ErrorCode Core::add_adjacencies( const EntityHandle entity_handle, Range& adjacencies, bool both_ways )
1778 {
1779  ErrorCode result = MB_SUCCESS;
1780 
1781  for( Range::iterator rit = adjacencies.begin(); rit != adjacencies.end(); ++rit )
1782  {
1783  result = aEntityFactory->add_adjacency( entity_handle, *rit, both_ways );MB_CHK_ERR( result );
1784  }
1785 
1786  return MB_SUCCESS;
1787 }
1788 
1789 ErrorCode Core::remove_adjacencies( const EntityHandle entity_handle,
1790  const EntityHandle* adjacencies,
1791  const int num_handles )
1792 {
1793  ErrorCode result = MB_SUCCESS;
1794 
1795  for( const EntityHandle* it = adjacencies; it != adjacencies + num_handles; it++ )
1796  {
1797  result = aEntityFactory->remove_adjacency( entity_handle, *it );MB_CHK_ERR( result );
1798  result = aEntityFactory->remove_adjacency( *it, entity_handle );MB_CHK_ERR( result );
1799  }
1800 
1801  return MB_SUCCESS;
1802 }
1803 
1806  const std::vector< EntityHandle >**& adjs_ptr,
1807  int& count )
1808 {
1809  // Make sure the entity should have a connectivity.
1810  EntityType entity_type = TYPE_FROM_HANDLE( *iter );
1811 
1812  // WARNING: This is very dependent on the ordering of the EntityType enum
1813  if( entity_type < MBVERTEX || entity_type > MBENTITYSET ) return MB_TYPE_OUT_OF_RANGE;
1814 
1815  EntitySequence* seq = NULL;
1816 
1817  // We know that connectivity is stored in an EntitySequence so jump straight
1818  // to the entity sequence
1819  ErrorCode rval = sequence_manager()->find( *iter, seq );
1820  if( !seq || rval != MB_SUCCESS ) return MB_ENTITY_NOT_FOUND;
1821 
1822  adjs_ptr = const_cast< const std::vector< EntityHandle >** >( seq->data()->get_adjacency_data() );
1823  if( !adjs_ptr ) return rval;
1824 
1825  adjs_ptr += *iter - seq->data()->start_handle();
1826 
1827  EntityHandle real_end = *( iter.end_of_block() );
1828  if( *end ) real_end = std::min( real_end, *end );
1829  count = real_end - *iter + 1;
1830 
1831  return MB_SUCCESS;
1832 }
1833 
1835  const int dimension,
1836  Range& entities,
1837  const bool recursive ) const
1838 {
1839  ErrorCode result = MB_SUCCESS;
1840  if( meshset )
1841  {
1842  const EntitySequence* seq;
1843  result = sequence_manager()->find( meshset, seq );MB_CHK_ERR( result );
1844  const MeshSetSequence* mseq = reinterpret_cast< const MeshSetSequence* >( seq );
1845  result = mseq->get_dimension( sequence_manager(), meshset, dimension, entities, recursive );MB_CHK_ERR( result );
1846  }
1847  else if( dimension > 3 )
1848  {
1849  sequence_manager()->get_entities( MBENTITYSET, entities );
1850  }
1851  else
1852  {
1853  for( EntityType this_type = CN::TypeDimensionMap[dimension].first;
1854  this_type <= CN::TypeDimensionMap[dimension].second; this_type++ )
1855  {
1856  sequence_manager()->get_entities( this_type, entities );
1857  }
1858  }
1859 
1860  return MB_SUCCESS;
1861 }
1862 
1864  const int dimension,
1865  std::vector< EntityHandle >& entities,
1866  const bool recursive ) const
1867 {
1868  ErrorCode result = MB_SUCCESS;
1869  if( meshset )
1870  {
1871  const EntitySequence* seq;
1872  result = sequence_manager()->find( meshset, seq );MB_CHK_ERR( result );
1873  const MeshSetSequence* mseq = reinterpret_cast< const MeshSetSequence* >( seq );
1874  result = mseq->get_dimension( sequence_manager(), meshset, dimension, entities, recursive );MB_CHK_ERR( result );
1875  }
1876  else if( dimension > 3 )
1877  {
1878  sequence_manager()->get_entities( MBENTITYSET, entities );
1879  }
1880  else
1881  {
1882  for( EntityType this_type = CN::TypeDimensionMap[dimension].first;
1883  this_type <= CN::TypeDimensionMap[dimension].second; this_type++ )
1884  {
1885  sequence_manager()->get_entities( this_type, entities );
1886  }
1887  }
1888 
1889  return MB_SUCCESS;
1890 }
1891 
1893  const EntityType entity_type,
1894  Range& entities,
1895  const bool recursive ) const
1896 {
1897  ErrorCode result = MB_SUCCESS;
1898  if( meshset )
1899  {
1900  const EntitySequence* seq;
1901  result = sequence_manager()->find( meshset, seq );MB_CHK_ERR( result );
1902  const MeshSetSequence* mseq = reinterpret_cast< const MeshSetSequence* >( seq );
1903  result = mseq->get_type( sequence_manager(), meshset, entity_type, entities, recursive );MB_CHK_ERR( result );
1904  }
1905  else
1906  {
1907  sequence_manager()->get_entities( entity_type, entities );
1908  }
1909 
1910  return MB_SUCCESS;
1911 }
1912 
1914  const EntityType entity_type,
1915  std::vector< EntityHandle >& entities,
1916  const bool recursive ) const
1917 {
1918  ErrorCode result = MB_SUCCESS;
1919  if( meshset )
1920  {
1921  const EntitySequence* seq;
1922  result = sequence_manager()->find( meshset, seq );MB_CHK_ERR( result );
1923  const MeshSetSequence* mseq = reinterpret_cast< const MeshSetSequence* >( seq );
1924  result = mseq->get_type( sequence_manager(), meshset, entity_type, entities, recursive );MB_CHK_ERR( result );
1925  }
1926  else
1927  {
1928  sequence_manager()->get_entities( entity_type, entities );
1929  }
1930 
1931  return MB_SUCCESS;
1932 }
1933 
1935  const EntityType entity_type,
1936  const Tag* tags,
1937  const void* const* values,
1938  const int num_tags,
1939  Range& entities,
1940  const int condition,
1941  const bool recursive ) const
1942 {
1943  ErrorCode result;
1944  Range range;
1945 
1946  result = get_entities_by_type( meshset, entity_type, range, recursive );MB_CHK_ERR( result );
1947  if( !entities.empty() && Interface::INTERSECT == condition ) range = intersect( entities, range );
1948 
1949  // For each tag:
1950  // if operation is INTERSECT remove from 'range' any non-tagged entities
1951  // if operation is UNION add to 'entities' any tagged entities
1952  for( int it = 0; it < num_tags && !range.empty(); it++ )
1953  {
1954  if( !valid_tag_handle( tags[it] ) ) return MB_TAG_NOT_FOUND;
1955 
1956  // Of the entities in 'range', put in 'tmp_range' the subset
1957  // that are tagged as requested for this tag.
1958  Range tmp_range;
1959 
1960  // get the entities with this tag/value combo
1961  if( NULL == values || NULL == values[it] )
1962  {
1963  result = tags[it]->get_tagged_entities( sequenceManager, tmp_range, entity_type, &range );MB_CHK_ERR( result );
1964  }
1965  else
1966  {
1967  result = tags[it]->find_entities_with_value( sequenceManager, mError, tmp_range, values[it], 0, entity_type,
1968  &range );MB_CHK_ERR( result );
1969  // if there is a default value, then we should return all entities
1970  // that are untagged
1971  if( tags[it]->equals_default_value( values[it] ) )
1972  {
1973  Range all_tagged, untagged;
1974  result = tags[it]->get_tagged_entities( sequenceManager, all_tagged, entity_type, &range );MB_CHK_ERR( result );
1975  // add to 'tmp_range' any untagged entities in 'range'
1976  tmp_range.merge( subtract( range, all_tagged ) );
1977  }
1978  }
1979 
1980  // The above calls should have already done the intersect for us.
1981  if( Interface::INTERSECT == condition )
1982  range.swap( tmp_range );
1983  else
1984  entities.merge( tmp_range );
1985  }
1986 
1987  if( Interface::INTERSECT == condition ) entities.swap( range );
1988 
1989  return MB_SUCCESS;
1990 }
1991 
1992 ErrorCode Core::get_entities_by_handle( const EntityHandle meshset, Range& entities, const bool recursive ) const
1993 {
1994  ErrorCode result = MB_SUCCESS;
1995  if( meshset )
1996  {
1997  const EntitySequence* seq;
1998  result = sequence_manager()->find( meshset, seq );MB_CHK_ERR( result );
1999  const MeshSetSequence* mseq = reinterpret_cast< const MeshSetSequence* >( seq );
2000  result = mseq->get_entities( sequence_manager(), meshset, entities, recursive );MB_CHK_ERR( result );
2001  }
2002  else
2003  {
2004  // iterate backwards so range insertion is quicker
2005  for( EntityType entity_type = MBENTITYSET; entity_type >= MBVERTEX; --entity_type )
2006  sequence_manager()->get_entities( entity_type, entities );
2007  }
2008 
2009  return MB_SUCCESS;
2010 }
2011 
2013  std::vector< EntityHandle >& entities,
2014  const bool recursive ) const
2015 {
2016  ErrorCode result;
2017  if( recursive || !meshset )
2018  {
2019  Range tmp_range;
2020  result = get_entities_by_handle( meshset, tmp_range, recursive );
2021  size_t offset = entities.size();
2022  entities.resize( offset + tmp_range.size() );
2023  std::copy( tmp_range.begin(), tmp_range.end(), entities.begin() + offset );
2024  }
2025  else
2026  {
2027  const EntitySequence* seq;
2028  result = sequence_manager()->find( meshset, seq );MB_CHK_ERR( result );
2029  const MeshSetSequence* mseq = reinterpret_cast< const MeshSetSequence* >( seq );
2030  result = mseq->get_entities( meshset, entities );MB_CHK_ERR( result );
2031  }
2032  return MB_SUCCESS;
2033 }
2034 
2035 //! get # entities of a given dimension
2037  const int dim,
2038  int& number,
2039  const bool recursive ) const
2040 {
2041  ErrorCode result = MB_SUCCESS;
2042 
2043  if( !meshset )
2044  {
2045  number = 0;
2046  for( EntityType this_type = CN::TypeDimensionMap[dim].first; this_type <= CN::TypeDimensionMap[dim].second;
2047  this_type++ )
2048  {
2049  number += sequence_manager()->get_number_entities( this_type );
2050  }
2051  }
2052  else
2053  {
2054  const EntitySequence* seq;
2055  result = sequence_manager()->find( meshset, seq );MB_CHK_ERR( result );
2056  const MeshSetSequence* mseq = reinterpret_cast< const MeshSetSequence* >( seq );
2057  result = mseq->num_dimension( sequence_manager(), meshset, dim, number, recursive );MB_CHK_ERR( result );
2058  }
2059 
2060  return MB_SUCCESS;
2061 }
2062 
2063 //! returns the number of entities with a given type and tag
2065  const EntityType entity_type,
2066  int& num_ent,
2067  const bool recursive ) const
2068 {
2069  ErrorCode result = MB_SUCCESS;
2070 
2071  if( recursive && entity_type == MBENTITYSET ) // will never return anything
2072  return MB_TYPE_OUT_OF_RANGE;
2073 
2074  if( meshset )
2075  {
2076  const EntitySequence* seq;
2077  result = sequence_manager()->find( meshset, seq );MB_CHK_ERR( result );
2078  const MeshSetSequence* mseq = reinterpret_cast< const MeshSetSequence* >( seq );
2079  result = mseq->num_type( sequence_manager(), meshset, entity_type, num_ent, recursive );MB_CHK_ERR( result );
2080  }
2081  else
2082  {
2083  num_ent = sequence_manager()->get_number_entities( entity_type );
2084  }
2085 
2086  return MB_SUCCESS;
2087 }
2088 
2090  const EntityType entity_type,
2091  const Tag* tag_handles,
2092  const void* const* values,
2093  const int num_tags,
2094  int& num_entities,
2095  int condition,
2096  const bool recursive ) const
2097 {
2098  Range dum_ents;
2099  ErrorCode result = get_entities_by_type_and_tag( meshset, entity_type, tag_handles, values, num_tags, dum_ents,
2100  condition, recursive );
2101  num_entities = dum_ents.size();
2102  return result;
2103 }
2104 
2105 ErrorCode Core::get_number_entities_by_handle( const EntityHandle meshset, int& num_ent, const bool recursive ) const
2106 {
2107  ErrorCode result;
2108  if( meshset )
2109  {
2110  const EntitySequence* seq;
2111  result = sequence_manager()->find( meshset, seq );MB_CHK_ERR( result );
2112  const MeshSetSequence* mseq = reinterpret_cast< const MeshSetSequence* >( seq );
2113  return mseq->num_entities( sequence_manager(), meshset, num_ent, recursive );
2114  }
2115 
2116  num_ent = 0;
2117  for( EntityType this_type = MBVERTEX; this_type < MBMAXTYPE; this_type++ )
2118  {
2119  int dummy = 0;
2120  result = get_number_entities_by_type( 0, this_type, dummy );
2121  if( result != MB_SUCCESS )
2122  {
2123  num_ent = 0;
2124  return result;
2125  }
2126  num_ent += dummy;
2127  }
2128 
2129  return MB_SUCCESS;
2130 }
2131 
2132 //! return the tag data for a given EntityHandle and Tag
2134  const EntityHandle* entity_handles,
2135  int num_entities,
2136  void* tag_data ) const
2137 {
2138  assert( valid_tag_handle( tag_handle ) );
2140  return tag_handle->get_data( sequenceManager, mError, entity_handles, num_entities, tag_data );
2141 }
2142 
2143 //! return the tag data for a given EntityHandle and Tag
2144 ErrorCode Core::tag_get_data( const Tag tag_handle, const Range& entity_handles, void* tag_data ) const
2145 {
2146  assert( valid_tag_handle( tag_handle ) );
2147  return tag_handle->get_data( sequenceManager, mError, entity_handles, tag_data );
2148 }
2149 
2150 //! set the data for given EntityHandles and Tag
2152  const EntityHandle* entity_handles,
2153  int num_entities,
2154  const void* tag_data )
2155 {
2156  assert( valid_tag_handle( tag_handle ) );
2158  return tag_handle->set_data( sequenceManager, mError, entity_handles, num_entities, tag_data );
2159 }
2160 
2161 //! set the data for given EntityHandles and Tag
2162 ErrorCode Core::tag_set_data( Tag tag_handle, const Range& entity_handles, const void* tag_data )
2163 {
2164  assert( valid_tag_handle( tag_handle ) );
2165  return tag_handle->set_data( sequenceManager, mError, entity_handles, tag_data );
2166 }
2167 
2168 //! return the tag data for a given EntityHandle and Tag
2170  const EntityHandle* entity_handles,
2171  int num_entities,
2172  const void** tag_data,
2173  int* tag_sizes ) const
2174 {
2175  assert( valid_tag_handle( tag_handle ) );
2177  ErrorCode result =
2178  tag_handle->get_data( sequenceManager, mError, entity_handles, num_entities, tag_data, tag_sizes );
2179  int typesize = TagInfo::size_from_data_type( tag_handle->get_data_type() );
2180  if( tag_sizes && typesize != 1 )
2181  for( int i = 0; i < num_entities; ++i )
2182  tag_sizes[i] /= typesize;
2183  return result;
2184 }
2185 
2186 //! return the tag data for a given EntityHandle and Tag
2188  const Range& entity_handles,
2189  const void** tag_data,
2190  int* tag_sizes ) const
2191 {
2192  assert( valid_tag_handle( tag_handle ) );
2193  ErrorCode result = tag_handle->get_data( sequenceManager, mError, entity_handles, tag_data, tag_sizes );
2194  int typesize = TagInfo::size_from_data_type( tag_handle->get_data_type() );
2195  if( tag_sizes && typesize != 1 )
2196  {
2197  int num_entities = entity_handles.size();
2198  for( int i = 0; i < num_entities; ++i )
2199  tag_sizes[i] /= typesize;
2200  }
2201  return result;
2202 }
2203 
2204 //! set the data for given EntityHandles and Tag
2206  const EntityHandle* entity_handles,
2207  int num_entities,
2208  void const* const* tag_data,
2209  const int* tag_sizes )
2210 {
2211  assert( valid_tag_handle( tag_handle ) );
2213  std::vector< int > tmp_sizes;
2214  int typesize = TagInfo::size_from_data_type( tag_handle->get_data_type() );
2215  if( typesize != 1 && tag_sizes )
2216  {
2217  tmp_sizes.resize( num_entities );
2218  for( int i = 0; i < num_entities; ++i )
2219  tmp_sizes[i] = tag_sizes[i] * typesize;
2220  tag_sizes = &tmp_sizes[0];
2221  }
2222  return tag_handle->set_data( sequenceManager, mError, entity_handles, num_entities, tag_data, tag_sizes );
2223 }
2224 
2225 //! set the data for given EntityHandles and Tag
2227  const Range& entity_handles,
2228  void const* const* tag_data,
2229  const int* tag_sizes )
2230 {
2231  assert( valid_tag_handle( tag_handle ) );
2232  std::vector< int > tmp_sizes;
2233  int typesize = TagInfo::size_from_data_type( tag_handle->get_data_type() );
2234  if( typesize != 1 && tag_sizes )
2235  {
2236  int num_entities = entity_handles.size();
2237  tmp_sizes.resize( num_entities );
2238  for( int i = 0; i < num_entities; ++i )
2239  tmp_sizes[i] = tag_sizes[i] * typesize;
2240  tag_sizes = &tmp_sizes[0];
2241  }
2242  return tag_handle->set_data( sequenceManager, mError, entity_handles, tag_data, tag_sizes );
2243 }
2244 
2245 //! set the data for given EntityHandles and Tag
2247  const EntityHandle* entity_handles,
2248  int num_entities,
2249  const void* tag_data,
2250  int tag_size )
2251 {
2252  assert( valid_tag_handle( tag_handle ) );
2254  return tag_handle->clear_data( sequenceManager, mError, entity_handles, num_entities, tag_data,
2255  tag_size * TagInfo::size_from_data_type( tag_handle->get_data_type() ) );
2256 }
2257 
2258 //! set the data for given EntityHandles and Tag
2259 ErrorCode Core::tag_clear_data( Tag tag_handle, const Range& entity_handles, const void* tag_data, int tag_size )
2260 {
2261  assert( valid_tag_handle( tag_handle ) );
2262  return tag_handle->clear_data( sequenceManager, mError, entity_handles, tag_data,
2263  tag_size * TagInfo::size_from_data_type( tag_handle->get_data_type() ) );
2264 }
2265 
2266 static bool is_zero_bytes( const void* mem, size_t size )
2267 {
2268  const char* iter = reinterpret_cast< const char* >( mem );
2269  const char* const end = iter + size;
2270  for( ; iter != end; ++iter )
2271  if( *iter ) return false;
2272  return true;
2273 }
2274 
2276  int size,
2277  DataType data_type,
2278  Tag& tag_handle,
2279  unsigned flags,
2280  const void* default_value,
2281  bool* created )
2282 {
2283  if( created ) *created = false;
2284 
2285  // we always work with sizes in bytes internally
2286  if( !( ( flags & MB_TAG_VARLEN ) && size == MB_VARIABLE_LENGTH ) )
2287  {
2288  if( flags & MB_TAG_BYTES )
2289  {
2290  if( size % TagInfo::size_from_data_type( data_type ) ) return MB_INVALID_SIZE;
2291  }
2292  else
2293  {
2294  size *= TagInfo::size_from_data_type( data_type );
2295  }
2296  }
2297 
2298  const TagType storage = static_cast< TagType >( flags & 3 );
2299 
2300  // search for an existing tag
2301  tag_handle = 0;
2302  if( name && *name )
2303  { // not anonymous tag
2304  for( std::list< Tag >::iterator i = tagList.begin(); i != tagList.end(); ++i )
2305  {
2306  if( ( *i )->get_name() == name )
2307  {
2308  tag_handle = *i;
2309  break;
2310  }
2311  }
2312  }
2313 
2314  if( tag_handle )
2315  {
2316  if( flags & MB_TAG_EXCL ) return MB_ALREADY_ALLOCATED;
2317  // user asked that we not check anything
2318  if( flags & MB_TAG_ANY ) return MB_SUCCESS;
2319  // user asked that we also match storage types
2320  if( ( flags & MB_TAG_STORE ) && tag_handle->get_storage_type() != storage ) return MB_TYPE_OUT_OF_RANGE;
2321  // check if data type matches
2322  const DataType extype = tag_handle->get_data_type();
2323  if( extype != data_type )
2324  {
2325  if( flags & MB_TAG_NOOPQ )
2326  return MB_TYPE_OUT_OF_RANGE;
2327  else if( extype != MB_TYPE_OPAQUE && data_type != MB_TYPE_OPAQUE )
2328  return MB_TYPE_OUT_OF_RANGE;
2329  }
2330 
2331  // Require that the size value be zero or MB_VARIABLE_LENGTH
2332  // for variable length tags. The caller passing such a size
2333  // value is sufficient to indicate that the caller is aware
2334  // that it is requesting a variable-length tag, so no need
2335  // to also require/check the MB_TAG_VARLEN bit in the flags.
2336  if( tag_handle->variable_length() )
2337  {
2338  if( size != 0 && size != MB_VARIABLE_LENGTH && !( flags & MB_TAG_VARLEN ) ) return MB_INVALID_SIZE;
2339  }
2340  // But /do/ fail if MB_TAG_VARLEN flag is set and tag is
2341  // not variable length.
2342  else if( flags & MB_TAG_VARLEN )
2343  return MB_TYPE_OUT_OF_RANGE;
2344  // check size for fixed-length tag
2345  else if( tag_handle->get_size() != size )
2346  return MB_INVALID_SIZE;
2347 
2348  // If user passed a default value, check that it matches.
2349  // If user did not pass a default value, assume they're OK
2350  // with the existing one.
2351  // If tag does not have a default value but the user passed
2352  // one, allow it only if the tag is dense and the passed value
2353  // is all zero bytes because dense tags have an implicit default
2354  // value of all zeros in some cases.
2355  if( default_value && !( flags & MB_TAG_DFTOK ) &&
2356  !( tag_handle->equals_default_value( default_value, size ) ||
2357  ( !tag_handle->get_default_value() && tag_handle->get_storage_type() == MB_TAG_DENSE &&
2358  is_zero_bytes( default_value, size ) ) ) )
2359  return MB_ALREADY_ALLOCATED;
2360 
2361  return MB_SUCCESS;
2362  }
2363 
2364  // MB_TAG_EXCL implies MB_TAG_CREAT
2365  if( !( flags & ( MB_TAG_CREAT | MB_TAG_EXCL ) ) ) return MB_TAG_NOT_FOUND;
2366 
2367  // if a non-opaque non-bit type was specified, then the size
2368  // must be multiple of the size of the type
2369  if( ( !( flags & MB_TAG_VARLEN ) || default_value ) &&
2370  ( size <= 0 || ( size % TagInfo::size_from_data_type( data_type ) ) != 0 ) )
2371  return MB_INVALID_SIZE;
2372 
2373  // if MB_TYPE_BIT may be used only with MB_TAG_BIT
2374  // if (storage != MB_TAG_BIT && data_type == MB_TYPE_BIT)
2375  // return MB_INVALID_ARG;
2376  if( data_type == MB_TYPE_BIT ) flags &= ~(unsigned)( MB_TAG_DENSE | MB_TAG_SPARSE );
2377 
2378  // create the tag
2379  switch( flags & ( MB_TAG_DENSE | MB_TAG_SPARSE | MB_TAG_MESH | MB_TAG_VARLEN ) )
2380  {
2381  case MB_TAG_DENSE | MB_TAG_VARLEN:
2382  tag_handle = VarLenDenseTag::create_tag( sequenceManager, mError, name, data_type, default_value, size );
2383  break;
2384  case MB_TAG_DENSE:
2385  tag_handle = DenseTag::create_tag( sequenceManager, mError, name, size, data_type, default_value );
2386  break;
2388  tag_handle = new VarLenSparseTag( name, data_type, default_value, size );
2389  break;
2390  case MB_TAG_SPARSE:
2391  tag_handle = new SparseTag( name, size, data_type, default_value );
2392  break;
2393  case MB_TAG_MESH | MB_TAG_VARLEN:
2394  tag_handle = new MeshTag( name, MB_VARIABLE_LENGTH, data_type, default_value, size );
2395  break;
2396  case MB_TAG_MESH:
2397  tag_handle = new MeshTag( name, size, data_type, default_value, size );
2398  break;
2399  case MB_TAG_BIT:
2400  if( MB_TYPE_BIT != data_type && MB_TYPE_OPAQUE != data_type ) return MB_TYPE_OUT_OF_RANGE;
2401  tag_handle = BitTag::create_tag( name, size, default_value );
2402  break;
2403  default: // some illegal combination (multiple storage types, variable-length bit tag, etc.)
2404  return MB_TYPE_OUT_OF_RANGE;
2405  }
2406 
2407  if( !tag_handle ) return MB_INVALID_SIZE;
2408 
2409  if( created ) *created = true;
2410  tagList.push_back( tag_handle );
2411  return MB_SUCCESS;
2412 }
2413 
2415  int size,
2416  DataType data_type,
2417  Tag& tag_handle,
2418  unsigned flags,
2419  const void* default_value ) const
2420 {
2421  // If caller specified MB_TAG_EXCL, then we must fail because
2422  // this const function can never create a tag. We need to test
2423  // this here because the non-const version of this function
2424  // assumes MB_TAG_CREAT if MB_TAG_EXCL is specified.
2425  if( flags & MB_TAG_EXCL )
2426  {
2427  // anonymous tag?
2428  if( !name || !*name ) return MB_TAG_NOT_FOUND;
2429 
2430  // search for an existing tag
2431  tag_handle = 0;
2432  for( std::list< Tag >::const_iterator i = tagList.begin(); i != tagList.end(); ++i )
2433  {
2434  if( ( *i )->get_name() == name )
2435  {
2436  tag_handle = *i;
2437  return MB_ALREADY_ALLOCATED;
2438  }
2439  }
2440 
2441  return MB_TAG_NOT_FOUND;
2442  }
2443 
2444  return const_cast< Core* >( this )->tag_get_handle( name, size, data_type, tag_handle,
2445  flags & ~(unsigned)MB_TAG_CREAT, default_value );
2446 }
2447 
2448 //! removes the tag from the entity
2449 ErrorCode Core::tag_delete_data( Tag tag_handle, const EntityHandle* entity_handles, int num_entities )
2450 {
2451  assert( valid_tag_handle( tag_handle ) );
2453  return tag_handle->remove_data( sequenceManager, mError, entity_handles, num_entities );
2454 }
2455 
2456 //! removes the tag from the entity
2457 ErrorCode Core::tag_delete_data( Tag tag_handle, const Range& entity_handles )
2458 {
2459  assert( valid_tag_handle( tag_handle ) );
2460  return tag_handle->remove_data( sequenceManager, mError, entity_handles );
2461 }
2462 
2463 //! removes the tag from MB
2465 {
2466  std::list< TagInfo* >::iterator i = std::find( tagList.begin(), tagList.end(), tag_handle );
2467  if( i == tagList.end() ) return MB_TAG_NOT_FOUND;
2468 
2469  ErrorCode rval = tag_handle->release_all_data( sequenceManager, mError, true );MB_CHK_ERR( rval );
2470 
2471  tagList.erase( i );
2472  delete tag_handle;
2473  return MB_SUCCESS;
2474 }
2475 
2477  Range::const_iterator iter,
2479  int& count,
2480  void*& data_ptr,
2481  bool allocate )
2482 {
2483  Range::const_iterator init = iter;
2484  assert( valid_tag_handle( tag_handle ) );
2485  ErrorCode result = tag_handle->tag_iterate( sequenceManager, mError, iter, end, data_ptr, allocate );
2486  if( MB_SUCCESS == result ) count = iter - init;
2487  return result;
2488 }
2489 
2490 //! gets the tag name string for the tag_handle
2491 ErrorCode Core::tag_get_name( const Tag tag_handle, std::string& tag_name ) const
2492 {
2493  if( !valid_tag_handle( tag_handle ) ) return MB_TAG_NOT_FOUND;
2494  tag_name = tag_handle->get_name();
2495  return MB_SUCCESS;
2496 }
2497 
2498 ErrorCode Core::tag_get_handle( const char* tag_name, Tag& tag_handle ) const
2499 {
2500  return tag_get_handle( tag_name, 0, MB_TYPE_OPAQUE, tag_handle, MB_TAG_ANY );
2501 }
2502 
2503 //! get size of tag in bytes
2504 ErrorCode Core::tag_get_bytes( const Tag tag_handle, int& tag_size ) const
2505 {
2506  if( !valid_tag_handle( tag_handle ) ) return MB_TAG_NOT_FOUND;
2507 
2508  if( tag_handle->variable_length() )
2509  {
2510  tag_size = MB_VARIABLE_LENGTH;
2511  return MB_VARIABLE_DATA_LENGTH;
2512  }
2513  else if( tag_handle->get_storage_type() == MB_TAG_BIT )
2514  {
2515  tag_size = 1;
2516  return MB_SUCCESS;
2517  }
2518  else
2519  {
2520  tag_size = tag_handle->get_size();
2521  return MB_SUCCESS;
2522  }
2523 }
2524 
2525 //! get size of tag in $values
2526 ErrorCode Core::tag_get_length( const Tag tag_handle, int& tag_size ) const
2527 {
2528  if( !valid_tag_handle( tag_handle ) ) return MB_TAG_NOT_FOUND;
2529 
2530  if( tag_handle->variable_length() )
2531  {
2532  tag_size = MB_VARIABLE_LENGTH;
2533  return MB_VARIABLE_DATA_LENGTH;
2534  }
2535  else
2536  {
2537  tag_size = tag_handle->get_size() / TagInfo::size_from_data_type( tag_handle->get_data_type() );
2538  return MB_SUCCESS;
2539  }
2540 }
2541 
2542 ErrorCode Core::tag_get_data_type( const Tag handle, DataType& data_type ) const
2543 {
2544  if( !valid_tag_handle( handle ) ) return MB_TAG_NOT_FOUND;
2545 
2546  data_type = handle->get_data_type();
2547  return MB_SUCCESS;
2548 }
2549 
2550 //! get default value of the tag
2551 ErrorCode Core::tag_get_default_value( const Tag tag_handle, void* def_value ) const
2552 {
2553  if( !valid_tag_handle( tag_handle ) ) return MB_TAG_NOT_FOUND;
2554 
2555  if( tag_handle->variable_length() ) return MB_VARIABLE_DATA_LENGTH;
2556 
2557  if( !tag_handle->get_default_value() ) return MB_ENTITY_NOT_FOUND;
2558 
2559  memcpy( def_value, tag_handle->get_default_value(), tag_handle->get_default_value_size() );
2560  return MB_SUCCESS;
2561 }
2562 
2563 ErrorCode Core::tag_get_default_value( Tag tag, const void*& ptr, int& size ) const
2564 {
2565  if( !valid_tag_handle( tag ) ) return MB_ENTITY_NOT_FOUND;
2566 
2567  if( !tag->get_default_value() ) return MB_ENTITY_NOT_FOUND;
2568 
2569  ptr = tag->get_default_value();
2571  return MB_SUCCESS;
2572 }
2573 
2574 //! get type of tag (sparse, dense, etc.; 0 = dense, 1 = sparse, 2 = bit, 3 = static)
2575 ErrorCode Core::tag_get_type( const Tag tag_handle, TagType& tag_type ) const
2576 {
2577  assert( valid_tag_handle( tag_handle ) );
2578  tag_type = tag_handle->get_storage_type();
2579  return MB_SUCCESS;
2580 }
2581 
2582 //! get handles for all tags defined
2583 ErrorCode Core::tag_get_tags( std::vector< Tag >& tag_handles ) const
2584 {
2585  std::copy( tagList.begin(), tagList.end(), std::back_inserter( tag_handles ) );
2586  return MB_SUCCESS;
2587 }
2588 
2589 //! Get handles for all tags defined on this entity
2590 ErrorCode Core::tag_get_tags_on_entity( const EntityHandle entity, std::vector< Tag >& tag_handles ) const
2591 {
2592  for( std::list< TagInfo* >::const_iterator i = tagList.begin(); i != tagList.end(); ++i )
2593  if( ( *i )->is_tagged( sequenceManager, entity ) ) tag_handles.push_back( *i );
2594  return MB_SUCCESS;
2595 }
2596 
2598 {
2599  const int negone = -1;
2600  if( 0 == materialTag )
2601  tag_get_handle( MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, materialTag, MB_TAG_CREAT | MB_TAG_SPARSE, &negone );
2602  return materialTag;
2603 }
2604 
2606 {
2607  const int negone = -1;
2608  if( 0 == neumannBCTag )
2609  tag_get_handle( NEUMANN_SET_TAG_NAME, 1, MB_TYPE_INTEGER, neumannBCTag, MB_TAG_CREAT | MB_TAG_SPARSE, &negone );
2610  return neumannBCTag;
2611 }
2612 
2614 {
2615  const int negone = -1;
2616  if( 0 == dirichletBCTag )
2617  tag_get_handle( DIRICHLET_SET_TAG_NAME, 1, MB_TYPE_INTEGER, dirichletBCTag, MB_TAG_CREAT | MB_TAG_SPARSE,
2618  &negone );
2619  return dirichletBCTag;
2620 }
2621 
2623 {
2624  const int negone = -1;
2625  if( 0 == globalIdTag )
2626  tag_get_handle( GLOBAL_ID_TAG_NAME, 1, MB_TYPE_INTEGER, globalIdTag, MB_TAG_CREAT | MB_TAG_DENSE, &negone );
2627  return globalIdTag;
2628 }
2629 
2631 {
2632  const int negone = -1;
2633  if( 0 == geomDimensionTag )
2634  tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, geomDimensionTag, MB_TAG_CREAT | MB_TAG_SPARSE,
2635  &negone );
2636  return geomDimensionTag;
2637 }
2638 
2639 //! creates an element based on the type and connectivity. returns a handle and error code
2640 ErrorCode Core::create_element( const EntityType entity_type,
2641  const EntityHandle* connectivity,
2642  const int num_nodes,
2643  EntityHandle& handle )
2644 {
2645  // make sure we have enough vertices for this entity type
2646  if( num_nodes < CN::VerticesPerEntity( entity_type ) ) return MB_FAILURE;
2647 
2648  ErrorCode status = sequence_manager()->create_element( entity_type, connectivity, num_nodes, handle );
2649  if( MB_SUCCESS == status ) status = aEntityFactory->notify_create_entity( handle, connectivity, num_nodes );
2650 
2651 #ifdef MOAB_HAVE_AHF
2652  mesh_modified = true;
2653 #endif
2654 
2655  return status;
2656 }
2657 
2658 //! creates a vertex based on coordinates, returns a handle and error code
2659 ErrorCode Core::create_vertex( const double coords[3], EntityHandle& handle )
2660 {
2661  // get an available vertex handle
2662  return sequence_manager()->create_vertex( coords, handle );
2663 }
2664 
2665 ErrorCode Core::create_vertices( const double* coordinates, const int nverts, Range& entity_handles )
2666 {
2667  // Create vertices
2668  ReadUtilIface* read_iface;
2669  ErrorCode result = Interface::query_interface( read_iface );MB_CHK_ERR( result );
2670 
2671  std::vector< double* > arrays;
2672  EntityHandle start_handle_out = 0;
2673  result = read_iface->get_node_coords( 3, nverts, MB_START_ID, start_handle_out, arrays );
2674  Interface::release_interface( read_iface );MB_CHK_ERR( result );
2675  // Cppcheck warning (false positive): variable arrays is assigned a value that is never used
2676  for( int i = 0; i < nverts; i++ )
2677  {
2678  arrays[0][i] = coordinates[3 * i];
2679  arrays[1][i] = coordinates[3 * i + 1];
2680  arrays[2][i] = coordinates[3 * i + 2];
2681  }
2682 
2683  entity_handles.clear();
2684  entity_handles.insert( start_handle_out, start_handle_out + nverts - 1 );
2685 
2686  return MB_SUCCESS;
2687 }
2688 
2689 //! merges two entities
2691  EntityHandle entity_to_remove,
2692  bool auto_merge,
2693  bool delete_removed_entity )
2694 {
2695  if( auto_merge ) return MB_FAILURE;
2696 
2697  // The two entities to merge must not be the same entity.
2698  if( entity_to_keep == entity_to_remove ) return MB_FAILURE;
2699 
2700  // The two entities to merge must be of the same type
2701  EntityType type_to_keep = TYPE_FROM_HANDLE( entity_to_keep );
2702 
2703  if( type_to_keep != TYPE_FROM_HANDLE( entity_to_remove ) ) return MB_TYPE_OUT_OF_RANGE;
2704 
2705  // Make sure both entities exist before trying to merge.
2706  EntitySequence* seq = 0;
2707  ErrorCode result, status;
2708  status = sequence_manager()->find( entity_to_keep, seq );
2709  if( seq == 0 || status != MB_SUCCESS ) return MB_ENTITY_NOT_FOUND;
2710  status = sequence_manager()->find( entity_to_remove, seq );
2711  if( seq == 0 || status != MB_SUCCESS ) return MB_ENTITY_NOT_FOUND;
2712 
2713  // If auto_merge is not set, all sub-entities should
2714  // be merged if the entities are to be merged.
2715  int ent_dim = CN::Dimension( type_to_keep );
2716  if( ent_dim > 0 )
2717  {
2718  std::vector< EntityHandle > conn, conn2;
2719 
2720  result = get_connectivity( &entity_to_keep, 1, conn );MB_CHK_ERR( result );
2721  result = get_connectivity( &entity_to_remove, 1, conn2 );MB_CHK_ERR( result );
2722 
2723  // Check to see if we can merge before pulling adjacencies.
2724  int dum1, dum2;
2725  if( !auto_merge &&
2726  ( conn.size() != conn2.size() || !CN::ConnectivityMatch( &conn[0], &conn2[0], conn.size(), dum1, dum2 ) ) )
2727  return MB_FAILURE;
2728  }
2729 
2730  result = aEntityFactory->merge_adjust_adjacencies( entity_to_keep, entity_to_remove );
2731 
2732  if( MB_SUCCESS == result && delete_removed_entity ) result = delete_entities( &entity_to_remove, 1 );
2733 
2734  return result;
2735 }
2736 
2737 //! deletes an entity range
2739 {
2740  ErrorCode result = MB_SUCCESS, temp_result;
2741  Range failed_ents;
2742 
2743  for( std::list< TagInfo* >::iterator i = tagList.begin(); i != tagList.end(); ++i )
2744  {
2745  temp_result = ( *i )->remove_data( sequenceManager, mError, range );
2746  // ok if the error is tag_not_found, some ents may not have every tag on them
2747  if( MB_SUCCESS != temp_result && MB_TAG_NOT_FOUND != temp_result ) result = temp_result;
2748  }
2749 
2750  for( Range::const_reverse_iterator rit = range.rbegin(); rit != range.rend(); ++rit )
2751  {
2752 
2753  // tell AEntityFactory that this element is going away
2754  temp_result = aEntityFactory->notify_delete_entity( *rit );
2755  if( MB_SUCCESS != temp_result )
2756  {
2757  result = temp_result;
2758  failed_ents.insert( *rit );
2759  continue;
2760  }
2761 
2762  if( TYPE_FROM_HANDLE( *rit ) == MBENTITYSET )
2763  {
2764  if( MeshSet* ptr = get_mesh_set( sequence_manager(), *rit ) )
2765  {
2766  int j, count;
2767  const EntityHandle* rel;
2768  ptr->clear( *rit, a_entity_factory() );
2769  rel = ptr->get_parents( count );
2770  for( j = 0; j < count; ++j )
2771  remove_child_meshset( rel[j], *rit );
2772  rel = ptr->get_children( count );
2773  for( j = 0; j < count; ++j )
2774  remove_parent_meshset( rel[j], *rit );
2775  }
2776  }
2777  }
2778 
2779  if( !failed_ents.empty() )
2780  {
2781  Range dum_range = subtract( range, failed_ents );
2782  // don't test for success, since we'll return failure in this case
2783  sequence_manager()->delete_entities( mError, dum_range );
2784  }
2785  else
2786  // now delete the entities
2787  result = sequence_manager()->delete_entities( mError, range );
2788 
2789  return result;
2790 }
2791 
2792 //! deletes an entity vector
2793 ErrorCode Core::delete_entities( const EntityHandle* entities, const int num_entities )
2794 {
2795  ErrorCode result = MB_SUCCESS, temp_result;
2796  Range failed_ents;
2797 
2798  for( std::list< TagInfo* >::iterator i = tagList.begin(); i != tagList.end(); ++i )
2799  {
2800  temp_result = ( *i )->remove_data( sequenceManager, mError, entities, num_entities );
2801  // ok if the error is tag_not_found, some ents may not have every tag on them
2802  if( MB_SUCCESS != temp_result && MB_TAG_NOT_FOUND != temp_result ) result = temp_result;
2803  }
2804 
2805  for( int i = 0; i < num_entities; i++ )
2806  {
2807 
2808  // tell AEntityFactory that this element is going away
2809  bool failed = false;
2810  temp_result = aEntityFactory->notify_delete_entity( entities[i] );
2811  if( MB_SUCCESS != temp_result )
2812  {
2813  result = temp_result;
2814  failed = true;
2815  }
2816 
2817  if( TYPE_FROM_HANDLE( entities[i] ) == MBENTITYSET )
2818  {
2819  if( MeshSet* ptr = get_mesh_set( sequence_manager(), entities[i] ) )
2820  {
2821  int j, count;
2822  const EntityHandle* rel;
2823  ptr->clear( entities[i], a_entity_factory() );
2824  rel = ptr->get_parents( count );
2825  for( j = 0; j < count; ++j )
2826  remove_child_meshset( rel[j], entities[i] );
2827  rel = ptr->get_children( count );
2828  for( j = 0; j < count; ++j )
2829  remove_parent_meshset( rel[j], entities[i] );
2830  }
2831  }
2832 
2833  if( failed )
2834  // don't test for success, since we'll return failure in this case
2835  sequence_manager()->delete_entity( mError, entities[i] );
2836  else
2837  {
2838  // now delete the entity
2839  temp_result = sequence_manager()->delete_entity( mError, entities[i] );
2840  if( MB_SUCCESS != temp_result ) result = temp_result;
2841  }
2842  }
2843 
2844  return result;
2845 }
2846 
2847 ErrorCode Core::list_entities( const EntityHandle* entities, const int num_entities ) const
2848 {
2849  Range temp_range;
2850  ErrorCode result = MB_SUCCESS;
2851  if( NULL == entities && num_entities == 0 )
2852  {
2853  // just list the numbers of each entity type
2854  int num_ents;
2855  std::cout << std::endl;
2856  std::cout << "Number of entities per type: " << std::endl;
2857  for( EntityType this_type = MBVERTEX; this_type < MBMAXTYPE; this_type++ )
2858  {
2859  result = get_number_entities_by_type( 0, this_type, num_ents );
2860  std::cout << CN::EntityTypeName( this_type ) << ": " << num_ents << std::endl;
2861  }
2862  std::cout << std::endl;
2863 
2864  return MB_SUCCESS;
2865  }
2866  else if( NULL == entities && num_entities < 0 )
2867  {
2868 
2869  // list all entities of all types
2870  std::cout << std::endl;
2871  for( EntityType this_type = MBVERTEX; this_type < MBMAXTYPE; this_type++ )
2872  {
2873  result = get_entities_by_type( 0, this_type, temp_range );
2874  }
2875 
2876  return list_entities( temp_range );
2877  }
2878  else if( NULL == entities && num_entities > 0 )
2879  {
2880 
2881  // list all entities of type == num_entities
2882  std::cout << std::endl;
2883  result = get_entities_by_type( 0, (EntityType)num_entities, temp_range );
2884 
2885  return list_entities( temp_range );
2886  }
2887  else
2888  {
2889  ErrorCode tmp_result;
2890  for( int i = 0; i < num_entities; i++ )
2891  {
2892  EntityType this_type = TYPE_FROM_HANDLE( entities[i] );
2893  std::cout << CN::EntityTypeName( this_type ) << " " << ID_FROM_HANDLE( entities[i] ) << ":" << endl;
2894 
2895  tmp_result = ( const_cast< Core* >( this ) )->list_entity( entities[i] );
2896  if( MB_SUCCESS != tmp_result ) result = tmp_result;
2897  }
2898  }
2899 
2900  return result;
2901 }
2902 
2903 ErrorCode Core::list_entities( const Range& temp_range ) const
2904 {
2905  ErrorCode result = MB_SUCCESS, tmp_result;
2906 
2907  for( Range::const_iterator rit = temp_range.begin(); rit != temp_range.end(); ++rit )
2908  {
2909  EntityType this_type = TYPE_FROM_HANDLE( *rit );
2910  std::cout << CN::EntityTypeName( this_type ) << " " << ID_FROM_HANDLE( *rit ) << ":" << endl;
2911 
2912  tmp_result = ( const_cast< Core* >( this ) )->list_entity( *rit );
2913  if( MB_SUCCESS != tmp_result ) result = tmp_result;
2914  }
2915 
2916  return result;
2917 }
2918 
2920 {
2921  ErrorCode result;
2922  HandleVec adj_vec;
2923 
2924  if( !is_valid( entity ) )
2925  {
2926  std::cout << "(invalid)" << std::endl;
2927  return MB_SUCCESS;
2928  }
2929 
2930  if( 0 != globalIdTag )
2931  {
2932  int dum;
2933  result = tag_get_data( globalIdTag, &entity, 1, &dum );
2934  if( MB_SUCCESS == result ) std::cout << "Global id = " << dum << std::endl;
2935  }
2936 
2937  // list entity
2938  EntityType this_type = TYPE_FROM_HANDLE( entity );
2939  if( this_type == MBVERTEX )
2940  {
2941  double coords[3];
2942  result = get_coords( &( entity ), 1, coords );MB_CHK_ERR( result );
2943  std::cout << "Coordinates: (" << coords[0] << ", " << coords[1] << ", " << coords[2] << ")" << std::endl;
2944  }
2945  else if( this_type == MBENTITYSET )
2946  this->print( entity, "" );
2947 
2948  std::cout << " Adjacencies:" << std::endl;
2949  bool some = false;
2950  int multiple = 0;
2951  for( int dim = 0; dim <= 3; dim++ )
2952  {
2953  if( dim == CN::Dimension( this_type ) ) continue;
2954  adj_vec.clear();
2955  // use const_cast here 'cuz we're in a const function and we're passing 'false' for
2956  // create_if_missing, so we know we won't change anything
2957  result = ( const_cast< Core* >( this ) )->get_adjacencies( &( entity ), 1, dim, false, adj_vec );
2958  if( MB_FAILURE == result ) continue;
2959  for( HandleVec::iterator adj_it = adj_vec.begin(); adj_it != adj_vec.end(); ++adj_it )
2960  {
2961  if( adj_it != adj_vec.begin() )
2962  std::cout << ", ";
2963  else
2964  std::cout << " ";
2965  std::cout << CN::EntityTypeName( TYPE_FROM_HANDLE( *adj_it ) ) << " " << ID_FROM_HANDLE( *adj_it );
2966  }
2967  if( !adj_vec.empty() )
2968  {
2969  std::cout << std::endl;
2970  some = true;
2971  }
2972  if( MB_MULTIPLE_ENTITIES_FOUND == result ) multiple += dim;
2973  }
2974  if( !some ) std::cout << "(none)" << std::endl;
2975  const EntityHandle* explicit_adjs;
2976  int num_exp;
2977  aEntityFactory->get_adjacencies( entity, explicit_adjs, num_exp );
2978  if( NULL != explicit_adjs && 0 != num_exp )
2979  {
2980  std::cout << " Explicit adjacencies: ";
2981  for( int i = 0; i < num_exp; i++ )
2982  {
2983  if( i != 0 ) std::cout << ", ";
2984  std::cout << CN::EntityTypeName( TYPE_FROM_HANDLE( explicit_adjs[i] ) ) << " "
2985  << ID_FROM_HANDLE( explicit_adjs[i] );
2986  }
2987  std::cout << std::endl;
2988  }
2989  if( multiple != 0 ) std::cout << " (MULTIPLE = " << multiple << ")" << std::endl;
2990 
2991  result = print_entity_tags( std::string(), entity, MB_TAG_DENSE );
2992 
2993  std::cout << std::endl;
2994 
2995  return result;
2996 }
2997 
2999  const bool mid_side,
3000  const bool mid_face,
3001  const bool mid_volume,
3002  Interface::HONodeAddedRemoved* function_object )
3003 {
3004  HigherOrderFactory fact( this, function_object );
3005  return fact.convert( meshset, mid_side, mid_face, mid_volume );
3006 }
3007 
3008 //! function to get the side number given two elements; returns
3009 //! MB_FAILURE if child not related to parent; does *not* create adjacencies
3010 //! between parent and child
3012  const EntityHandle child,
3013  int& sd_number,
3014  int& sense,
3015  int& offset ) const
3016 {
3017  // get the connectivity of parent and child
3018  const EntityHandle *parent_conn = NULL, *child_conn = NULL;
3019  int num_parent_vertices = 0, num_child_vertices = 0;
3020  ErrorCode result = get_connectivity( parent, parent_conn, num_parent_vertices, true );
3021  if( MB_NOT_IMPLEMENTED == result )
3022  {
3023  static std::vector< EntityHandle > tmp_connect( CN::MAX_NODES_PER_ELEMENT );
3024  result = get_connectivity( parent, parent_conn, num_parent_vertices, true, &tmp_connect );
3025  }
3026  if( MB_SUCCESS != result ) return result;
3027 
3028  if( TYPE_FROM_HANDLE( child ) == MBVERTEX )
3029  {
3030  int child_index = std::find( parent_conn, parent_conn + num_parent_vertices, child ) - parent_conn;
3031  if( child_index == num_parent_vertices )
3032  {
3033  sd_number = -1;
3034  sense = 0;
3035  return MB_FAILURE;
3036  }
3037  else
3038  {
3039  sd_number = child_index;
3040  sense = 1;
3041  return MB_SUCCESS;
3042  }
3043  }
3044 
3045  if( TYPE_FROM_HANDLE( parent ) == MBPOLYHEDRON )
3046  {
3047  // find the child in the parent_conn connectivity list, and call it a day ..
3048  // it should work only for faces within a conn list
3049  for( int i = 0; i < num_parent_vertices; i++ )
3050  if( child == parent_conn[i] )
3051  {
3052  sd_number = i;
3053  sense = 1; // always
3054  offset = 0;
3055  return MB_SUCCESS;
3056  }
3057  return MB_FAILURE;
3058  }
3059  result = get_connectivity( child, child_conn, num_child_vertices, true );MB_CHK_ERR( result );
3060 
3061  // call handle vector-based function
3062  if( TYPE_FROM_HANDLE( parent ) != MBPOLYGON && TYPE_FROM_HANDLE( parent ) != MBPOLYHEDRON )
3063  {
3064 
3065  // find indices into parent_conn for each entry in child_conn
3066  int child_conn_indices[10];
3067  assert( (unsigned)num_child_vertices <= sizeof( child_conn_indices ) / sizeof( child_conn_indices[0] ) );
3068  for( int i = 0; i < num_child_vertices; ++i )
3069  {
3070  child_conn_indices[i] =
3071  std::find( parent_conn, parent_conn + num_parent_vertices, child_conn[i] ) - parent_conn;
3072  if( child_conn_indices[i] >= num_parent_vertices )
3073  {
3074  sd_number = -1;
3075  return MB_FAILURE;
3076  }
3077  }
3078 
3079  int temp_result = CN::SideNumber( TYPE_FROM_HANDLE( parent ), child_conn_indices, num_child_vertices,
3080  CN::Dimension( TYPE_FROM_HANDLE( child ) ), sd_number, sense, offset );
3081  return ( 0 == temp_result ? MB_SUCCESS : MB_FAILURE );
3082  }
3083  else if( TYPE_FROM_HANDLE( parent ) == MBPOLYGON )
3084  {
3085  // find location of 1st vertex; this works even for padded vertices
3086  const EntityHandle* first_v = std::find( parent_conn, parent_conn + num_parent_vertices, child_conn[0] );
3087  if( first_v == parent_conn + num_parent_vertices ) return MB_ENTITY_NOT_FOUND;
3088  sd_number = first_v - parent_conn;
3089  offset = sd_number;
3090  if( TYPE_FROM_HANDLE( child ) == MBVERTEX )
3091  {
3092  sense = 0;
3093  return MB_SUCCESS;
3094  }
3095  else if( TYPE_FROM_HANDLE( child ) == MBPOLYGON )
3096  {
3097  bool match = CN::ConnectivityMatch( parent_conn, child_conn, num_parent_vertices, sense, offset );
3098  sd_number = 0;
3099  if( match )
3100  return MB_SUCCESS;
3101  else
3102  return MB_ENTITY_NOT_FOUND;
3103  }
3104  else if( TYPE_FROM_HANDLE( child ) == MBEDGE )
3105  {
3106  // determine the actual number of vertices, for the padded case
3107  // the padded case could be like ABCDEFFF; num_parent_vertices=8,
3108  // actual_num_parent_vertices=6
3109  int actual_num_parent_vertices = num_parent_vertices;
3110  while( actual_num_parent_vertices >= 3 &&
3111  ( parent_conn[actual_num_parent_vertices - 2] == parent_conn[actual_num_parent_vertices - 1] ) )
3112  actual_num_parent_vertices--;
3113 
3114  if( parent_conn[( sd_number + 1 ) % num_parent_vertices] == child_conn[1] )
3115  sense = 1;
3116  else if( parent_conn[( sd_number + num_parent_vertices - 1 ) % num_parent_vertices] ==
3117  child_conn[1] ) // this will also cover edge AF for padded case, side will
3118  // be 0, sense -1
3119  sense = -1;
3120  // if edge FA in above example, we should return sd_number = 5, sense 1
3121  else if( ( sd_number == actual_num_parent_vertices - 1 ) && ( child_conn[1] == parent_conn[0] ) )
3122  sense = 1;
3123  else
3124  return MB_ENTITY_NOT_FOUND;
3125  return MB_SUCCESS;
3126  }
3127  }
3128 
3129  return MB_FAILURE;
3130 }
3131 
3132 //! given an entity and the connectivity and type of one of its subfacets, find the
3133 //! high order node on that subfacet, if any
3135  const EntityHandle* subfacet_conn,
3136  const EntityType subfacet_type,
3137  EntityHandle& hon ) const
3138 {
3139  hon = 0;
3140 
3141  EntityType parent_type = TYPE_FROM_HANDLE( parent_handle );
3142 
3143  // get the parent's connectivity
3144  const EntityHandle* parent_conn = NULL;
3145  int num_parent_vertices = 0;
3146  ErrorCode result = get_connectivity( parent_handle, parent_conn, num_parent_vertices, false );MB_CHK_ERR( result );
3147 
3148  // find whether this entity has ho nodes
3149  int mid_nodes[4];
3150  CN::HasMidNodes( parent_type, num_parent_vertices, mid_nodes );
3151 
3152  // check whether this entity has mid nodes on this dimension subfacet;
3153  // use dimension-1 because vertices don't have mid nodes
3154  if( !mid_nodes[CN::Dimension( subfacet_type )] ) return MB_SUCCESS;
3155 
3156  // ok, we have mid nodes; now must compute expected index in connectivity array;
3157  // ho nodes stored for edges, faces then entity
3158 
3159  // offset starts with # corner vertices
3160  int offset = CN::VerticesPerEntity( parent_type );
3161  int i;
3162 
3163  for( i = 0; i < CN::Dimension( subfacet_type ) - 1; i++ )
3164  // for each dimension lower than that of the subfacet we're looking for,
3165  // if this entity has midnodes in that dimension, increment offset by #
3166  // of subfacets of that dimension; use dimension-1 in loop because
3167  // canon numbering table only has 2 positions, for edges and faces;
3168  if( mid_nodes[i + 1] ) offset += CN::mConnectivityMap[parent_type][i].num_sub_elements;
3169 
3170  // now add the index of this subfacet; only need to if it's not the highest dimension
3171  if( subfacet_type != parent_type )
3172  {
3173 
3174  // find indices into parent_conn for each entry in child_conn
3175  unsigned subfacet_size = CN::VerticesPerEntity( subfacet_type );
3176  int subfacet_indices[10];
3177  assert( subfacet_size <= sizeof( subfacet_indices ) / sizeof( subfacet_indices[0] ) );
3178  for( unsigned j = 0; j < subfacet_size; ++j )
3179  {
3180  subfacet_indices[j] =
3181  std::find( parent_conn, parent_conn + num_parent_vertices, subfacet_conn[j] ) - parent_conn;
3182  if( subfacet_indices[j] >= num_parent_vertices )
3183  {
3184  return MB_FAILURE;
3185  }
3186  }
3187 
3188  int dum, side_no, temp_offset;
3189  int temp_result =
3190  CN::SideNumber( parent_type, subfacet_indices, subfacet_size, subfacet_type, side_no, dum, temp_offset );
3191  if( temp_result != 0 ) return MB_FAILURE;
3192 
3193  offset += side_no;
3194  }
3195 
3196  // offset shouldn't be off the end of the connectivity vector
3197  if( offset >= num_parent_vertices ) return MB_INDEX_OUT_OF_RANGE;
3198 
3199  hon = parent_conn[offset];
3200 
3201  return MB_SUCCESS;
3202 }
3203 
3204 //! given an entity and a target dimension & side number, get that entity
3206  const int dim,
3207  const int sd_number,
3208  EntityHandle& target_entity ) const
3209 {
3210  // get a handle on the connectivity
3211  const EntityHandle* verts;
3212  int num_verts;
3213  ErrorCode result = get_connectivity( source_entity, verts, num_verts );MB_CHK_ERR( result );
3214 
3215  // special case for vertices
3216  if( dim == 0 )
3217  {
3218  if( sd_number < num_verts )
3219  {
3220  target_entity = verts[sd_number];
3221  return MB_SUCCESS;
3222  }
3223 
3224  else
3225  return MB_INDEX_OUT_OF_RANGE;
3226  }
3227 
3228  // get the vertices comprising the target entity
3229  Range side_verts, target_ents;
3230  const EntityType source_type = TYPE_FROM_HANDLE( source_entity );
3231  // first get the indices
3232  std::vector< int > vertex_indices;
3233 
3234  int temp_result = CN::AdjacentSubEntities( source_type, &sd_number, 1, dim, 0, vertex_indices );
3235  if( 0 != temp_result ) return MB_FAILURE;
3236  // now get the actual vertices
3237  for( unsigned int i = 0; i < vertex_indices.size(); i++ )
3238  side_verts.insert( verts[vertex_indices[i]] );
3239 
3240  // now look for an entity of the correct type
3241  // use const_cast here 'cuz we're in a const function and we're passing 'false' for
3242  // create_if_missing, so we know we won't change anything
3243  result = ( const_cast< Core* >( this ) )->get_adjacencies( side_verts, dim, false, target_ents );
3244  if( MB_SUCCESS != result && MB_MULTIPLE_ENTITIES_FOUND != result ) return result;
3245 
3246  if( !target_ents.empty() && TYPE_FROM_HANDLE( *( target_ents.begin() ) ) != MBVERTEX &&
3247  TYPE_FROM_HANDLE( *( target_ents.begin() ) ) !=
3248  CN::mConnectivityMap[source_type][dim - 1].target_type[sd_number] )
3249  return MB_ENTITY_NOT_FOUND;
3250 
3251  if( !target_ents.empty() ) target_entity = *( target_ents.begin() );
3252 
3253  return result;
3254 }
3255 
3256 //-------------------------Set Functions---------------------//
3257 
3258 ErrorCode Core::create_meshset( const unsigned int setoptions, EntityHandle& ms_handle, int )
3259 {
3260  return sequence_manager()->create_mesh_set( setoptions, ms_handle );
3261 }
3262 
3263 ErrorCode Core::get_meshset_options( const EntityHandle ms_handle, unsigned int& setoptions ) const
3264 {
3265  if( !ms_handle )
3266  { // root set
3267  setoptions = MESHSET_SET | MESHSET_TRACK_OWNER;
3268  return MB_SUCCESS;
3269  }
3270 
3271  const MeshSet* set = get_mesh_set( sequence_manager(), ms_handle );
3272  if( !set ) return MB_ENTITY_NOT_FOUND;
3273 
3274  setoptions = set->flags();
3275  return MB_SUCCESS;
3276 }
3277 
3278 ErrorCode Core::set_meshset_options( const EntityHandle ms_handle, const unsigned int setoptions )
3279 {
3280  MeshSet* set = get_mesh_set( sequence_manager(), ms_handle );
3281  if( !set ) return MB_ENTITY_NOT_FOUND;
3282 
3283  return set->set_flags( setoptions, ms_handle, a_entity_factory() );
3284 }
3285 
3286 ErrorCode Core::clear_meshset( const EntityHandle* ms_handles, const int num_meshsets )
3287 {
3288  ErrorCode result = MB_SUCCESS;
3289  for( int i = 0; i < num_meshsets; ++i )
3290  {
3291  MeshSet* set = get_mesh_set( sequence_manager(), ms_handles[i] );
3292  if( set )
3293  set->clear( ms_handles[i], a_entity_factory() );
3294  else
3295  result = MB_ENTITY_NOT_FOUND;
3296  }
3297 
3298  return result;
3299 }
3300 
3302 {
3303  ErrorCode result = MB_SUCCESS;
3304  for( Range::iterator i = ms_handles.begin(); i != ms_handles.end(); ++i )
3305  {
3306  MeshSet* set = get_mesh_set( sequence_manager(), *i );
3307  if( set )
3308  set->clear( *i, a_entity_factory() );
3309  else
3310  result = MB_ENTITY_NOT_FOUND;
3311  }
3312 
3313  return result;
3314 }
3315 
3317 {
3318  MeshSet* set1 = get_mesh_set( sequence_manager(), meshset1 );
3319  MeshSet* set2 = get_mesh_set( sequence_manager(), meshset2 );
3320  if( !set1 || !set2 ) return MB_ENTITY_NOT_FOUND;
3321 
3322  return set1->subtract( set2, meshset1, a_entity_factory() );
3323 }
3324 
3326 {
3327  MeshSet* set1 = get_mesh_set( sequence_manager(), meshset1 );
3328  MeshSet* set2 = get_mesh_set( sequence_manager(), meshset2 );
3329  if( !set1 || !set2 ) return MB_ENTITY_NOT_FOUND;
3330 
3331  return set1->intersect( set2, meshset1, a_entity_factory() );
3332 }
3333 
3335 {
3336  MeshSet* set1 = get_mesh_set( sequence_manager(), meshset1 );
3337  MeshSet* set2 = get_mesh_set( sequence_manager(), meshset2 );
3338  if( !set1 || !set2 ) return MB_ENTITY_NOT_FOUND;
3339 
3340  return set1->unite( set2, meshset1, a_entity_factory() );
3341 }
3342 
3344 {
3345  MeshSet* set = get_mesh_set( sequence_manager(), meshset );
3346  if( set )
3347  return set->add_entities( entities, meshset, a_entity_factory() );
3348  else
3349  return MB_ENTITY_NOT_FOUND;
3350 }
3351 
3352 ErrorCode Core::add_entities( EntityHandle meshset, const EntityHandle* entities, const int num_entities )
3353 {
3354  MeshSet* set = get_mesh_set( sequence_manager(), meshset );
3355  if( set )
3356  return set->add_entities( entities, num_entities, meshset, a_entity_factory() );
3357  else
3358  return MB_ENTITY_NOT_FOUND;
3359 }
3360 
3361 //! remove a range of entities from a meshset
3363 {
3364  MeshSet* set = get_mesh_set( sequence_manager(), meshset );
3365  if( set )
3366  return set->remove_entities( entities, meshset, a_entity_factory() );
3367  else
3368  return MB_ENTITY_NOT_FOUND;
3369 }
3370 
3371 //! remove a vector of entities from a meshset
3372 ErrorCode Core::remove_entities( EntityHandle meshset, const EntityHandle* entities, const int num_entities )
3373 {
3374  MeshSet* set = get_mesh_set( sequence_manager(), meshset );
3375  if( set )
3376  return set->remove_entities( entities, num_entities, meshset, a_entity_factory() );
3377  else
3378  return MB_ENTITY_NOT_FOUND;
3379 }
3380 
3381 //! return true if all entities are contained in set
3383  const EntityHandle* entities,
3384  int num_entities,
3385  const int operation_type )
3386 {
3387  if( !meshset ) // root
3388  return true;
3389  else if( MeshSet* set = get_mesh_set( sequence_manager(), meshset ) )
3390  return set->contains_entities( entities, num_entities, operation_type );
3391  else
3392  return false;
3393 }
3394 
3395 // replace entities in a meshset
3397  const EntityHandle* old_entities,
3398  const EntityHandle* new_entities,
3399  int num_entities )
3400 {
3401  MeshSet* set = get_mesh_set( sequence_manager(), meshset );
3402  if( set )
3403  return set->replace_entities( meshset, old_entities, new_entities, num_entities, a_entity_factory() );
3404  else
3405  return MB_ENTITY_NOT_FOUND;
3406 }
3407 
3409  std::vector< EntityHandle >& parents,
3410  const int num_hops ) const
3411 {
3412  if( 0 == meshset ) return MB_ENTITY_NOT_FOUND;
3413 
3414  const EntitySequence* seq;
3415  ErrorCode rval = sequence_manager()->find( meshset, seq );
3416  if( MB_SUCCESS != rval ) return MB_ENTITY_NOT_FOUND;
3417  const MeshSetSequence* mseq = reinterpret_cast< const MeshSetSequence* >( seq );
3418 
3419  return mseq->get_parents( sequence_manager(), meshset, parents, num_hops );
3420 }
3421 
3422 ErrorCode Core::get_parent_meshsets( const EntityHandle meshset, Range& parents, const int num_hops ) const
3423 {
3424  if( 0 == meshset ) return MB_ENTITY_NOT_FOUND;
3425 
3426  std::vector< EntityHandle > parent_vec;
3427  ErrorCode result = get_parent_meshsets( meshset, parent_vec, num_hops );MB_CHK_ERR( result );
3428  std::sort( parent_vec.begin(), parent_vec.end() );
3429  std::copy( parent_vec.rbegin(), parent_vec.rend(), range_inserter( parents ) );
3430  return MB_SUCCESS;
3431 }
3432 
3434  std::vector< EntityHandle >& children,
3435  const int num_hops ) const
3436 {
3437  if( 0 == meshset ) return MB_ENTITY_NOT_FOUND;
3438 
3439  const EntitySequence* seq;
3440  ErrorCode rval = sequence_manager()->find( meshset, seq );
3441  if( MB_SUCCESS != rval ) return MB_ENTITY_NOT_FOUND;
3442  const MeshSetSequence* mseq = reinterpret_cast< const MeshSetSequence* >( seq );
3443 
3444  return mseq->get_children( sequence_manager(), meshset, children, num_hops );
3445 }
3446 
3447 ErrorCode Core::get_child_meshsets( const EntityHandle meshset, Range& children, const int num_hops ) const
3448 {
3449  if( 0 == meshset ) return MB_ENTITY_NOT_FOUND;
3450 
3451  std::vector< EntityHandle > child_vec;
3452  ErrorCode result = get_child_meshsets( meshset, child_vec, num_hops );MB_CHK_ERR( result );
3453  std::sort( child_vec.begin(), child_vec.end() );
3454  std::copy( child_vec.rbegin(), child_vec.rend(), range_inserter( children ) );
3455  return MB_SUCCESS;
3456 }
3457 
3459  std::vector< EntityHandle >& children,
3460  const int num_hops ) const
3461 {
3462  if( 0 == meshset )
3463  {
3464  return get_entities_by_type( meshset, MBENTITYSET, children );
3465  }
3466 
3467  const EntitySequence* seq;
3468  ErrorCode rval = sequence_manager()->find( meshset, seq );
3469  if( MB_SUCCESS != rval ) return MB_ENTITY_NOT_FOUND;
3470  const MeshSetSequence* mseq = reinterpret_cast< const MeshSetSequence* >( seq );
3471 
3472  return mseq->get_contained_sets( sequence_manager(), meshset, children, num_hops );
3473 }
3474 
3475 ErrorCode Core::get_contained_meshsets( const EntityHandle meshset, Range& children, const int num_hops ) const
3476 {
3477  if( 0 == meshset )
3478  {
3479  return get_entities_by_type( meshset, MBENTITYSET, children );
3480  }
3481 
3482  std::vector< EntityHandle > child_vec;
3483  ErrorCode result = get_contained_meshsets( meshset, child_vec, num_hops );MB_CHK_ERR( result );
3484  std::sort( child_vec.begin(), child_vec.end() );
3485  std::copy( child_vec.rbegin(), child_vec.rend(), range_inserter( children ) );
3486  return MB_SUCCESS;
3487 }
3488 
3489 ErrorCode Core::num_parent_meshsets( const EntityHandle meshset, int* number, const int num_hops ) const
3490 {
3491  if( 0 == meshset ) return MB_ENTITY_NOT_FOUND;
3492 
3493  const EntitySequence* seq;
3494  ErrorCode rval = sequence_manager()->find( meshset, seq );
3495  if( MB_SUCCESS != rval ) return MB_ENTITY_NOT_FOUND;
3496  const MeshSetSequence* mseq = reinterpret_cast< const MeshSetSequence* >( seq );
3497 
3498  return mseq->num_parents( sequence_manager(), meshset, *number, num_hops );
3499 }
3500 
3501 ErrorCode Core::num_child_meshsets( const EntityHandle meshset, int* number, const int num_hops ) const
3502 {
3503  if( 0 == meshset ) return MB_ENTITY_NOT_FOUND;
3504 
3505  const EntitySequence* seq;
3506  ErrorCode rval = sequence_manager()->find( meshset, seq );
3507  if( MB_SUCCESS != rval ) return MB_ENTITY_NOT_FOUND;
3508  const MeshSetSequence* mseq = reinterpret_cast< const MeshSetSequence* >( seq );
3509 
3510  return mseq->num_children( sequence_manager(), meshset, *number, num_hops );
3511 }
3512 
3513 ErrorCode Core::num_contained_meshsets( const EntityHandle meshset, int* number, const int num_hops ) const
3514 {
3515  if( 0 == meshset )
3516  {
3517  return get_number_entities_by_type( 0, MBENTITYSET, *number );
3518  }
3519 
3520  const EntitySequence* seq;
3521  ErrorCode rval = sequence_manager()->find( meshset, seq );
3522  if( MB_SUCCESS != rval ) return MB_ENTITY_NOT_FOUND;
3523  const MeshSetSequence* mseq = reinterpret_cast< const MeshSetSequence* >( seq );
3524 
3525  return mseq->num_contained_sets( sequence_manager(), meshset, *number, num_hops );
3526 }
3527 
3529 {
3530  MeshSet* set_ptr = get_mesh_set( sequence_manager(), meshset );
3531  MeshSet* parent_ptr = get_mesh_set( sequence_manager(), parent_meshset );
3532  if( !set_ptr || !parent_ptr ) return MB_ENTITY_NOT_FOUND;
3533 
3534  set_ptr->add_parent( parent_meshset );
3535  return MB_SUCCESS;
3536 }
3537 
3538 ErrorCode Core::add_parent_meshsets( EntityHandle meshset, const EntityHandle* parents, int count )
3539 {
3540  MeshSet* set_ptr = get_mesh_set( sequence_manager(), meshset );
3541  if( !set_ptr ) return MB_ENTITY_NOT_FOUND;
3542 
3543  for( int i = 0; i < count; ++i )
3544  if( !get_mesh_set( sequence_manager(), parents[i] ) ) return MB_ENTITY_NOT_FOUND;
3545 
3546  for( int i = 0; i < count; ++i )
3547  set_ptr->add_parent( parents[i] );
3548  return MB_SUCCESS;
3549 }
3550 
3552 {
3553  MeshSet* set_ptr = get_mesh_set( sequence_manager(), meshset );
3554  MeshSet* child_ptr = get_mesh_set( sequence_manager(), child_meshset );
3555  if( !set_ptr || !child_ptr ) return MB_ENTITY_NOT_FOUND;
3556 
3557  set_ptr->add_child( child_meshset );
3558  return MB_SUCCESS;
3559 }
3560 
3562 {
3563  MeshSet* set_ptr = get_mesh_set( sequence_manager(), meshset );
3564  if( !set_ptr ) return MB_ENTITY_NOT_FOUND;
3565 
3566  for( int i = 0; i < count; ++i )
3567  if( !get_mesh_set( sequence_manager(), children[i] ) ) return MB_ENTITY_NOT_FOUND;
3568 
3569  for( int i = 0; i < count; ++i )
3570  set_ptr->add_child( children[i] );
3571  return MB_SUCCESS;
3572 }
3573 
3575 {
3576  MeshSet* parent_ptr = get_mesh_set( sequence_manager(), parent );
3577  MeshSet* child_ptr = get_mesh_set( sequence_manager(), child );
3578  if( !parent_ptr || !child_ptr ) return MB_ENTITY_NOT_FOUND;
3579 
3580  parent_ptr->add_child( child );
3581  child_ptr->add_parent( parent );
3582  return MB_SUCCESS;
3583 }
3584 
3586 {
3587  MeshSet* parent_ptr = get_mesh_set( sequence_manager(), parent );
3588  MeshSet* child_ptr = get_mesh_set( sequence_manager(), child );
3589  if( !parent_ptr || !child_ptr ) return MB_ENTITY_NOT_FOUND;
3590 
3591  parent_ptr->remove_child( child );
3592  child_ptr->remove_parent( parent );
3593  return MB_SUCCESS;
3594 }
3595 
3597 {
3598  MeshSet* set_ptr = get_mesh_set( sequence_manager(), meshset );
3599  if( !set_ptr ) return MB_ENTITY_NOT_FOUND;
3600  set_ptr->remove_parent( parent_meshset );
3601  return MB_SUCCESS;
3602 }
3603 
3605 {
3606  MeshSet* set_ptr = get_mesh_set( sequence_manager(), meshset );
3607  if( !set_ptr ) return MB_ENTITY_NOT_FOUND;
3608  set_ptr->remove_child( child_meshset );
3609  return MB_SUCCESS;
3610 }
3611 
3612 ErrorCode Core::get_last_error( std::string& info ) const
3613 {
3615  return MB_SUCCESS;
3616 }
3617 
3618 std::string Core::get_error_string( const ErrorCode code ) const
3619 {
3620  return (unsigned)code <= (unsigned)MB_FAILURE ? ErrorCodeStr[code] : "INVALID ERROR CODE";
3621 }
3622 
3623 void Core::print( const EntityHandle ms_handle, const char* prefix, bool first_call ) const
3624 {
3625  // get the entities
3626  Range entities;
3627 
3628  if( 0 != ms_handle )
3629  {
3630  get_entities_by_handle( ms_handle, entities );
3631  std::cout << prefix << "MBENTITYSET " << ID_FROM_HANDLE( ms_handle ) << std::endl;
3632  }
3633  else
3634  {
3635  get_entities_by_dimension( 0, 3, entities );
3636  if( entities.empty() ) get_entities_by_dimension( 0, 2, entities );
3637  if( entities.empty() ) get_entities_by_dimension( 0, 1, entities );
3638  get_entities_by_dimension( 0, 0, entities );
3639  get_entities_by_type( 0, MBENTITYSET, entities );
3640  std::cout << prefix << "--: " << std::endl;
3641  }
3642 
3643  std::string indent_prefix = prefix;
3644  indent_prefix += " ";
3645  entities.print( indent_prefix.c_str() );
3646 
3647  if( !first_call || !ms_handle ) return;
3648 
3649  // print parent/children
3650  Range temp;
3651  this->get_parent_meshsets( ms_handle, temp );
3652  std::cout << " Parent sets: ";
3653  if( temp.empty() )
3654  std::cout << "(none)" << std::endl;
3655  else
3656  {
3657  for( Range::iterator rit = temp.begin(); rit != temp.end(); ++rit )
3658  {
3659  if( rit != temp.begin() ) std::cout << ", ";
3660  std::cout << ID_FROM_HANDLE( *rit );
3661  }
3662  std::cout << std::endl;
3663  }
3664 
3665  temp.clear();
3666  this->get_child_meshsets( ms_handle, temp );
3667  std::cout << " Child sets: ";
3668  if( temp.empty() )
3669  std::cout << "(none)" << std::endl;
3670  else
3671  {
3672  for( Range::iterator rit = temp.begin(); rit != temp.end(); ++rit )
3673  {
3674  if( rit != temp.begin() ) std::cout << ", ";
3675  std::cout << ID_FROM_HANDLE( *rit );
3676  }
3677  std::cout << std::endl;
3678  }
3679 
3680  // print all sparse tags
3681  print_entity_tags( indent_prefix, ms_handle, MB_TAG_SPARSE );
3682 }
3683 
3684 ErrorCode Core::print_entity_tags( std::string indent_prefix, const EntityHandle handle, TagType tp ) const
3685 {
3686  std::vector< Tag > set_tags;
3687  ErrorCode result = this->tag_get_tags_on_entity( handle, set_tags );
3688  std::cout << indent_prefix << ( tp == MB_TAG_SPARSE ? "Sparse tags:" : "Dense tags:" ) << std::endl;
3689  indent_prefix += " ";
3690 
3691  for( std::vector< Tag >::iterator vit = set_tags.begin(); vit != set_tags.end(); ++vit )
3692  {
3693  TagType this_type;
3694  result = this->tag_get_type( *vit, this_type );
3695  if( MB_SUCCESS != result || tp != this_type ) continue;
3696  DataType this_data_type;
3697  result = this->tag_get_data_type( *vit, this_data_type );
3698  if( MB_SUCCESS != result ) continue;
3699  int this_size;
3700  result = this->tag_get_length( *vit, this_size );
3701  if( MB_SUCCESS != result ) continue;
3702  // use double since this is largest single-valued tag
3703  std::vector< double > dbl_vals( this_size );
3704  std::vector< int > int_vals( this_size );
3705  std::vector< EntityHandle > hdl_vals( this_size );
3706  std::string tag_name;
3707  result = this->tag_get_name( *vit, tag_name );
3708  if( MB_SUCCESS != result ) continue;
3709  switch( this_data_type )
3710  {
3711  case MB_TYPE_INTEGER:
3712  result = this->tag_get_data( *vit, &handle, 1, &int_vals[0] );
3713  if( MB_SUCCESS != result ) continue;
3714  std::cout << indent_prefix << tag_name << " = ";
3715  if( this_size < 10 )
3716  for( int i = 0; i < this_size; i++ )
3717  std::cout << int_vals[i] << " ";
3718  else
3719  std::cout << int_vals[0] << "... (mult values)";
3720  std::cout << std::endl;
3721  break;
3722  case MB_TYPE_DOUBLE:
3723  result = this->tag_get_data( *vit, &handle, 1, &dbl_vals[0] );
3724  if( MB_SUCCESS != result ) continue;
3725  std::cout << indent_prefix << tag_name << " = ";
3726  if( this_size < 10 )
3727  for( int i = 0; i < this_size; i++ )
3728  std::cout << dbl_vals[i] << " ";
3729  else
3730  std::cout << dbl_vals[0] << "... (mult values)";
3731  std::cout << std::endl;
3732  break;
3733  case MB_TYPE_HANDLE:
3734  result = this->tag_get_data( *vit, &handle, 1, &hdl_vals[0] );
3735  if( MB_SUCCESS != result ) continue;
3736  std::cout << indent_prefix << tag_name << " = ";
3737  if( this_size < 10 )
3738  for( int i = 0; i < this_size; i++ )
3739  std::cout << hdl_vals[i] << " ";
3740  else
3741  std::cout << hdl_vals[0] << "... (mult values)";
3742  std::cout << std::endl;
3743  break;
3744  case MB_TYPE_OPAQUE:
3745  if( NAME_TAG_SIZE == this_size )
3746  {
3747  char dum_tag[NAME_TAG_SIZE];
3748  result = this->tag_get_data( *vit, &handle, 1, &dum_tag );
3749  if( MB_SUCCESS != result ) continue;
3750  // insert NULL just in case there isn't one
3751  dum_tag[NAME_TAG_SIZE - 1] = '\0';
3752  std::cout << indent_prefix << tag_name << " = " << dum_tag << std::endl;
3753  }
3754  break;
3755  case MB_TYPE_BIT:
3756  break;
3757  }
3758  }
3759 
3760  return MB_SUCCESS;
3761 }
3762 
3764 {
3765  // run through all entities, checking adjacencies and reverse-evaluating them
3766  Range all_ents;
3767  ErrorCode result = get_entities_by_handle( 0, all_ents );MB_CHK_ERR( result );
3768 
3769  for( Range::iterator rit = all_ents.begin(); rit != all_ents.end(); ++rit )
3770  {
3771  result = check_adjacencies( &( *rit ), 1 );MB_CHK_ERR( result );
3772  }
3773 
3774  return MB_SUCCESS;
3775 }
3776 
3777 ErrorCode Core::check_adjacencies( const EntityHandle* ents, int num_ents )
3778 {
3779 
3780  ErrorCode result = MB_SUCCESS, tmp_result;
3781  std::ostringstream oss;
3782 
3783  for( int i = 0; i < num_ents; i++ )
3784  {
3785  EntityHandle this_ent = ents[i];
3786  std::ostringstream ent_str;
3787  ent_str << CN::EntityTypeName( TYPE_FROM_HANDLE( this_ent ) ) << " " << ID_FROM_HANDLE( this_ent ) << ": ";
3788  int this_dim = dimension_from_handle( this_ent );
3789 
3790  if( !is_valid( this_ent ) )
3791  {
3792  std::cerr << ent_str.str() << "Not a valid entity." << std::endl;
3793  result = MB_FAILURE;
3794  }
3795 
3796  else
3797  {
3798  if( TYPE_FROM_HANDLE( this_ent ) == MBENTITYSET ) continue;
3799 
3800  // get adjacencies for this entity
3801  Range adjs;
3802  for( int dim = 0; dim <= 3; dim++ )
3803  {
3804  if( dim == this_dim ) continue;
3805  tmp_result = get_adjacencies( &this_ent, 1, dim, false, adjs, Interface::UNION );
3806  if( MB_SUCCESS != tmp_result )
3807  {
3808  oss << ent_str.str() << "Failed to get adjacencies for dimension " << dim << "." << std::endl;
3809  result = tmp_result;
3810  }
3811  }
3812  if( !oss.str().empty() )
3813  {
3814  std::cerr << oss.str();
3815  oss.str( "" );
3816  }
3817 
3818  // now check and reverse-evaluate them
3819  for( Range::iterator rit = adjs.begin(); rit != adjs.end(); ++rit )
3820  {
3821  EntitySequence* seq = 0;
3822  tmp_result = sequence_manager()->find( *rit, seq );
3823  if( seq == 0 || tmp_result != MB_SUCCESS )
3824  {
3825  oss << ent_str.str() << "Adjacent entity " << CN::EntityTypeName( TYPE_FROM_HANDLE( *rit ) ) << " "
3826  << ID_FROM_HANDLE( *rit ) << " is invalid." << std::endl;
3827  result = tmp_result;
3828  }
3829  else
3830  {
3831  Range rev_adjs;
3832  tmp_result = get_adjacencies( &( *rit ), 1, this_dim, false, rev_adjs );
3833  if( MB_SUCCESS != tmp_result )
3834  {
3835  oss << ent_str.str() << "Failed to get reverse adjacency from "
3836  << CN::EntityTypeName( TYPE_FROM_HANDLE( *rit ) ) << " " << ID_FROM_HANDLE( *rit );
3837  if( MB_MULTIPLE_ENTITIES_FOUND == tmp_result )
3838  oss << " (MULTIPLE)" << std::endl;
3839  else
3840  oss << " (" << tmp_result << ")" << std::endl;
3841  result = tmp_result;
3842  }
3843  else if( rev_adjs.find( this_ent ) == rev_adjs.end() )
3844  {
3845  oss << ent_str.str() << "Failed to find adjacency to this entity from "
3846  << CN::EntityTypeName( TYPE_FROM_HANDLE( *rit ) ) << " " << ID_FROM_HANDLE( *rit ) << "."
3847  << std::endl;
3848  result = tmp_result;
3849  }
3850  }
3851  if( !oss.str().empty() )
3852  {
3853  std::cerr << oss.str();
3854  oss.str( "" );
3855  }
3856  }
3857  }
3858  }
3859 
3860  return result;
3861 }
3862 
3863 bool Core::is_valid( const EntityHandle this_ent ) const
3864 {
3865  const EntitySequence* seq = 0;
3866  ErrorCode result = sequence_manager()->find( this_ent, seq );
3867  return seq != 0 && result == MB_SUCCESS;
3868 }
3869 
3871  EntityType ent_type,
3872  int ent_dim,
3873  int chunk_size,
3874  bool check_valid,
3875  SetIterator*& set_iter )
3876 {
3877  // check the type of set
3878  unsigned int setoptions;
3879  ErrorCode rval = MB_SUCCESS;
3880  if( meshset )
3881  {
3882  rval = get_meshset_options( meshset, setoptions );MB_CHK_ERR( rval );
3883  }
3884 
3885  if( !meshset || ( setoptions & MESHSET_SET ) )
3886  set_iter = new( std::nothrow ) RangeSetIterator( this, meshset, chunk_size, ent_type, ent_dim, check_valid );
3887  else
3888  set_iter = new( std::nothrow ) VectorSetIterator( this, meshset, chunk_size, ent_type, ent_dim, check_valid );
3889 
3890  setIterators.push_back( set_iter );
3891  return MB_SUCCESS;
3892 }
3893 
3894 /** \brief Remove the set iterator from the instance's list
3895  * This function is called from the SetIterator destructor, and should not be called directly
3896  * from anywhere else.
3897  * \param set_iter Set iterator being removed
3898  */
3900 {
3901  std::vector< SetIterator* >::iterator vit = std::find( setIterators.begin(), setIterators.end(), set_iter );
3902  if( vit == setIterators.end() )
3903  {
3904  MB_SET_ERR( MB_FAILURE, "Didn't find that iterator" );
3905  }
3906 
3907  setIterators.erase( vit );
3908 
3909  return MB_SUCCESS;
3910 }
3911 
3912 /** \brief Get all set iterators associated with the set passed in
3913  * \param meshset Meshset for which iterators are requested
3914  * \param set_iters Set iterators for the set
3915  */
3916 ErrorCode Core::get_set_iterators( EntityHandle meshset, std::vector< SetIterator* >& set_iters )
3917 {
3918  for( std::vector< SetIterator* >::const_iterator vit = setIterators.begin(); vit != setIterators.end(); ++vit )
3919  if( ( *vit )->ent_set() == meshset ) set_iters.push_back( *vit );
3920  return MB_SUCCESS;
3921 }
3922 
3924  type_memstorage* total_storage,
3925  type_memstorage* total_amortized_storage,
3926  type_memstorage* entity_storage,
3927  type_memstorage* amortized_entity_storage,
3928  type_memstorage* adjacency_storage,
3929  type_memstorage* amortized_adjacency_storage,
3930  const Tag* tag_array,
3931  unsigned num_tags,
3932  type_memstorage* tag_storage,
3933  type_memstorage* amortized_tag_storage )
3934 {
3935  // Figure out which values we need to calculate
3936  type_memstorage i_entity_storage, ia_entity_storage, i_adjacency_storage, ia_adjacency_storage, i_tag_storage,
3937  ia_tag_storage;
3938  type_memstorage *total_tag_storage = 0, *amortized_total_tag_storage = 0;
3939  if( !tag_array )
3940  {
3941  total_tag_storage = tag_storage;
3942  amortized_total_tag_storage = amortized_tag_storage;
3943  }
3944  if( total_storage || total_amortized_storage )
3945  {
3946  if( !entity_storage ) entity_storage = &i_entity_storage;
3947  if( !amortized_entity_storage ) amortized_entity_storage = &ia_entity_storage;
3948  if( !adjacency_storage ) adjacency_storage = &i_adjacency_storage;
3949  if( !amortized_adjacency_storage ) amortized_adjacency_storage = &ia_adjacency_storage;
3950  }
3951  else
3952  {
3953  if( entity_storage || amortized_entity_storage )
3954  {
3955  if( !amortized_entity_storage )
3956  amortized_entity_storage = &ia_entity_storage;
3957  else if( !entity_storage )
3958  entity_storage = &i_entity_storage;
3959  }
3960  if( adjacency_storage || amortized_adjacency_storage )
3961  {
3962  if( !amortized_adjacency_storage )
3963  amortized_adjacency_storage = &ia_adjacency_storage;
3964  else if( !adjacency_storage )
3965  adjacency_storage = &i_adjacency_storage;
3966  }
3967  }
3968  if( !total_tag_storage && total_storage ) total_tag_storage = &i_tag_storage;
3969  if( !amortized_total_tag_storage && total_amortized_storage ) amortized_total_tag_storage = &ia_tag_storage;
3970 
3971  // get entity storage
3972  if( amortized_entity_storage )
3973  {
3974  if( ents )
3975  sequenceManager->get_memory_use( *ents, *entity_storage, *amortized_entity_storage );
3976  else
3977  sequenceManager->get_memory_use( *entity_storage, *amortized_entity_storage );
3978  }
3979 
3980  // get adjacency storage
3981  if( amortized_adjacency_storage )
3982  {
3983  if( ents )
3984  aEntityFactory->get_memory_use( *ents, *adjacency_storage, *amortized_adjacency_storage );
3985  else
3986 #ifdef MOAB_HAVE_AHF
3987  ahfRep->get_memory_use( *adjacency_storage, *amortized_adjacency_storage );
3988 #else
3989  aEntityFactory->get_memory_use( *adjacency_storage, *amortized_adjacency_storage );
3990 #endif
3991  }
3992 
3993  // get storage for requested list of tags
3994  if( tag_array )
3995  {
3996  for( unsigned i = 0; i < num_tags; ++i )
3997  {
3998  if( !valid_tag_handle( tag_array[i] ) ) continue;
3999 
4000  unsigned long total = 0, per_ent = 0;
4001  tag_array[i]->get_memory_use( sequenceManager, total, per_ent );
4002 
4003  if( ents )
4004  {
4005  size_t count = 0, count2 = 0;
4006  tag_array[i]->num_tagged_entities( sequenceManager, count, MBMAXTYPE, ents );
4007  if( tag_storage ) tag_storage[i] = count * per_ent;
4008  if( amortized_tag_storage )
4009  {
4010  tag_array[i]->num_tagged_entities( sequenceManager, count2 );
4011  if( count2 )
4012  amortized_tag_storage[i] = static_cast< type_memstorage >( total * count * 1.0 / count2 );
4013  }
4014  }
4015  else
4016  {
4017  size_t count = 0;
4018  if( tag_storage )
4019  {
4020  tag_array[i]->num_tagged_entities( sequenceManager, count );
4021  tag_storage[i] = count * per_ent;
4022  }
4023  if( amortized_tag_storage ) amortized_tag_storage[i] = total;
4024  }
4025  }
4026  }
4027 
4028  // get storage for all tags
4029  if( total_tag_storage || amortized_total_tag_storage )
4030  {
4031  if( amortized_total_tag_storage ) *amortized_total_tag_storage = 0;
4032  if( total_tag_storage ) *total_tag_storage = 0;
4033 
4034  std::vector< Tag > tags;
4035  tag_get_tags( tags );
4036  for( std::list< TagInfo* >::const_iterator i = tagList.begin(); i != tagList.end(); ++i )
4037  {
4038  unsigned long total = 0, per_ent = 0;
4039  ( *i )->get_memory_use( sequenceManager, total, per_ent );
4040 
4041  if( ents )
4042  {
4043  size_t count = 0, count2 = 0;
4044  ( *i )->num_tagged_entities( sequenceManager, count, MBMAXTYPE, ents );
4045  if( total_tag_storage ) *total_tag_storage += count * per_ent;
4046  if( amortized_total_tag_storage )
4047  {
4048  ( *i )->num_tagged_entities( sequenceManager, count2 );
4049  if( count2 )
4050  *amortized_total_tag_storage += static_cast< type_memstorage >( total * count * 1.0 / count2 );
4051  }
4052  }
4053  else
4054  {
4055  size_t count = 0;
4056  if( total_tag_storage )
4057  {
4058  ( *i )->num_tagged_entities( sequenceManager, count );
4059  *total_tag_storage += count * per_ent;
4060  }
4061  if( amortized_total_tag_storage ) *amortized_total_tag_storage += total;
4062  }
4063  }
4064  }
4065 
4066  // calculate totals
4067  if( total_storage ) *total_storage = *entity_storage + *adjacency_storage + *total_tag_storage;
4068 
4069  if( total_amortized_storage )
4070  *total_amortized_storage =
4071  *amortized_entity_storage + *amortized_adjacency_storage + *amortized_total_tag_storage;
4072 }
4073 
4075  unsigned long num_ents,
4076  type_memstorage* total_storage,
4077  type_memstorage* total_amortized_storage,
4078  type_memstorage* entity_storage,
4079  type_memstorage* amortized_entity_storage,
4080  type_memstorage* adjacency_storage,
4081  type_memstorage* amortized_adjacency_storage,
4082  const Tag* tag_array,
4083  unsigned num_tags,
4084  type_memstorage* tag_storage,
4085  type_memstorage* amortized_tag_storage )
4086 {
4087  Range range;
4088 
4089  // If non-empty entity list, call range version of function
4090  if( ent_array )
4091  {
4092  if( num_ents > 20 )
4093  {
4094  std::vector< EntityHandle > list( num_ents );
4095  std::copy( ent_array, ent_array + num_ents, list.begin() );
4096  std::sort( list.begin(), list.end() );
4097  Range::iterator j = range.begin();
4098  for( std::vector< EntityHandle >::reverse_iterator i = list.rbegin(); i != list.rend(); ++i )
4099  j = range.insert( j, *i, *i );
4100  }
4101  else
4102  {
4103  std::copy( ent_array, ent_array + num_ents, range_inserter( range ) );
4104  }
4105  }
4106 
4107  estimated_memory_use_internal( ent_array ? &range : 0, total_storage, total_amortized_storage, entity_storage,
4108  amortized_entity_storage, adjacency_storage, amortized_adjacency_storage, tag_array,
4109  num_tags, tag_storage, amortized_tag_storage );
4110 }
4111 
4113  type_memstorage* total_storage,
4114  type_memstorage* total_amortized_storage,
4115  type_memstorage* entity_storage,
4116  type_memstorage* amortized_entity_storage,
4117  type_memstorage* adjacency_storage,
4118  type_memstorage* amortized_adjacency_storage,
4119  const Tag* tag_array,
4120  unsigned num_tags,
4121  type_memstorage* tag_storage,
4122  type_memstorage* amortized_tag_storage )
4123 {
4124  estimated_memory_use_internal( &ents, total_storage, total_amortized_storage, entity_storage,
4125  amortized_entity_storage, adjacency_storage, amortized_adjacency_storage, tag_array,
4126  num_tags, tag_storage, amortized_tag_storage );
4127 }
4128 
4130 {
4131  ErrorCode rval;
4133  const TypeSequenceManager& verts = sequence_manager()->entity_map( MBVERTEX );
4134  if( !verts.empty() )
4135  printf( " Vertex ID X Y Z Adjacencies \n"
4136  " ---------- -------- -------- -------- -----------...\n" );
4137  const EntityHandle* adj;
4138  int nadj;
4139  for( i = verts.begin(); i != verts.end(); ++i )
4140  {
4141  const VertexSequence* seq = static_cast< const VertexSequence* >( *i );
4142  printf( "(Sequence [%d,%d] in SequenceData [%d,%d])\n", (int)ID_FROM_HANDLE( seq->start_handle() ),
4143  (int)ID_FROM_HANDLE( seq->end_handle() ), (int)ID_FROM_HANDLE( seq->data()->start_handle() ),
4144  (int)ID_FROM_HANDLE( seq->data()->end_handle() ) );
4145 
4146  double c[3];
4147  for( EntityHandle h = seq->start_handle(); h <= seq->end_handle(); ++h )
4148  {
4149  rval = seq->get_coordinates( h, c );
4150  if( MB_SUCCESS == rval )
4151  printf( " %10d %8g %8g %8g", (int)ID_FROM_HANDLE( h ), c[0], c[1], c[2] );
4152  else
4153  printf( " %10d < ERROR %4d >", (int)ID_FROM_HANDLE( h ), (int)rval );
4154 
4155  rval = a_entity_factory()->get_adjacencies( h, adj, nadj );
4156  if( MB_SUCCESS != rval )
4157  {
4158  printf( " <ERROR %d>\n", (int)rval );
4159  continue;
4160  }
4161  EntityType pt = MBMAXTYPE;
4162  for( int j = 0; j < nadj; ++j )
4163  {
4164  if( TYPE_FROM_HANDLE( adj[j] ) != pt )
4165  {
4166  pt = TYPE_FROM_HANDLE( adj[j] );
4167  printf( " %s", pt >= MBMAXTYPE ? "INVALID TYPE" : CN::EntityTypeName( pt ) );
4168  }
4169  printf( " %d", (int)ID_FROM_HANDLE( adj[j] ) );
4170  }
4171  printf( "\n" );
4172  }
4173  }
4174 
4175  for( EntityType t = MBEDGE; t < MBENTITYSET; ++t )
4176  {
4177  const TypeSequenceManager& elems = sequence_manager()->entity_map( t );
4178  if( elems.empty() ) continue;
4179 
4180  int clen = 0;
4181  for( i = elems.begin(); i != elems.end(); ++i )
4182  {
4183  int n = static_cast< const ElementSequence* >( *i )->nodes_per_element();
4184  if( n > clen ) clen = n;
4185  }
4186 
4187  clen *= 5;
4188  if( clen < (int)strlen( "Connectivity" ) ) clen = strlen( "Connectivity" );
4189  std::vector< char > dashes( clen, '-' );
4190  dashes.push_back( '\0' );
4191  printf( " %7s ID %-*s Adjacencies\n", CN::EntityTypeName( t ), clen, "Connectivity" );
4192  printf( " ---------- %s -----------...\n", &dashes[0] );
4193 
4194  std::vector< EntityHandle > storage;
4195  const EntityHandle* conn;
4196  int nconn;
4197  for( i = elems.begin(); i != elems.end(); ++i )
4198  {
4199  const ElementSequence* seq = static_cast< const ElementSequence* >( *i );
4200  printf( "(Sequence [%d,%d] in SequenceData [%d,%d])\n", (int)ID_FROM_HANDLE( seq->start_handle() ),
4201  (int)ID_FROM_HANDLE( seq->end_handle() ), (int)ID_FROM_HANDLE( seq->data()->start_handle() ),
4202  (int)ID_FROM_HANDLE( seq->data()->end_handle() ) );
4203 
4204  for( EntityHandle h = seq->start_handle(); h <= seq->end_handle(); ++h )
4205  {
4206  printf( " %10d", (int)ID_FROM_HANDLE( h ) );
4207  rval = get_connectivity( h, conn, nconn, false, &storage );
4208  if( MB_SUCCESS != rval )
4209  printf( " <ERROR %2d>%*s", (int)rval, clen - 10, "" );
4210  else
4211  {
4212  for( int j = 0; j < nconn; ++j )
4213  printf( " %4d", (int)ID_FROM_HANDLE( conn[j] ) );
4214  printf( "%*s", clen - 5 * nconn, "" );
4215  }
4216 
4217  rval = a_entity_factory()->get_adjacencies( h, adj, nadj );
4218  if( MB_SUCCESS != rval )
4219  {
4220  printf( " <ERROR %d>\n", (int)rval );
4221  continue;
4222  }
4223  EntityType pt = MBMAXTYPE;
4224  for( int j = 0; j < nadj; ++j )
4225  {
4226  if( TYPE_FROM_HANDLE( adj[j] ) != pt )
4227  {
4228  pt = TYPE_FROM_HANDLE( adj[j] );
4229  printf( " %s", pt >= MBMAXTYPE ? "INVALID TYPE" : CN::EntityTypeName( pt ) );
4230  }
4231  printf( " %d", (int)ID_FROM_HANDLE( adj[j] ) );
4232  }
4233  printf( "\n" );
4234  }
4235  }
4236  }
4237 }
4238 
4240  const HomCoord& coord_max,
4241  EntityType entity_type,
4242  EntityID start_id_hint,
4243  EntityHandle& first_handle_out,
4244  EntitySequence*& sequence_out )
4245 {
4246  // NR: Previously, the structured element sequences were created via direct call to
4247  // the sequence manager instead of using the same from the ScdInterface which
4248  // creates the associated scd bounding box after element sequence creation.
4249 
4250  if( !scdInterface ) scdInterface = new ScdInterface( this );
4251  ScdBox* newBox = NULL;
4252  ErrorCode rval = scdInterface->create_scd_sequence( coord_min, coord_max, entity_type,
4253  /*starting_id*/ (int)start_id_hint, newBox );MB_CHK_ERR( rval );
4254 
4255  if( MBVERTEX == entity_type )
4256  first_handle_out = newBox->get_vertex( coord_min );
4257  else
4258  first_handle_out = newBox->get_element( coord_min );
4259  return sequence_manager()->find( first_handle_out, sequence_out );
4260 }
4261 
4263  EntitySequence* elem_seq,
4264  const HomCoord& p1,
4265  const HomCoord& q1,
4266  const HomCoord& p2,
4267  const HomCoord& q2,
4268  const HomCoord& p3,
4269  const HomCoord& q3,
4270  bool bb_input,
4271  const HomCoord* bb_min,
4272  const HomCoord* bb_max )
4273 {
4274  return sequence_manager()->add_vsequence( vert_seq, elem_seq, p1, q1, p2, q2, p3, q3, bb_input, bb_min, bb_max );
4275 }
4276 
4277 } // namespace moab