[MOAB-dev] r4641 - MOAB/trunk/itaps/imesh

[jvporter at wisc.edu]

Author: jvporter
Date: 2011-03-24 18:40:14 -0500 (Thu, 24 Mar 2011)
New Revision: 4641

Modified:
   MOAB/trunk/itaps/imesh/iMesh_MOAB.cpp
Log:
Fix list-list boolean intersection. The new algorithm is O(m+n) if TR1 is
available or O(n log n) if not.


Modified: MOAB/trunk/itaps/imesh/iMesh_MOAB.cpp
===================================================================
--- MOAB/trunk/itaps/imesh/iMesh_MOAB.cpp	2011-03-24 22:33:51 UTC (rev 4640)
+++ MOAB/trunk/itaps/imesh/iMesh_MOAB.cpp	2011-03-24 23:40:14 UTC (rev 4641)
@@ -14,6 +14,14 @@
 #include "moab_mpi.h"
 #endif
 
+#define STRINGIFY_(X) #X
+#define STRINGIFY(X) STRINGIFY_(X)
+#ifdef HAVE_UNORDERED_MAP
+# include STRINGIFY(HAVE_UNORDERED_MAP)
+#else
+# include <map>
+#endif
+
 #include <iostream>
 #include <cassert>
 #include <cctype>
@@ -2671,9 +2679,63 @@
       result = MOABI->unite_meshset(temp_set, set1);
     }
     else {
-      result = MOABI->unite_meshset(temp_set, set1);
-      if (MB_SUCCESS == result)
-        result = MOABI->intersect_meshset(temp_set, set2);
+      if (isList1 && isList2) {
+          // ITAPS has very specific rules about the behavior of intersection on
+          // list-type sets. Since MOAB doesn't (and likely will never) work
+          // exactly this way, we implement our own algorithm here.
+
+#ifdef HAVE_UNORDERED_MAP
+        typedef UNORDERED_MAP_NS::unordered_map<EntityHandle, size_t> lookup_t;
+#else
+        typedef std::map<EntityHandle, size_t> lookup_t;
+#endif
+
+          // First, build a lookup table for the second set.
+        lookup_t lookup;
+        {
+          std::vector<EntityHandle> contents2;
+          result = MOABI->get_entities_by_handle(set2, contents2);
+          CHKERR(result,"iMesh_intersect: ERROR intersect failed.");
+
+          for (std::vector<EntityHandle>::iterator i = contents2.begin();
+               i != contents2.end(); ++i) {
+#ifdef HAVE_UNORDERED_MAP
+            lookup_t::iterator j = lookup.find(*i);
+#else
+            lookup_t::iterator j = lookup.lower_bound(*i);
+#endif