[MOAB-dev] r1478 - MOAB/trunk

tautges at mcs.anl.gov tautges at mcs.anl.gov
Wed Dec 19 11:59:01 CST 2007


Author: tautges
Date: 2007-12-19 11:59:01 -0600 (Wed, 19 Dec 2007)
New Revision: 1478

Modified:
   MOAB/trunk/GeomTopoTool.cpp
   MOAB/trunk/Tqdcfr.cpp
   MOAB/trunk/Tqdcfr.hpp
   MOAB/trunk/WriteSLAC.cpp
Log:
Fixed bug in cub reader for meshes with non-contiguous element numbering.

Changed algorithm used to restore geometric topology from elements, at the cost of memory.  For a 
large geometry with relatively few elements the cost is maybe 200MB; could be higher for a large mesh,
so I'll eventually change this to be controllable.

Fixed a few compile warnings in WriteSLAC.

Reading in a large cub file (> 5k volumes, > 20k surfaces, ~100k hexes) takes < 20 seconds, unoptimized.


Modified: MOAB/trunk/GeomTopoTool.cpp
===================================================================
--- MOAB/trunk/GeomTopoTool.cpp	2007-12-19 16:15:03 UTC (rev 1477)
+++ MOAB/trunk/GeomTopoTool.cpp	2007-12-19 17:59:01 UTC (rev 1478)
@@ -50,11 +50,31 @@
   result = separate_by_dimension(geom_sets, entities, geom_tag);
   if (MB_SUCCESS != result)
     return result;
+
+  std::vector<MBEntityHandle> parents;
+  MBRange tmp_parents;
   
-  std::vector<MBEntityHandle> dp1ents;
-
     // loop over dimensions
   for (int dim = 2; dim >= 0; dim--) {
+      // mark entities of next higher dimension with their owners; regenerate tag
+      // each dimension so prev dim's tag data goes away
+    MBTag owner_tag;
+    MBEntityHandle dum_val = 0;
+    result = mdbImpl->tag_create("__owner_tag", sizeof(MBEntityHandle), MB_TAG_DENSE,
+                                 MB_TYPE_HANDLE, owner_tag, &dum_val);
+    if (MB_SUCCESS != result) continue;
+    MBRange dp1ents;
+    std::vector<MBEntityHandle> owners;
+    for (MBRange::iterator rit = entities[dim+1].begin(); rit != entities[dim+1].end(); rit++) {
+      dp1ents.clear();
+      result = mdbImpl->get_entities_by_dimension(*rit, dim+1, dp1ents);
+      if (MB_SUCCESS != result) continue;
+      owners.resize(dp1ents.size());
+      std::fill(owners.begin(), owners.end(), *rit);
+      result = mdbImpl->tag_set_data(owner_tag, dp1ents, &owners[0]);
+      if (MB_SUCCESS != result) continue;
+    }
+    
     for (MBRange::iterator d_it = entities[dim].begin(); 
          d_it != entities[dim].end(); d_it++) {
       MBRange dents;
@@ -68,23 +88,25 @@
                                         false, dp1ents);
       if (MB_SUCCESS != result || dp1ents.empty()) continue;
 
-    //   . for each geom entity if dim d+1, if it contains any of the ents,
-    //     add it to list of parents
-
-      MBRange parents;      
-      for (MBRange::iterator git = entities[dim+1].begin(); 
-           git != entities[dim+1].end(); git++) {
-        if (mdbImpl->contains_entities(*git, &dp1ents[0], 
-                                       dp1ents.size()))
-          parents.insert(*git);
-      }
+        // get owner tags
+      parents.resize(dp1ents.size());
+      result = mdbImpl->tag_get_data(owner_tag, dp1ents, &parents[0]);
+      assert(MB_TAG_NOT_FOUND != result);
+      if (MB_SUCCESS != result) continue;
       
-    //   . make parent/child links with parents
-      for (MBRange::iterator pit = parents.begin(); pit != parents.end(); pit++) {
+        // compress to a range to remove duplicates
+      tmp_parents.clear();
+      std::copy(parents.begin(), parents.end(), mb_range_inserter(tmp_parents));
+      for (MBRange::iterator pit = tmp_parents.begin(); pit != tmp_parents.end(); pit++) {
         result = mdbImpl->add_parent_child(*pit, *d_it);
+        if (MB_SUCCESS != result) return result;
       }
     }
     
+      // now delete owner tag on this dimension, automatically removes tag data
+    result = mdbImpl->tag_delete(owner_tag);
+    if (MB_SUCCESS != result) return result;
+    
   } // dim
 
   return result;

