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