1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
#include <iostream>
#include "meshkit/cleanup.hpp"
//#include "moab/AdaptiveKDTree.hpp"
#include "moab/OrientedBoxTreeTool.hpp"

using namespace moab;
namespace cleanup {
  // The obbtrees are no longer valid because the triangles have been altered.
  //  -Surface and volume sets are tagged with tags holding the obb tree
  //   root handles.
  //  -Surface/volume set handles are added to the root meshset.
  // Somehow, delete the old tree without deleting the
  // surface and volume sets, then build a new tree.
  ErrorCode remove_obb_tree(bool verbose) {<--- The function 'remove_obb_tree' is never used.
    ErrorCode result;
    Range obb_entities;
    Tag obbTag;
    result = MBI()->tag_get_handle( "OBB_TREE", sizeof(EntityHandle),
	       		  MB_TYPE_HANDLE, obbTag, MB_TAG_DENSE, NULL, 0 );
    if(verbose)
    {
    if(gen::error(MB_SUCCESS != result, "could not get OBB tree handle")) return result;
    }
    else
    {
    if(gen::error(MB_SUCCESS != result, "")) return result;
    }
    // This gets the surface/volume sets. I don't want to delete the sets.
    // I want to remove the obbTag that contains the tree root handle and
    // delete the tree.
    result = MBI()->get_entities_by_type_and_tag( 0, MBENTITYSET,
						    &obbTag, 0, 1, obb_entities );
    assert(MB_SUCCESS == result);
    std::cout << "  found " << obb_entities.size() << " OBB entities" << std::endl;
    //gen::print_range( obb_entities );
    //result = MBI()->delete_entities( obb_entities );
 
    // find tree roots
    Range trees;
    OrientedBoxTreeTool tool( MBI() );
    Tag rootTag;
    for(Range::iterator i=obb_entities.begin(); i!=obb_entities.end(); i++) {<--- Prefer prefix ++/-- operators for non-primitive types.
      EntityHandle root;
      result = MBI()->tag_get_data( obbTag, &(*i), 1, &root );
      if(gen::error(MB_SUCCESS!=result, "coule not get OBB tree data")) return result;
      //assert(MB_SUCCESS == result);
      tool.delete_tree( root );
    }
    result = MBI()->tag_delete( obbTag ); // use this for DENSE tags
    assert(MB_SUCCESS == result);


    result = MBI()->tag_get_handle ( "OBB", sizeof(double), 
                                  MB_TYPE_DOUBLE, rootTag, MB_TAG_SPARSE, 0, 0);
    assert(MB_SUCCESS==result || MB_ALREADY_ALLOCATED==result);
    /*    result = MBI()->get_entities_by_type_and_tag( 0, MBENTITYSET, &rootTag, 
                                                   NULL, 1, trees );
    if(MB_SUCCESS != result) std::cout << "result=" << result << std::endl;
    assert(MB_SUCCESS == result);
    //tool.find_all_trees( trees );
    std::cout << trees.size() << " tree(s) contained in file" << std::endl;
    //gen::print_range( trees );
  
    // delete the trees
    for (Range::iterator i = trees.begin(); i != trees.end(); ++i) {
      result = tool.delete_tree( *i );
      //if(MB_SUCCESS != result) std::cout << "result=" << result << std::endl;
      //assert(MB_SUCCESS == result);
    }
    */ 
    // Were all of the trees deleted? Perhaps some of the roots we found were
    // child roots that got deleted with their parents.
    trees.clear();
    result = MBI()->get_entities_by_type_and_tag( 0, MBENTITYSET, &rootTag,
						    NULL, 1, trees );
    assert(MB_SUCCESS == result);
    std::cout << "  " << trees.size() << " OBB tree(s) contained in file" << std::endl;
    return MB_SUCCESS;
  }