Modified: MOAB/trunk/Tqdcfr.cpp
===================================================================
--- MOAB/trunk/Tqdcfr.cpp	2007-12-19 16:15:03 UTC (rev 1477)
+++ MOAB/trunk/Tqdcfr.cpp	2007-12-19 17:59:01 UTC (rev 1478)
@@ -862,8 +862,10 @@
         ent_list->push_back(int_buf[i]+currNodeIdOffset);
     }
     else {
-      for (int i = 0; i < id_buf_size; i++)
+      for (int i = 0; i < id_buf_size; i++) {
+        assert(0 != (*cubMOABVertexMap)[int_buf[i]]);
         ent_list->push_back((*cubMOABVertexMap)[int_buf[i]]);
+      }
     }
   }
   else {
@@ -897,7 +899,10 @@
                                Tqdcfr::ModelEntry *model,
                                Tqdcfr::GeomHeader *entity) 
 {
-  if (entity->nodeCt == 0) return MB_SUCCESS;
+  if (entity->nodeCt == 0) {
+    if (debug) std::cout << "(no nodes) ";
+    return MB_SUCCESS;
+  }
   
     // get the ids & coords in separate calls to minimize memory usage
     // position the file
@@ -931,8 +936,9 @@
   result = mdbImpl->tag_set_data(globalIdTag, dum_range, &int_buf[0]);
   if (MB_SUCCESS != result) return result;
 
-    // check for id contiguity
-  long unsigned int node_offset = mdbImpl->id_from_handle( node_handle);
+    // check for id contiguity; long because we subtract 1st index later, which
+    // might make this offset negative
+  long node_offset = mdbImpl->id_from_handle( node_handle);
 
   int max_id = -1;
   int contig;
@@ -948,26 +954,29 @@
   
     if (contig && -1 == currNodeIdOffset)
       currNodeIdOffset = node_offset;
-    else if ((contig && (long unsigned int) currNodeIdOffset != node_offset) ||
+    else if ((contig && currNodeIdOffset != node_offset) ||
              !contig) {
       // node offsets no longer valid - need to build cub id - vertex handle map
-      cubMOABVertexMap = new std::vector<MBEntityHandle>(node_offset+entity->nodeCt);
+      MBRange vrange;
+      result = mdbImpl->get_entities_by_type(0, MBVERTEX, vrange); RR;
+#define MAX(a,b) (a > b ? a : b)
+      int map_size = MAX((node_offset+entity->nodeCt), ((int)*vrange.rbegin()+1));
+      cubMOABVertexMap = new std::vector<MBEntityHandle>(map_size);
+      std::fill(cubMOABVertexMap->begin(), cubMOABVertexMap->end(), 0);
 
-      if (-1 != currNodeIdOffset && currNodeIdOffset != (int) node_offset) {
-          // now fill the missing values for vertices which already exist
-        MBRange vrange;
-        result = mdbImpl->get_entities_by_type(0, MBVERTEX, vrange); RR;
-        MBRange::const_iterator rit = vrange.lower_bound(vrange.begin(), vrange.end(),
-                                                         currNodeIdOffset);
-        for (; *rit <= node_offset; rit++)
+        // now fill the missing values for vertices which already exist
+      for (MBRange::iterator rit = vrange.begin(); rit != vrange.end(); rit++)
           (*cubMOABVertexMap)[*rit] = *rit;
-      }
     }
   }
   
   if (NULL != cubMOABVertexMap) {
       // expand the size if necessary
-    if (max_id > (int) cubMOABVertexMap->size()-1) cubMOABVertexMap->resize(max_id+1);
+    if (max_id > (int) cubMOABVertexMap->size()-1) {
+      long old_size = cubMOABVertexMap->size();
+      cubMOABVertexMap->resize(max_id+1);
+      std::fill(&(*cubMOABVertexMap)[old_size], &(*cubMOABVertexMap)[0]+cubMOABVertexMap->size(), 0);
+    }
     
       // now set the new values
     std::vector<int>::iterator vit;
@@ -978,6 +987,10 @@
     }
   }
 