  ErrorCode delete_small_edge_and_tris( const EntityHandle vert0,<--- The function 'delete_small_edge_and_tris' is never used.
                                          EntityHandle &vert1,
                                          const double tol ) {
    // If the verts are the same, this is not meaningful.
    if(vert0 == vert1) return MB_SUCCESS;
    ErrorCode result = MB_SUCCESS;

    // If the edge is small, delete it and the adjacent tris.
    if(tol > gen::dist_between_verts(vert0, vert1)) {
      // get tris to delete
      Range tris;
      EntityHandle verts[2] = {vert0, vert1};
      result = MBI()->get_adjacencies( verts, 2, 2, false, tris );                   
      assert(MB_SUCCESS == result);
      result = MBI()->delete_entities( tris );
      assert(MB_SUCCESS == result);
      std::cout << "delete_small_edge_and_tris: deleted " << tris.size() 
                << " tris." << std::endl;
      // now merge the verts, keeping the first one
      // IN FUTURE, AVERAGE THE LOCATIONS???????????
      result = MBI()->merge_entities( vert0, vert1, false, true);
      assert(MB_SUCCESS == result);
      vert1 = vert0;
    }
    return result;
  }

  ErrorCode delete_small_edges(const Range &surfaces, const double FACET_TOL) {
    // PROBLEM: THIS IS INVALID BECAUSE TRIS CAN HAVE LONG EDGES BUT
    // SMALL AREA. All three pts are in a line. This is the nature of
    // faceting vs. meshing.
    /* Remove small triangles by removing edges that are too small. 
    Remove small edges by merging their endpoints together, creating
    degenerate triangles. Delete the degenerate triangles. */
    ErrorCode result;
    for(Range::const_iterator i=surfaces.begin(); i!=surfaces.end(); i++) {<--- Prefer prefix ++/-- operators for non-primitive types.
      std::cout << "surf_id=" << gen::geom_id_by_handle(*i) << std::endl;

      // get all tris
      Range tris;
      result = MBI()->get_entities_by_type( *i, MBTRI, tris );
      assert(MB_SUCCESS == result);

      // Check to ensure there area no degenerate tris
      for(Range::iterator j=tris.begin(); j!=tris.end(); j++) {<--- Prefer prefix ++/-- operators for non-primitive types.
        // get endpts
        const EntityHandle *endpts;
        int n_verts;                                      
        result = MBI()->get_connectivity( *j, endpts, n_verts);
        assert(MB_SUCCESS == result);
        assert(3 == n_verts);
        assert( endpts[0]!=endpts[1] && endpts[1]!=endpts[2] );
      }


      // get the skin first, because my find_skin does not check before creating edges.
      Range skin_edges;
      //result = gen::find_skin( tris, 1, skin_edges, false );
      Skinner tool(MBI());
      result = tool.find_skin( 0 , tris, 1, skin_edges, false );
      assert(MB_SUCCESS == result);

      // create the edges
      Range edges;
      result = MBI()->get_adjacencies( tris, 1, true, edges, Interface::UNION );
      if(MB_SUCCESS != result) {
	std::cout << "result=" << result << std::endl;
      }
      assert(MB_SUCCESS == result);

      // get the internal edges
      Range internal_edges = subtract(edges, skin_edges);

      for(Range::iterator j=internal_edges.begin(); j!=internal_edges.end(); j++) {<--- Prefer prefix ++/-- operators for non-primitive types.
        int n_internal_edges = internal_edges.size();
	std::cout << "edge=" << *j << std::endl;
        MBI()->list_entity( *j );
        assert(MB_SUCCESS == result);

        // get endpts
        const EntityHandle *endpts;
        int n_verts;                                      
        result = MBI()->get_connectivity( *j, endpts, n_verts);
        assert(MB_SUCCESS == result);
        assert(2 == n_verts);

        // does another edge exist w the same endpts? Why would it?
        Range duplicate_edges;
        result = MBI()->get_adjacencies( endpts, 2, 1, true, duplicate_edges );
        assert(MB_SUCCESS == result);
        if(1 < duplicate_edges.size()) MBI()->list_entities( duplicate_edges );
        assert(1 == duplicate_edges.size());

        // if the edge length is less than MERGE_TOL do nothing
        if(FACET_TOL < gen::dist_between_verts( endpts[0], endpts[1] )) continue; 
 
        // quick check
        for(Range::iterator k=internal_edges.begin(); k!=internal_edges.end(); k++) {<--- Prefer prefix ++/-- operators for non-primitive types.
          const EntityHandle *epts;
          int n_vts;                                      
          result = MBI()->get_connectivity( *k, epts, n_vts);
          assert(MB_SUCCESS == result);
          assert(2 == n_vts);
	  // The skin edges/verts cannot be moved, therefore both endpoints cannot 
	  // be on the skin. If they are, continue.
	  Range adj_edges0;
	  result = MBI()->get_adjacencies( &epts[0], 1, 1, true, adj_edges0 );
	  assert(MB_SUCCESS == result);
	  if(3 > adj_edges0.size()) {
	    std::cout << "adj_edges0.size()=" << adj_edges0.size() 
		      << " epts[0]=" << epts[0] << std::endl;
	    MBI()->list_entity( epts[0] );
	    //MBI()->write_mesh( "test_output.h5m" );
	    assert(MB_SUCCESS == result);
	  }
	  assert(3 <= adj_edges0.size());
	  Range adj_skin_edges0 = intersect( adj_edges0, skin_edges );
//	  bool endpt0_is_skin;
//	  if(adj_skin_edges0.empty()) endpt0_is_skin = false;
//	  else endpt0_is_skin = true;

	  Range adj_edges1;
	  result = MBI()->get_adjacencies( &epts[1], 1, 1, true, adj_edges1 );
	  assert(MB_SUCCESS == result);
	  if(3 > adj_edges1.size()) {
	    std::cout << "adj_edges1.size()=" << adj_edges1.size() 
		      << " epts[1]=" << epts[1] << std::endl;
	    MBI()->list_entity( epts[1] );
	    //MBI()->write_mesh( "test_output.h5m" );
	    assert(MB_SUCCESS == result);
	  }
	  assert(3 <= adj_edges1.size());
        }


        // The skin edges/verts cannot be moved, therefore both endpoints cannot 
        // be on the skin. If they are, continue.
        Range adj_edges0;
        result = MBI()->get_adjacencies( &endpts[0], 1, 1, true, adj_edges0 );
        assert(MB_SUCCESS == result);
        if(3 > adj_edges0.size()) {
          std::cout << "adj_edges0.size()=" << adj_edges0.size() 
                    << " endpts[0]=" << endpts[0] << std::endl;
          MBI()->list_entity( endpts[0] );
          //MBI()->write_mesh( "test_output.h5m" );
          assert(MB_SUCCESS == result);
        }
        assert(3 <= adj_edges0.size());
        Range adj_skin_edges0 = intersect( adj_edges0, skin_edges );
        bool endpt0_is_skin;
        if(adj_skin_edges0.empty()) endpt0_is_skin = false;
        else endpt0_is_skin = true;

        Range adj_edges1;
        result = MBI()->get_adjacencies( &endpts[1], 1, 1, true, adj_edges1 );
        assert(MB_SUCCESS == result);
        if(3 > adj_edges1.size()) {
          std::cout << "adj_edges1.size()=" << adj_edges1.size() 
                    << " endpts[1]=" << endpts[1] << std::endl;
          MBI()->list_entity( endpts[1] );
          //MBI()->write_mesh( "test_output.h5m" );
          assert(MB_SUCCESS == result);
        }
        assert(3 <= adj_edges1.size());
        Range adj_skin_edges1 = intersect( adj_edges1, skin_edges );
        bool endpt1_is_skin;
        if(adj_skin_edges1.empty()) endpt1_is_skin = false;
        else endpt1_is_skin = true;
        if(endpt0_is_skin && endpt1_is_skin) continue;
        
        // Keep the skin endpt, and delete the other endpt
        EntityHandle keep_endpt, delete_endpt;
        if(endpt0_is_skin) {
          keep_endpt   = endpts[0];
          delete_endpt = endpts[1];
        } else {
          keep_endpt   = endpts[1];
          delete_endpt = endpts[0];
        }

        // get the adjacent tris
	std::vector<EntityHandle> adj_tris;
        result = MBI()->get_adjacencies( &(*j), 1, 2, false, adj_tris );
        assert(MB_SUCCESS == result);
	        if(2 != adj_tris.size()) {
	std::cout << "adj_tris.size()=" << adj_tris.size() << std::endl;
	for(unsigned int i=0; i<adj_tris.size(); i++) gen::print_triangle( adj_tris[i], true );
        } 
        assert(2 == adj_tris.size());
 
        // When merging away an edge, a tri and 2 edges will be deleted.
        // Get each triangle's edge other edge the will be deleted.
        Range tri0_delete_edge_verts;
        result = MBI()->get_adjacencies( &adj_tris[0], 1, 0, true, tri0_delete_edge_verts );
        assert(MB_SUCCESS == result);
        assert(3 == tri0_delete_edge_verts.size());
        tri0_delete_edge_verts.erase( keep_endpt );
        Range tri0_delete_edge;
        result = MBI()->get_adjacencies( tri0_delete_edge_verts, 1, true, tri0_delete_edge );
        assert(MB_SUCCESS == result);
        assert(1 == tri0_delete_edge.size());      
 
        Range tri1_delete_edge_verts;
        result = MBI()->get_adjacencies( &adj_tris[1], 1, 0, true, tri1_delete_edge_verts );
        assert(MB_SUCCESS == result);
        assert(3 == tri1_delete_edge_verts.size());
        tri1_delete_edge_verts.erase( keep_endpt );
        Range tri1_delete_edge;
        result = MBI()->get_adjacencies( tri1_delete_edge_verts, 1, true, tri1_delete_edge );
        assert(MB_SUCCESS == result);
        assert(1 == tri1_delete_edge.size());      

        // When an edge is merged, it will be deleted and its to adjacent tris
        // will be deleted because they are degenerate. We cannot alter the skin.
        // How many skin edges does tri0 have?
	/*        Range tri_edges;
        result = MBI()->get_adjacencies( &adj_tris[0], 1, 1, false, tri_edges );
        assert(MB_SUCCESS == result);
        assert(3 == tri_edges.size());
        Range tri0_internal_edges = intersect(tri_edges, internal_edges);

        // Cannot merge the edge away if the tri has more than one skin edge.
        // Otherwise we would delete a skin edge. We already know the edges in 
        // edge_set are not on the skin.
        if(2 > tri0_internal_edges.size()) continue;

        // check the other tri
        tri_edges.clear();
        result = MBI()->get_adjacencies( &adj_tris[1], 1, 1, false, tri_edges );
        assert(MB_SUCCESS == result);
        assert(3 == tri_edges.size());
        Range tri1_internal_edges = intersect(tri_edges, internal_edges);
        if(2 > tri1_internal_edges.size()) continue;

        // Check to make sure that the internal edges are not on the skin
        Range temp;
        temp = intersect( tri0_internal_edges, skin_edges );
        assert(temp.empty());
        temp = intersect( tri1_internal_edges, skin_edges );
        assert(temp.empty());

        // We know that the edge will be merged away. Find the keep_vert and
        // delete_vert. The delete_vert should never be a skin vertex because the
        // skin must not move.
        Range delete_vert;
        result = MBI()->get_adjacencies( tri0_internal_edges, 0, false, delete_vert);
        assert(MB_SUCCESS == result);
        assert(1 == delete_vert.size());        
	*/

        // *********************************************************************
        // Test to see if the merge would create inverted tris. Several special
        // cases to avoid all result in inverted tris.
        // *********************************************************************
        // get all the tris adjacent to the point the will be moved.
        Range altered_tris;
        result = MBI()->get_adjacencies( &delete_endpt, 1, 2, false, altered_tris );
        assert(MB_SUCCESS == result);
        bool inverted_tri = false;
        for(Range::const_iterator k=altered_tris.begin(); k!=altered_tris.end(); ++k) {
          const EntityHandle *conn;
          int n_verts;
          result = MBI()->get_connectivity( *k, conn, n_verts );
          assert(MB_SUCCESS == result);
          assert(3 == tris.size());
          EntityHandle new_conn[3];
          for(unsigned int i=0; i<3; ++i) {
            new_conn[i] = (conn[i]==delete_endpt) ? keep_endpt : conn[i];
          }
          double area;
          result = gen::triangle_area( new_conn, area );
          assert(MB_SUCCESS == result);
          if(0 > area) {
	    std::cout << "inverted tri detected, area=" << area << std::endl;
            inverted_tri = true;
            break;
          }
        }
        if(inverted_tri) continue;      

        // If we got here, then the merge will occur. Delete the edge the will be
        // made degenerate and the edge that will be made duplicate.
	std::cout << "A merge will occur" << std::endl;
      	gen::print_triangle( adj_tris[0], true );
	gen::print_triangle( adj_tris[1], true );
        //tri0_internal_edges.erase( *j );
        //tri1_internal_edges.erase( *j );
        internal_edges.erase( tri0_delete_edge.front() );
        internal_edges.erase( tri1_delete_edge.front() );
	std::cout << "merged verts=" << keep_endpt << " " << delete_endpt << std::endl;
	MBI()->list_entity( keep_endpt );
	MBI()->list_entity( delete_endpt );
        result = MBI()->merge_entities( keep_endpt, delete_endpt, false, true );          
	assert(MB_SUCCESS == result);
        result = MBI()->delete_entities( tri0_delete_edge );
        assert(MB_SUCCESS == result);
        result = MBI()->delete_entities( tri1_delete_edge );
        assert(MB_SUCCESS == result);
        result = MBI()->delete_entities( &(*j), 1 );
        assert(MB_SUCCESS == result);
	std::cout << "deleted edges=" << *j << " " << tri0_delete_edge.front()
                  << " " << tri1_delete_edge.front() << std::endl;
	          
	// delete degenerate tris                          
	result = MBI()->delete_entities( &adj_tris[0], 2 );                 
	assert(MB_SUCCESS == result);
	MBI()->list_entity( keep_endpt );

        // remove the edge from the range
        j = internal_edges.erase(*j) - 1; 
	std::cout << "next iter=" << *j << std::endl;        

        Range new_tris;
        result = MBI()->get_entities_by_type( *i, MBTRI, new_tris );
        assert(MB_SUCCESS == result);
        Range new_skin_edges;
        result = tool.find_skin( *i, new_tris, 1, new_skin_edges, false );
        assert(MB_SUCCESS == result);
        assert(skin_edges.size() == new_skin_edges.size());
        for(unsigned int k=0; k<skin_edges.size(); k++) {
          if(skin_edges[k] != new_skin_edges[k]) {
            MBI()->list_entity( skin_edges[k] );
            MBI()->list_entity( new_skin_edges[k] );
          }
          assert(skin_edges[k] == new_skin_edges[k]);
        }
        if (n_internal_edges != (int) internal_edges.size()+3)
        assert(n_internal_edges == (int) internal_edges.size()+3);
      }
  
      // cleanup edges
      result = MBI()->get_entities_by_type( 0, MBEDGE, edges );
      assert(MB_SUCCESS == result);
      result = MBI()->delete_entities( edges ); 
      assert(MB_SUCCESS == result);
   
    }
    return MB_SUCCESS;
  } 
  
  // Lots of edges have been created but are no longer needed.
  // Delete edges that are not in curves. These should be the only edges
  // that remain. This incredibly speeds up the watertight_check tool (100x?).
  ErrorCode cleanup_edges( Range curve_meshsets ) {<--- The function 'cleanup_edges' is never used.
    ErrorCode result;
    Range edges, edges_to_keep;
    for(Range::iterator i=curve_meshsets.begin(); i!=curve_meshsets.end(); i++) {<--- Prefer prefix ++/-- operators for non-primitive types.
      result = MBI()->get_entities_by_dimension( *i, 1, edges );
      assert(MB_SUCCESS == result);
      edges_to_keep.merge( edges );
    }

    Range all_edges;
    result = MBI()->get_entities_by_dimension( 0, 1, all_edges );
    assert(MB_SUCCESS == result);

    // delete the edges that are not in curves.
    //Range edges_to_delete = all_edges.subtract( edges_to_keep );
    Range edges_to_delete = subtract( all_edges, edges_to_keep );
    std::cout << "deleting " << edges_to_delete.size() << " unused edges" << std::endl;
    result = MBI()->delete_entities( edges_to_delete );
    assert(MB_SUCCESS == result);
    return result;
  }
 
  
}