+  if (debug)
+    std::cout << "(currOffset=" << currNodeIdOffset << ",nodeh=" << node_handle
+              << ",nodect=" << entity->nodeCt << ")  ";
+
     // set the dimension to at least zero (entity has at least nodes) on the geom tag
   int max_dim = 0;
   result = mdbImpl->tag_set_data(geomTag, &(entity->setHandle), 1, &max_dim);
@@ -1094,7 +1107,14 @@
     for (i = 0; i < total_conn; i++) {
       if (NULL == cubMOABVertexMap)
         new_handle = CREATE_HANDLE(MBVERTEX, currNodeIdOffset+tmp_conn[i], dum_err);
-      else new_handle = (*cubMOABVertexMap)[tmp_conn[i]];
+      else {
+        if (debug) {
+          if (0 == i) std::cout << "Conn";
+          std::cout << ", " << tmp_conn[i];
+        }
+        assert(0 != (*cubMOABVertexMap)[tmp_conn[i]]);
+        new_handle = (*cubMOABVertexMap)[tmp_conn[i]];
+      }
       assert(MB_SUCCESS == 
              mdbImpl->handle_from_id(MBVERTEX, mdbImpl->id_from_handle(new_handle), 
                                      dum_handle));

Modified: MOAB/trunk/Tqdcfr.hpp
===================================================================
--- MOAB/trunk/Tqdcfr.hpp	2007-12-19 16:15:03 UTC (rev 1477)
+++ MOAB/trunk/Tqdcfr.hpp	2007-12-19 17:59:01 UTC (rev 1478)
@@ -273,8 +273,8 @@
   FileTOC fileTOC;
   ModelEntry *modelEntries;
   MetaDataContainer modelMetaData;
-  int currNodeIdOffset;
-  int currElementIdOffset[MBMAXTYPE];
+  long currNodeIdOffset;
+  long currElementIdOffset[MBMAXTYPE];
   MBTag globalIdTag, cubIdTag, geomTag, uniqueIdTag, blockTag, nsTag, ssTag,
     attribVectorTag, entityNameTag, categoryTag;
   std::map<int, MBEntityHandle> uidSetMap;

Modified: MOAB/trunk/WriteSLAC.cpp
===================================================================
--- MOAB/trunk/WriteSLAC.cpp	2007-12-19 16:15:03 UTC (rev 1477)
+++ MOAB/trunk/WriteSLAC.cpp	2007-12-19 17:59:01 UTC (rev 1478)
@@ -691,7 +691,7 @@
     // first write the interior hexes
   NcVar *hex_conn = NULL;
   if (mesh_info.bdy_hexes.size() != 0 || mesh_info.num_int_hexes != 0) {
-    char *hex_name = "hexahedron_interior";
+    const char *hex_name = "hexahedron_interior";
     hex_conn = ncFile->get_var(hex_name);
     if (NULL == hex_conn) return MB_FAILURE;
   }
@@ -731,7 +731,7 @@
 
   NcVar *tet_conn = NULL;
   if (mesh_info.bdy_tets.size() != 0 || mesh_info.num_int_tets != 0) {
-    char *tet_name = "tetrahedron_interior";
+    const char *tet_name = "tetrahedron_interior";
     tet_conn = ncFile->get_var(tet_name);
     if (NULL == tet_conn) return MB_FAILURE;
   }
@@ -769,7 +769,7 @@
   
     // now the exterior hexes
   if (mesh_info.bdy_hexes.size() != 0) {
-    char *hex_name = "hexahedron_exterior";
+    const char *hex_name = "hexahedron_exterior";
     hex_conn = ncFile->get_var(hex_name);
     if (NULL == hex_conn) return MB_FAILURE;
 
@@ -820,7 +820,7 @@
 
     // now the exterior tets
   if (mesh_info.bdy_tets.size() != 0) {
-    char *tet_name = "tetrahedron_exterior";
+    const char *tet_name = "tetrahedron_exterior";
     tet_conn = ncFile->get_var(tet_name);
     if (NULL == tet_conn) return MB_FAILURE;
 




More information about the moab-dev mailing list