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
 448
 449
 450
 451
 452
 453
 454
 455
 456
 457
 458
 459
 460
 461
 462
 463
 464
 465
 466
 467
 468
 469
 470
 471
 472
 473
 474
 475
 476
 477
 478
 479
 480
 481
 482
 483
 484
 485
 486
 487
 488
 489
 490
 491
 492
 493
 494
 495
 496
 497
 498
 499
 500
 501
 502
 503
 504
 505
 506
 507
 508
 509
 510
 511
 512
 513
 514
 515
 516
 517
 518
 519
 520
 521
 522
 523
 524
 525
 526
 527
 528
 529
 530
 531
 532
 533
 534
 535
 536
 537
 538
 539
 540
 541
 542
 543
 544
 545
 546
 547
 548
 549
 550
 551
 552
 553
 554
 555
 556
 557
 558
 559
 560
 561
 562
 563
 564
 565
 566
 567
 568
 569
 570
 571
 572
 573
 574
 575
 576
 577
 578
 579
 580
 581
 582
 583
 584
 585
 586
 587
 588
 589
 590
 591
 592
 593
 594
 595
 596
 597
 598
 599
 600
 601
 602
 603
 604
 605
 606
 607
 608
 609
 610
 611
 612
 613
 614
 615
 616
 617
 618
 619
 620
 621
 622
 623
 624
 625
 626
 627
 628
 629
 630
 631
 632
 633
 634
 635
 636
 637
 638
 639
 640
 641
 642
 643
 644
 645
 646
 647
 648
 649
 650
 651
 652
 653
 654
 655
 656
 657
 658
 659
 660
 661
 662
 663
 664
 665
 666
 667
 668
 669
 670
 671
 672
 673
 674
 675
 676
 677
 678
 679
 680
 681
 682
 683
 684
 685
 686
 687
 688
 689
 690
 691
 692
 693
 694
 695
 696
 697
 698
 699
 700
 701
 702
 703
 704
 705
 706
 707
 708
 709
 710
 711
 712
 713
 714
 715
 716
 717
 718
 719
 720
 721
 722
 723
 724
 725
 726
 727
 728
 729
 730
 731
 732
 733
 734
 735
 736
 737
 738
 739
 740
 741
 742
 743
 744
 745
 746
 747
 748
 749
 750
 751
 752
 753
 754
 755
 756
 757
 758
 759
 760
 761
 762
 763
 764
 765
 766
 767
 768
 769
 770
 771
 772
 773
 774
 775
 776
 777
 778
 779
 780
 781
 782
 783
 784
 785
 786
 787
 788
 789
 790
 791
 792
 793
 794
 795
 796
 797
 798
 799
 800
 801
 802
 803
 804
 805
 806
 807
 808
 809
 810
 811
 812
 813
 814
 815
 816
 817
 818
 819
 820
 821
 822
 823
 824
 825
 826
 827
 828
 829
 830
 831
 832
 833
 834
 835
 836
 837
 838
 839
 840
 841
 842
 843
 844
 845
 846
 847
 848
 849
 850
 851
 852
 853
 854
 855
 856
 857
 858
 859
 860
 861
 862
 863
 864
 865
 866
 867
 868
 869
 870
 871
 872
 873
 874
 875
 876
 877
 878
 879
 880
 881
 882
 883
 884
 885
 886
 887
 888
 889
 890
 891
 892
 893
 894
 895
 896
 897
 898
 899
 900
 901
 902
 903
 904
 905
 906
 907
 908
 909
 910
 911
 912
 913
 914
 915
 916
 917
 918
 919
 920
 921
 922
 923
 924
 925
 926
 927
 928
 929
 930
 931
 932
 933
 934
 935
 936
 937
 938
 939
 940
 941
 942
 943
 944
 945
 946
 947
 948
 949
 950
 951
 952
 953
 954
 955
 956
 957
 958
 959
 960
 961
 962
 963
 964
 965
 966
 967
 968
 969
 970
 971
 972
 973
 974
 975
 976
 977
 978
 979
 980
 981
 982
 983
 984
 985
 986
 987
 988
 989
 990
 991
 992
 993
 994
 995
 996
 997
 998
 999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374
1375
1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401
1402
1403
1404
1405
1406
1407
1408
1409
1410
1411
1412
1413
1414
1415
1416
1417
1418
1419
1420
1421
1422
1423
1424
1425
1426
1427
1428
1429
1430
1431
1432
1433
1434
1435
1436
1437
1438
1439
1440
1441
1442
1443
1444
1445
1446
1447
1448
1449
1450
1451
1452
1453
1454
1455
1456
1457
1458
1459
1460
1461
1462
1463
1464
1465
1466
1467
1468
1469
1470
1471
1472
1473
1474
1475
1476
1477
1478
1479
1480
1481
1482
1483
1484
1485
1486
1487
1488
1489
1490
1491
1492
1493
1494
1495
1496
1497
1498
1499
1500
1501
1502
1503
1504
1505
1506
1507
1508
1509
1510
1511
1512
1513
/*
 * Intx2Mesh.cpp
 *
 *  Created on: Oct 2, 2012
 */

#include "moab/IntxMesh/Intx2Mesh.hpp"
#ifdef MOAB_HAVE_MPI
#include "moab/ParallelComm.hpp"
#include "MBParallelConventions.h"
#include "moab/ParallelMergeMesh.hpp"
#endif /* MOAB_HAVE_MPI */
#include "MBTagConventions.hpp"
// this is for DBL_MAX
#include <cfloat>
#include <queue>
#include <sstream>
#include "moab/GeomUtil.hpp"
#include "moab/AdaptiveKDTree.hpp"

namespace moab
{

#ifdef ENABLE_DEBUG
int Intx2Mesh::dbg_1 = 0;
#endif

Intx2Mesh::Intx2Mesh( Interface* mbimpl )
    : mb( mbimpl ), mbs1( 0 ), mbs2( 0 ), outSet( 0 ), gid( 0 ), TgtFlagTag( 0 ), tgtParentTag( 0 ), srcParentTag( 0 ),
      countTag( 0 ), srcNeighTag( 0 ), tgtNeighTag( 0 ), neighTgtEdgeTag( 0 ), orgSendProcTag( 0 ), tgtConn( NULL ),
      srcConn( NULL ), epsilon_1( 0.0 ), epsilon_area( 0.0 ), box_error( 0.0 ), localRoot( 0 ), my_rank( 0 )
#ifdef MOAB_HAVE_MPI
      ,
      parcomm( NULL ), remote_cells( NULL ), remote_cells_with_tracers( NULL )
#endif
      ,
      max_edges_1( 0 ), max_edges_2( 0 ), counting( 0 )
{
    gid = mbimpl->globalId_tag();
}

Intx2Mesh::~Intx2Mesh()
{
    // TODO Auto-generated destructor stub
#ifdef MOAB_HAVE_MPI
    if( remote_cells )
    {
        delete remote_cells;
        remote_cells = NULL;
    }
#endif
}
ErrorCode Intx2Mesh::FindMaxEdgesInSet( EntityHandle eset, int& max_edges )
{
    Range cells;
    ErrorCode rval = mb->get_entities_by_dimension( eset, 2, cells );MB_CHK_ERR( rval );

    max_edges = 0;  // can be 0 for point clouds
    for( Range::iterator cit = cells.begin(); cit != cells.end(); cit++ )
    {
        EntityHandle cell = *cit;
        const EntityHandle* conn4;
        int nnodes = 3;
        rval       = mb->get_connectivity( cell, conn4, nnodes );MB_CHK_SET_ERR( rval, "can't get connectivity of a cell" );
        if( nnodes > max_edges ) max_edges = nnodes;
    }
    // if in parallel, communicate the actual max_edges; it is not needed for tgt mesh (to be
    // global) but it is better to be consistent
#ifdef MOAB_HAVE_MPI
    if( parcomm )
    {
        int local_max_edges = max_edges;
        // now reduce max_edges over all processors
        int mpi_err =
            MPI_Allreduce( &local_max_edges, &max_edges, 1, MPI_INT, MPI_MAX, parcomm->proc_config().proc_comm() );
        if( MPI_SUCCESS != mpi_err ) return MB_FAILURE;
    }
#endif

    return MB_SUCCESS;
}
ErrorCode Intx2Mesh::FindMaxEdges( EntityHandle set1, EntityHandle set2 )
{
    ErrorCode rval = FindMaxEdgesInSet( set1, max_edges_1 );MB_CHK_SET_ERR( rval, "can't determine max_edges in set 1" );
    rval = FindMaxEdgesInSet( set2, max_edges_2 );MB_CHK_SET_ERR( rval, "can't determine max_edges in set 2" );

    return MB_SUCCESS;
}

ErrorCode Intx2Mesh::createTags()
{
    if( tgtParentTag ) mb->tag_delete( tgtParentTag );
    if( srcParentTag ) mb->tag_delete( srcParentTag );
    if( countTag ) mb->tag_delete( countTag );

    unsigned char def_data_bit = 0;  // unused by default
    // maybe the tgt tag is better to be deleted every time, and recreated;
    // or is it easy to set all values to something again? like 0?
    ErrorCode rval = mb->tag_get_handle( "tgtFlag", 1, MB_TYPE_BIT, TgtFlagTag, MB_TAG_CREAT, &def_data_bit );MB_CHK_SET_ERR( rval, "can't get tgt flag tag" );

    // create tgt edges if they do not exist yet; so when they are looked upon, they are found
    // this is the only call that is potentially NlogN, in the whole method
    rval = mb->get_adjacencies( rs2, 1, true, TgtEdges, Interface::UNION );MB_CHK_SET_ERR( rval, "can't get adjacent tgt edges" );

    // now, create a map from each edge to a list of potential new nodes on a tgt edge
    // this memory has to be cleaned up
    // change it to a vector, and use the index in range of tgt edges
    int indx = 0;
    extraNodesVec.resize( TgtEdges.size() );
    for( Range::iterator eit = TgtEdges.begin(); eit != TgtEdges.end(); ++eit, indx++ )
    {
        std::vector< EntityHandle >* nv = new std::vector< EntityHandle >;
        extraNodesVec[indx]             = nv;
    }

    int defaultInt = -1;

    rval = mb->tag_get_handle( "TargetParent", 1, MB_TYPE_INTEGER, tgtParentTag, MB_TAG_DENSE | MB_TAG_CREAT,
                               &defaultInt );MB_CHK_SET_ERR( rval, "can't create positive tag" );

    rval = mb->tag_get_handle( "SourceParent", 1, MB_TYPE_INTEGER, srcParentTag, MB_TAG_DENSE | MB_TAG_CREAT,
                               &defaultInt );MB_CHK_SET_ERR( rval, "can't create negative tag" );

    rval = mb->tag_get_handle( "Counting", 1, MB_TYPE_INTEGER, countTag, MB_TAG_DENSE | MB_TAG_CREAT, &defaultInt );MB_CHK_SET_ERR( rval, "can't create Counting tag" );

    // for each cell in set 1, determine its neigh in set 1 (could be null too)
    // for each cell in set 2, determine its neigh in set 2 (if on boundary, could be 0)
    rval = DetermineOrderedNeighbors( mbs1, max_edges_1, srcNeighTag );MB_CHK_SET_ERR( rval, "can't determine neighbors for set 1" );
    rval = DetermineOrderedNeighbors( mbs2, max_edges_2, tgtNeighTag );MB_CHK_SET_ERR( rval, "can't determine neighbors for set 2" );

    // for tgt cells, save a dense tag with the bordering edges, so we do not have to search for
    // them each time edges were for sure created before (tgtEdges)
    std::vector< EntityHandle > zeroh( max_edges_2, 0 );
    // if we have a tag with this name, it could be of a different size, so delete it
    rval = mb->tag_get_handle( "__tgtEdgeNeighbors", neighTgtEdgeTag );
    if( rval == MB_SUCCESS && neighTgtEdgeTag ) mb->tag_delete( neighTgtEdgeTag );
    rval = mb->tag_get_handle( "__tgtEdgeNeighbors", max_edges_2, MB_TYPE_HANDLE, neighTgtEdgeTag,
                               MB_TAG_DENSE | MB_TAG_CREAT, &zeroh[0] );MB_CHK_SET_ERR( rval, "can't create tgt edge neighbors tag" );
    for( Range::iterator rit = rs2.begin(); rit != rs2.end(); rit++ )
    {
        EntityHandle tgtCell = *rit;
        int num_nodes        = 0;
        rval                 = mb->get_connectivity( tgtCell, tgtConn, num_nodes );MB_CHK_SET_ERR( rval, "can't get  tgt conn" );
        // account for padded polygons
        while( tgtConn[num_nodes - 2] == tgtConn[num_nodes - 1] && num_nodes > 3 )
            num_nodes--;

        int i = 0;
        for( i = 0; i < num_nodes; i++ )
        {
            EntityHandle v[2] = { tgtConn[i],
                                  tgtConn[( i + 1 ) % num_nodes] };  // this is fine even for padded polygons
            std::vector< EntityHandle > adj_entities;
            rval = mb->get_adjacencies( v, 2, 1, false, adj_entities, Interface::INTERSECT );
            if( rval != MB_SUCCESS || adj_entities.size() < 1 ) return rval;  // get out , big error
            zeroh[i] = adj_entities[0];                                       // should be only one edge between 2 nodes
            // also, even if number of edges is less than max_edges_2, they will be ignored, even if
            // the tag is dense
        }
        // now set the value of the tag
        rval = mb->tag_set_data( neighTgtEdgeTag, &tgtCell, 1, &( zeroh[0] ) );MB_CHK_SET_ERR( rval, "can't set edge tgt tag" );
    }
    return MB_SUCCESS;
}

ErrorCode Intx2Mesh::DetermineOrderedNeighbors( EntityHandle inputSet, int max_edges, Tag& neighTag )
{
    Range cells;
    ErrorCode rval = mb->get_entities_by_dimension( inputSet, 2, cells );MB_CHK_SET_ERR( rval, "can't get cells in set" );

    std::vector< EntityHandle > neighbors( max_edges );
    std::vector< EntityHandle > zeroh( max_edges, 0 );
    // nameless tag, as the name is not important; we will have 2 related tags, but one on tgt mesh,
    // one on src mesh
    rval = mb->tag_get_handle( "", max_edges, MB_TYPE_HANDLE, neighTag, MB_TAG_DENSE | MB_TAG_CREAT, &zeroh[0] );MB_CHK_SET_ERR( rval, "can't create neighbors tag" );

    for( Range::iterator cit = cells.begin(); cit != cells.end(); cit++ )
    {
        EntityHandle cell = *cit;
        int nnodes        = 3;
        // will get the nnodes ordered neighbors;
        // first cell is for nodes 0, 1, second to 1, 2, third to 2, 3, last to nnodes-1,
        const EntityHandle* conn4;
        rval = mb->get_connectivity( cell, conn4, nnodes );MB_CHK_SET_ERR( rval, "can't get connectivity of a cell" );
        int nsides = nnodes;
        // account for possible padded polygons
        while( conn4[nsides - 2] == conn4[nsides - 1] && nsides > 3 )
            nsides--;

        for( int i = 0; i < nsides; i++ )
        {
            EntityHandle v[2];
            v[0] = conn4[i];
            v[1] = conn4[( i + 1 ) % nsides];
            // get all cells adjacent to these 2 vertices on the edge
            std::vector< EntityHandle > adjcells;
            std::vector< EntityHandle > cellsInSet;
            rval = mb->get_adjacencies( v, 2, 2, false, adjcells, Interface::INTERSECT );MB_CHK_SET_ERR( rval, "can't adjacency to 2 verts" );
            // now look for the cells contained in the input set;
            // the input set should be a correct mesh, not overlapping cells, and manifold
            size_t siz = adjcells.size();
            for( size_t j = 0; j < siz; j++ )
                if( mb->contains_entities( inputSet, &( adjcells[j] ), 1 ) ) cellsInSet.push_back( adjcells[j] );
            siz = cellsInSet.size();

            if( siz > 2 )
            {
                std::cout << "non manifold mesh, error" << mb->list_entities( &( cellsInSet[0] ), cellsInSet.size() )
                          << "\n";MB_CHK_SET_ERR( MB_FAILURE, "non-manifold input mesh set" );  // non-manifold
            }
            if( siz == 1 )
            {
                // it must be the border of the input mesh;
                neighbors[i] = 0;  // we are guaranteed that ids are !=0; this is marking a border
                // borders do not appear for a sphere in serial, but they do appear for
                // parallel processing anyway
                continue;
            }
            // here siz ==2, it is either the first or second
            if( cell == cellsInSet[0] )
                neighbors[i] = cellsInSet[1];
            else
                neighbors[i] = cellsInSet[0];
        }
        // fill the rest with 0
        for( int i = nsides; i < max_edges; i++ )
            neighbors[i] = 0;
        // now simply set the neighbors tag; the last few positions will not be used, but for
        // simplicity will keep them all (MAXEDGES)
        rval = mb->tag_set_data( neighTag, &cell, 1, &neighbors[0] );MB_CHK_SET_ERR( rval, "can't set neigh tag" );
    }
    return MB_SUCCESS;
}

// slow interface; this will not do the advancing front trick
// some are triangles, some are quads, some are polygons ...
ErrorCode Intx2Mesh::intersect_meshes_kdtree( EntityHandle mbset1, EntityHandle mbset2, EntityHandle& outputSet )
{
    ErrorCode rval;
    mbs1   = mbset1;  // set 1 is departure, and it is completely covering the euler set on proc
    mbs2   = mbset2;
    outSet = outputSet;
    rval   = mb->get_entities_by_dimension( mbs1, 2, rs1 );MB_CHK_ERR( rval );
    rval = mb->get_entities_by_dimension( mbs2, 2, rs2 );MB_CHK_ERR( rval );
    // from create tags, copy relevant ones
    if( tgtParentTag ) mb->tag_delete( tgtParentTag );
    if( srcParentTag ) mb->tag_delete( srcParentTag );
    if( countTag ) mb->tag_delete( countTag );

    // create tgt edges if they do not exist yet; so when they are looked upon, they are found
    // this is the only call that is potentially NlogN, in the whole method
    rval = mb->get_adjacencies( rs2, 1, true, TgtEdges, Interface::UNION );MB_CHK_SET_ERR( rval, "can't get adjacent tgt edges" );

    int indx = 0;
    extraNodesVec.resize( TgtEdges.size() );
    for( Range::iterator eit = TgtEdges.begin(); eit != TgtEdges.end(); ++eit, indx++ )
    {
        std::vector< EntityHandle >* nv = new std::vector< EntityHandle >;
        extraNodesVec[indx]             = nv;
    }

    int defaultInt = -1;

    rval = mb->tag_get_handle( "TargetParent", 1, MB_TYPE_INTEGER, tgtParentTag, MB_TAG_DENSE | MB_TAG_CREAT,
                               &defaultInt );MB_CHK_SET_ERR( rval, "can't create positive tag" );

    rval = mb->tag_get_handle( "SourceParent", 1, MB_TYPE_INTEGER, srcParentTag, MB_TAG_DENSE | MB_TAG_CREAT,
                               &defaultInt );MB_CHK_SET_ERR( rval, "can't create negative tag" );

    rval = mb->tag_get_handle( "Counting", 1, MB_TYPE_INTEGER, countTag, MB_TAG_DENSE | MB_TAG_CREAT, &defaultInt );MB_CHK_SET_ERR( rval, "can't create Counting tag" );

    // for tgt cells, save a dense tag with the bordering edges, so we do not have to search for
    // them each time edges were for sure created before (tgtEdges)
    // if we have a tag with this name, it could be of a different size, so delete it
    rval = mb->tag_get_handle( "__tgtEdgeNeighbors", neighTgtEdgeTag );
    if( rval == MB_SUCCESS && neighTgtEdgeTag ) mb->tag_delete( neighTgtEdgeTag );
    std::vector< EntityHandle > zeroh( max_edges_2, 0 );
    rval = mb->tag_get_handle( "__tgtEdgeNeighbors", max_edges_2, MB_TYPE_HANDLE, neighTgtEdgeTag,
                               MB_TAG_DENSE | MB_TAG_CREAT, &zeroh[0] );MB_CHK_SET_ERR( rval, "can't create tgt edge neighbors tag" );
    for( Range::iterator rit = rs2.begin(); rit != rs2.end(); rit++ )
    {
        EntityHandle tgtCell = *rit;
        int num_nodes        = 0;
        rval                 = mb->get_connectivity( tgtCell, tgtConn, num_nodes );MB_CHK_SET_ERR( rval, "can't get  tgt conn" );
        // account for padded polygons
        while( tgtConn[num_nodes - 2] == tgtConn[num_nodes - 1] && num_nodes > 3 )
            num_nodes--;

        int i = 0;
        for( i = 0; i < num_nodes; i++ )
        {
            EntityHandle v[2] = { tgtConn[i],
                                  tgtConn[( i + 1 ) % num_nodes] };  // this is fine even for padded polygons
            std::vector< EntityHandle > adj_entities;
            rval = mb->get_adjacencies( v, 2, 1, false, adj_entities, Interface::INTERSECT );
            if( rval != MB_SUCCESS || adj_entities.size() < 1 ) return rval;  // get out , big error
            zeroh[i] = adj_entities[0];                                       // should be only one edge between 2 nodes
            // also, even if number of edges is less than max_edges_2, they will be ignored, even if
            // the tag is dense
        }
        // now set the value of the tag
        rval = mb->tag_set_data( neighTgtEdgeTag, &tgtCell, 1, &( zeroh[0] ) );MB_CHK_SET_ERR( rval, "can't set edge tgt tag" );
    }

    // create the kd tree on source cells, and intersect all targets in an expensive loop
    // build a kd tree with the rs1 (source) cells
    AdaptiveKDTree kd( mb );
    EntityHandle tree_root = 0;
    rval                   = kd.build_tree( rs1, &tree_root );MB_CHK_ERR( rval );

    // find out max edge on source mesh;
    double max_length = 0;
    {
        std::vector< double > coords;
        coords.resize( 3 * max_edges_1 );
        for( Range::iterator it = rs1.begin(); it != rs1.end(); it++ )
        {
            const EntityHandle* conn = NULL;
            int nnodes;
            rval = mb->get_connectivity( *it, conn, nnodes );MB_CHK_SET_ERR( rval, "can't get connectivity" );
            while( conn[nnodes - 2] == conn[nnodes - 1] && nnodes > 3 )
                nnodes--;
            rval = mb->get_coords( conn, nnodes, &coords[0] );MB_CHK_SET_ERR( rval, "can't get coordinates" );
            for( int j = 0; j < nnodes; j++ )
            {
                int next = ( j + 1 ) % nnodes;
                double leng;
                leng = ( coords[3 * j] - coords[3 * next] ) * ( coords[3 * j] - coords[3 * next] ) +
                       ( coords[3 * j + 1] - coords[3 * next + 1] ) * ( coords[3 * j + 1] - coords[3 * next + 1] ) +
                       ( coords[3 * j + 2] - coords[3 * next + 2] ) * ( coords[3 * j + 2] - coords[3 * next + 2] );
                leng = sqrt( leng );
                if( leng > max_length ) max_length = leng;
            }
        }
    }
    // maximum sag on a spherical mesh make sense only for intx on a sphere, with radius 1 :(
    double tolerance = 1.e-15;
    if( max_length < 1. )
    {
        // basically, the sag for an arc of length max_length on a circle of radius 1
        tolerance = 1. - sqrt( 1 - max_length * max_length / 4 );
        if( box_error < tolerance ) box_error = tolerance;
        tolerance = 3 * tolerance;  // we use it for gnomonic plane too, projected sag could be =* sqrt(2.)
        // be more generous, use 1.5 ~= sqrt(2.)

        if( !my_rank )
        {
            std::cout << " max edge length: " << max_length << "  tolerance for kd tree: " << tolerance << "\n";
            std::cout << " box overlap tolerance: " << box_error << "\n";
        }
    }
    for( Range::iterator it = rs2.begin(); it != rs2.end(); ++it )
    {
        EntityHandle tcell = *it;
        // find vertex positions
        const EntityHandle* conn = NULL;
        int nnodes               = 0;
        rval                     = mb->get_connectivity( tcell, conn, nnodes );MB_CHK_ERR( rval );
        // find leaves close to those positions
        double areaTgtCell   = setup_tgt_cell( tcell, nnodes );  // this is the area in the gnomonic plane
        double recoveredArea = 0;
        std::vector< double > positions;
        positions.resize( nnodes * 3 );
        rval = mb->get_coords( conn, nnodes, &positions[0] );MB_CHK_ERR( rval );

        // distance to search will be based on average edge length
        double av_len = 0;
        for( int k = 0; k < nnodes; k++ )
        {
            int ik      = ( k + 1 ) % nnodes;
            double len1 = 0;
            for( int j = 0; j < 3; j++ )
            {
                double len2 = positions[3 * k + j] - positions[3 * ik + j];
                len1 += len2 * len2;
            }
            av_len += sqrt( len1 );
        }
        if( nnodes > 0 ) av_len /= nnodes;
        // find leaves within a distance from each vertex of target
        // in those leaves, collect all cells; we will try for an intx in there
        Range close_source_cells;
        std::vector< EntityHandle > leaves;
        for( int i = 0; i < nnodes; i++ )
        {

            leaves.clear();
            rval = kd.distance_search( &positions[3 * i], av_len, leaves, tolerance, epsilon_1 );MB_CHK_ERR( rval );

            for( std::vector< EntityHandle >::iterator j = leaves.begin(); j != leaves.end(); ++j )
            {
                Range tmp;
                rval = mb->get_entities_by_dimension( *j, 2, tmp );MB_CHK_ERR( rval );

                close_source_cells.merge( tmp.begin(), tmp.end() );
            }
        }
#ifdef VERBOSE
        if( close_source_cells.empty() )
        {
            std::cout << " there are no close source cells to target cell " << tcell << " id from handle "
                      << mb->id_from_handle( tcell ) << "\n";
        }
#endif
        for( Range::iterator it2 = close_source_cells.begin(); it2 != close_source_cells.end(); ++it2 )
        {
            EntityHandle startSrc = *it2;
            double area           = 0;
            // if area is > 0 , we have intersections
            double P[10 * MAXEDGES];  // max 8 intx points + 8 more in the polygon
            //
            int nP = 0;
            int nb[MAXEDGES], nr[MAXEDGES];  // sides 3 or 4? also, check boxes first
            int nsTgt, nsSrc;
            rval = computeIntersectionBetweenTgtAndSrc( tcell, startSrc, P, nP, area, nb, nr, nsSrc, nsTgt, true );MB_CHK_ERR( rval );
            if( area > 0 )
            {
                if( nP > 1 )
                {  // this will also construct triangles/polygons in the new mesh, if needed
                    rval = findNodes( tcell, nnodes, startSrc, nsSrc, P, nP );MB_CHK_ERR( rval );
                }
                recoveredArea += area;
            }
        }
        recoveredArea = ( recoveredArea - areaTgtCell ) / areaTgtCell;  // replace now with recovery fract
    }
    // before cleaning up , we need to settle the position of the intersection points
    // on the boundary edges
    // this needs to be collective, so we should maybe wait something
#ifdef MOAB_HAVE_MPI
    rval = resolve_intersection_sharing();MB_CHK_SET_ERR( rval, "can't correct position, Intx2Mesh.cpp \n" );
#endif

    this->clean();
    return MB_SUCCESS;
}
// main interface; this will do the advancing front trick
// some are triangles, some are quads, some are polygons ...
ErrorCode Intx2Mesh::intersect_meshes( EntityHandle mbset1, EntityHandle mbset2, EntityHandle& outputSet )
{
    ErrorCode rval;
    mbs1   = mbset1;  // set 1 is departure, and it is completely covering the euler set on proc
    mbs2   = mbset2;
    outSet = outputSet;
#ifdef VERBOSE
    std::stringstream ffs, fft;
    ffs << "source_rank0" << my_rank << ".vtk";
    rval = mb->write_mesh( ffs.str().c_str(), &mbset1, 1 );MB_CHK_ERR( rval );
    fft << "target_rank0" << my_rank << ".vtk";
    rval = mb->write_mesh( fft.str().c_str(), &mbset2, 1 );MB_CHK_ERR( rval );

#endif
    // really, should be something from t1 and t2; src is 1 (lagrange), tgt is 2 (euler)

    EntityHandle startSrc = 0, startTgt = 0;

    rval = mb->get_entities_by_dimension( mbs1, 2, rs1 );MB_CHK_ERR( rval );
    rval = mb->get_entities_by_dimension( mbs2, 2, rs2 );MB_CHK_ERR( rval );
    // std::cout << "rs1.size() = " << rs1.size() << " and rs2.size() = "  << rs2.size() << "\n";
    // std::cout.flush();

    createTags();  // will also determine max_edges_1, max_edges_2 (for src and tgt meshes)

    Range rs22 = rs2;  // a copy of the initial range; we will remove from it elements as we
                       // advance ; rs2 is needed for marking the polygon to the tgt parent

    // create the local kdd tree with source elements; will use it to search
    // more efficiently for the seeds in advancing front;
    // some of the target cells will not be covered by source cells, and they need to be eliminated
    // early from contention

    // build a kd tree with the rs1 (source) cells
    AdaptiveKDTree kd( mb );
    EntityHandle tree_root = 0;
    rval                   = kd.build_tree( rs1, &tree_root );MB_CHK_ERR( rval );

    while( !rs22.empty() )
    {
#if defined( ENABLE_DEBUG ) || defined( VERBOSE )
        if( rs22.size() < rs2.size() )
        {
            std::cout << " possible not connected arrival mesh; my_rank: " << my_rank << " counting: " << counting
                      << "\n";
            std::stringstream ffo;
            ffo << "file0" << counting << "rank0" << my_rank << ".vtk";
            rval = mb->write_mesh( ffo.str().c_str(), &outSet, 1 );MB_CHK_ERR( rval );
        }
#endif
        bool seedFound = false;
        for( Range::iterator it = rs22.begin(); it != rs22.end(); ++it )
        {
            startTgt  = *it;
            int found = 0;
            // find vertex positions
            const EntityHandle* conn = NULL;
            int nnodes               = 0;
            rval                     = mb->get_connectivity( startTgt, conn, nnodes );MB_CHK_ERR( rval );
            // find leaves close to those positions
            std::vector< double > positions;
            positions.resize( nnodes * 3 );
            rval = mb->get_coords( conn, nnodes, &positions[0] );MB_CHK_ERR( rval );
            // find leaves within a distance from each vertex of target
            // in those leaves, collect all cells; we will try for an intx in there, instead of
            // looping over all rs1 cells, as before
            Range close_source_cells;
            std::vector< EntityHandle > leaves;
            for( int i = 0; i < nnodes; i++ )
            {

                leaves.clear();
                rval = kd.distance_search( &positions[3 * i], epsilon_1, leaves, epsilon_1, epsilon_1 );MB_CHK_ERR( rval );

                for( std::vector< EntityHandle >::iterator j = leaves.begin(); j != leaves.end(); ++j )
                {
                    Range tmp;
                    rval = mb->get_entities_by_dimension( *j, 2, tmp );MB_CHK_ERR( rval );

                    close_source_cells.merge( tmp.begin(), tmp.end() );
                }
            }

            for( Range::iterator it2 = close_source_cells.begin(); it2 != close_source_cells.end() && !found; ++it2 )
            {
                startSrc    = *it2;
                double area = 0;
                // if area is > 0 , we have intersections
                double P[10 * MAXEDGES];  // max 8 intx points + 8 more in the polygon
                //
                int nP = 0;
                int nb[MAXEDGES], nr[MAXEDGES];  // sides 3 or 4? also, check boxes first
                int nsTgt, nsSrc;
                rval =
                    computeIntersectionBetweenTgtAndSrc( startTgt, startSrc, P, nP, area, nb, nr, nsSrc, nsTgt, true );MB_CHK_ERR( rval );
                if( area > 0 )
                {
                    found     = 1;
                    seedFound = true;
                    break;  // found 2 elements that intersect; these will be the seeds
                }
            }
            if( found )
                break;
            else
            {
#if defined( VERBOSE )
                std::cout << " on rank " << my_rank << " target cell " << ID_FROM_HANDLE( startTgt )
                          << " not intx with any source\n";
#endif
                rs22.erase( startTgt );
            }
        }
        if( !seedFound ) continue;  // continue while(!rs22.empty())

        std::queue< EntityHandle > srcQueue;  // these are corresponding to Ta,
        srcQueue.push( startSrc );
        std::queue< EntityHandle > tgtQueue;
        tgtQueue.push( startTgt );

        Range toResetSrcs;  // will be used to reset src flags for every tgt quad
        // processed

        /*if (my_rank==0)
          dbg_1 = 1;*/
        unsigned char used = 1;
        // mark the start tgt quad as used, so it will not come back again
        rval = mb->tag_set_data( TgtFlagTag, &startTgt, 1, &used );MB_CHK_ERR( rval );
        while( !tgtQueue.empty() )
        {
            // flags for the side : 0 means a src cell not found on side
            // a paired src not found yet for the neighbors of tgt
            Range nextSrc[MAXEDGES];  // there are new ranges of possible next src cells for
                                      // seeding the side j of tgt cell

            EntityHandle currentTgt = tgtQueue.front();
            tgtQueue.pop();
            int nsidesTgt;                                                   // will be initialized now
            double areaTgtCell   = setup_tgt_cell( currentTgt, nsidesTgt );  // this is the area in the gnomonic plane
            double recoveredArea = 0;
            // get the neighbors of tgt, and if they are solved already, do not bother with that
            // side of tgt
            EntityHandle tgtNeighbors[MAXEDGES] = { 0 };
            rval                                = mb->tag_get_data( tgtNeighTag, &currentTgt, 1, tgtNeighbors );MB_CHK_SET_ERR( rval, "can't get neighbors of current tgt" );
#ifdef ENABLE_DEBUG
            if( dbg_1 )
            {
                std::cout << "Next: neighbors for current tgt ";
                for( int kk = 0; kk < nsidesTgt; kk++ )
                {
                    if( tgtNeighbors[kk] > 0 )
                        std::cout << mb->id_from_handle( tgtNeighbors[kk] ) << " ";
                    else
                        std::cout << 0 << " ";
                }
                std::cout << std::endl;
            }
#endif
            // now get the status of neighbors; if already solved, make them 0, so not to bother
            // anymore on that side of tgt
            for( int j = 0; j < nsidesTgt; j++ )
            {
                EntityHandle tgtNeigh = tgtNeighbors[j];
                unsigned char status  = 1;
                if( tgtNeigh == 0 ) continue;
                rval = mb->tag_get_data( TgtFlagTag, &tgtNeigh, 1, &status );MB_CHK_ERR( rval );                     // status 0 is unused
                if( 1 == status ) tgtNeighbors[j] = 0;  // so will not look anymore on this side of tgt
            }

#ifdef ENABLE_DEBUG
            if( dbg_1 )
            {
                std::cout << "reset sources: ";
                for( Range::iterator itr = toResetSrcs.begin(); itr != toResetSrcs.end(); ++itr )
                    std::cout << mb->id_from_handle( *itr ) << " ";
                std::cout << std::endl;
            }
#endif
            EntityHandle currentSrc = srcQueue.front();
            // tgt and src queues are parallel; for clarity we should have kept in the queue pairs
            // of entity handle std::pair<EntityHandle, EntityHandle>; so just one queue, with
            // pairs;
            //  at every moment, the queue contains pairs of cells that intersect, and they form the
            //  "advancing front"
            srcQueue.pop();
            toResetSrcs.clear();  // empty the range of used srcs, will have to be set unused again,
            // at the end of tgt element processing
            toResetSrcs.insert( currentSrc );
            // mb2->set_tag_data
            std::queue< EntityHandle > localSrc;
            localSrc.push( currentSrc );
#ifdef VERBOSE
            int countingStart = counting;
#endif
            // will advance-front search in the neighborhood of tgt cell, until we finish processing
            // all
            //   possible src cells; localSrc queue will contain all possible src cells that cover
            //   the current tgt cell
            while( !localSrc.empty() )
            {
                //
                EntityHandle srcT = localSrc.front();
                localSrc.pop();
                double P[10 * MAXEDGES], area;  //
                int nP           = 0;
                int nb[MAXEDGES] = { 0 };
                int nr[MAXEDGES] = { 0 };

                int nsidesSrc;  ///
                // area is in 2d, points are in 3d (on a sphere), back-projected, or in a plane
                // intersection points could include the vertices of initial elements
                // nb [j] = 0 means no intersection on the side j for element src (markers)
                // nb [j] = 1 means that the side j (from j to j+1) of src poly intersects the
                // tgt poly.  A potential next poly in the tgt queue is the tgt poly that is
                // adjacent to this side
                rval = computeIntersectionBetweenTgtAndSrc( /* tgt */ currentTgt, srcT, P, nP, area, nb, nr, nsidesSrc,
                                                            nsidesTgt );MB_CHK_ERR( rval );
                if( nP > 0 )
                {
#ifdef ENABLE_DEBUG
                    if( dbg_1 )
                    {
                        for( int k = 0; k < 3; k++ )
                        {
                            std::cout << " nb, nr: " << k << " " << nb[k] << " " << nr[k] << "\n";
                        }
                    }
#endif

                    // intersection found: output P and original triangles if nP > 2
                    EntityHandle neighbors[MAXEDGES] = { 0 };
                    rval                             = mb->tag_get_data( srcNeighTag, &srcT, 1, neighbors );
                    if( rval != MB_SUCCESS )
                    {
                        std::cout << " can't get the neighbors for src element " << mb->id_from_handle( srcT );
                        return MB_FAILURE;
                    }

                    // add neighbors to the localSrc queue, if they are not marked
                    for( int nn = 0; nn < nsidesSrc; nn++ )
                    {
                        EntityHandle neighbor = neighbors[nn];
                        if( neighbor > 0 && nb[nn] > 0 )  // advance across src boundary nn
                        {
                            if( toResetSrcs.find( neighbor ) == toResetSrcs.end() )
                            {
                                localSrc.push( neighbor );
#ifdef ENABLE_DEBUG
                                if( dbg_1 )
                                {
                                    std::cout << " local src elem " << mb->id_from_handle( neighbor )
                                              << " for tgt:" << mb->id_from_handle( currentTgt ) << "\n";
                                    mb->list_entities( &neighbor, 1 );
                                }
#endif
                                toResetSrcs.insert( neighbor );
                            }
                        }
                    }
                    // n(find(nc>0))=ac;        % ac is starting candidate for neighbor
                    for( int nn = 0; nn < nsidesTgt; nn++ )
                    {
                        if( nr[nn] > 0 && tgtNeighbors[nn] > 0 )
                            nextSrc[nn].insert( srcT );  // potential src cell that can intersect
                                                         // the tgt neighbor nn
                    }
                    if( nP > 1 )
                    {  // this will also construct triangles/polygons in the new mesh, if needed
                        rval = findNodes( currentTgt, nsidesTgt, srcT, nsidesSrc, P, nP );MB_CHK_ERR( rval );
                    }

                    recoveredArea += area;
                }
#ifdef ENABLE_DEBUG
                else if( dbg_1 )
                {
                    std::cout << " tgt, src, do not intersect: " << mb->id_from_handle( currentTgt ) << " "
                              << mb->id_from_handle( srcT ) << "\n";
                }
#endif
            }                                                               // end while (!localSrc.empty())
            recoveredArea = ( recoveredArea - areaTgtCell ) / areaTgtCell;  // replace now with recovery fraction
#if defined( ENABLE_DEBUG ) || defined( VERBOSE )
            if( fabs( recoveredArea ) > epsilon_1 )
            {
#ifdef VERBOSE
                std::cout << " tgt area: " << areaTgtCell << " recovered :" << recoveredArea * ( 1 + areaTgtCell )
                          << " fraction error recovery:" << recoveredArea
                          << " tgtID: " << mb->id_from_handle( currentTgt ) << " countingStart:" << countingStart
                          << "\n";
#endif
            }
#endif
            // here, we are finished with tgtCurrent, take it out of the rs22 range (tgt, arrival
            // mesh)
            rs22.erase( currentTgt );
            // also, look at its neighbors, and add to the seeds a next one

            for( int j = 0; j < nsidesTgt; j++ )
            {
                EntityHandle tgtNeigh = tgtNeighbors[j];
                if( tgtNeigh == 0 || nextSrc[j].size() == 0 )  // if tgt is bigger than src, there could be no src
                                                               // to advance on that side
                    continue;
                int nsidesTgt2 = 0;
                setup_tgt_cell( tgtNeigh,
                                nsidesTgt2 );  // find possible intersection with src cell from nextSrc
                for( Range::iterator nit = nextSrc[j].begin(); nit != nextSrc[j].end(); ++nit )
                {
                    EntityHandle nextB = *nit;
                    // we identified tgt quad n[j] as possibly intersecting with neighbor j of the
                    // src quad
                    double P[10 * MAXEDGES], area;  //
                    int nP           = 0;
                    int nb[MAXEDGES] = { 0 };
                    int nr[MAXEDGES] = { 0 };

                    int nsidesSrc;  ///
                    rval = computeIntersectionBetweenTgtAndSrc(
                        /* tgt */ tgtNeigh, nextB, P, nP, area, nb, nr, nsidesSrc, nsidesTgt2 );MB_CHK_ERR( rval );
                    if( area > 0 )
                    {
                        tgtQueue.push( tgtNeigh );
                        srcQueue.push( nextB );
#ifdef ENABLE_DEBUG
                        if( dbg_1 )
                            std::cout << "new polys pushed: src, tgt:" << mb->id_from_handle( tgtNeigh ) << " "
                                      << mb->id_from_handle( nextB ) << std::endl;
#endif
                        rval = mb->tag_set_data( TgtFlagTag, &tgtNeigh, 1, &used );MB_CHK_ERR( rval );
                        break;  // so we are done with this side of tgt, we have found a proper next
                                // seed
                    }
                }
            }

        }  // end while (!tgtQueue.empty())
    }
#ifdef ENABLE_DEBUG
    if( dbg_1 )
    {
        for( int k = 0; k < 6; k++ )
            mout_1[k].close();
    }
#endif
    // before cleaning up , we need to settle the position of the intersection points
    // on the boundary edges
    // this needs to be collective, so we should maybe wait something
#ifdef MOAB_HAVE_MPI
    rval = resolve_intersection_sharing();MB_CHK_SET_ERR( rval, "can't correct position, Intx2Mesh.cpp \n" );
#endif

    this->clean();
    return MB_SUCCESS;
}

// clean some memory allocated
void Intx2Mesh::clean()
{
    //
    int indx = 0;
    for( Range::iterator eit = TgtEdges.begin(); eit != TgtEdges.end(); ++eit, indx++ )
    {
        delete extraNodesVec[indx];
    }
    // extraNodesMap.clear();
    extraNodesVec.clear();
    // also, delete some bit tags, used to mark processed tgts and srcs
    mb->tag_delete( TgtFlagTag );
    counting = 0;  // reset counting to original value
}
// this method will reduce number of nodes, collapse edges that are of length 0
// so a polygon like 428 431 431 will become a line 428 431
// or something like 428 431 431 531 -> 428 431 531
void Intx2Mesh::correct_polygon( EntityHandle* nodes, int& nP )
{
    int i = 0;
    while( i < nP )
    {
        int nextIndex = ( i + 1 ) % nP;
        if( nodes[i] == nodes[nextIndex] )
        {
#ifdef ENABLE_DEBUG
            // we need to reduce nP, and collapse nodes
            if( dbg_1 )
            {
                std::cout << " nodes duplicated in list: ";
                for( int j = 0; j < nP; j++ )
                    std::cout << nodes[j] << " ";
                std::cout << "\n";
                std::cout << " node " << nodes[i] << " at index " << i << " is duplicated"
                          << "\n";
            }
#endif
            // this will work even if we start from 1 2 3 1; when i is 3, we find nextIndex is 0,
            // then next thing does nothing
            //  (nP-1 is 3, so k is already >= nP-1); it will result in nodes -> 1, 2, 3
            for( int k = i; k < nP - 1; k++ )
                nodes[k] = nodes[k + 1];
            nP--;  // decrease the number of nodes; also, decrease i, just if we may need to check
                   // again
            i--;
        }
        i++;
    }
    return;
}
#ifdef MOAB_HAVE_MPI
ErrorCode Intx2Mesh::build_processor_euler_boxes( EntityHandle euler_set, Range& local_verts )
{
    localEnts.clear();
    ErrorCode rval = mb->get_entities_by_dimension( euler_set, 2, localEnts );ERRORR( rval, "can't get ents by dimension" );

    rval                = mb->get_connectivity( localEnts, local_verts );
    int num_local_verts = (int)local_verts.size();ERRORR( rval, "can't get local vertices" );

    assert( parcomm != NULL );

    // get the position of local vertices, and decide local boxes (allBoxes...)
    double bmin[3] = { DBL_MAX, DBL_MAX, DBL_MAX };
    double bmax[3] = { -DBL_MAX, -DBL_MAX, -DBL_MAX };

    std::vector< double > coords( 3 * num_local_verts );
    rval = mb->get_coords( local_verts, &coords[0] );ERRORR( rval, "can't get coords of vertices " );

    for( int i = 0; i < num_local_verts; i++ )
    {
        for( int k = 0; k < 3; k++ )
        {
            double val = coords[3 * i + k];
            if( val < bmin[k] ) bmin[k] = val;
            if( val > bmax[k] ) bmax[k] = val;
        }
    }
    int numprocs = parcomm->proc_config().proc_size();
    allBoxes.resize( 6 * numprocs );

    my_rank = parcomm->proc_config().proc_rank();
    for( int k = 0; k < 3; k++ )
    {
        allBoxes[6 * my_rank + k]     = bmin[k];
        allBoxes[6 * my_rank + 3 + k] = bmax[k];
    }

    // now communicate to get all boxes
    int mpi_err;
#if( MPI_VERSION >= 2 )
    // use "in place" option
    mpi_err = MPI_Allgather( MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, &allBoxes[0], 6, MPI_DOUBLE,
                             parcomm->proc_config().proc_comm() );
#else
    {
        std::vector< double > allBoxes_tmp( 6 * parcomm->proc_config().proc_size() );
        mpi_err  = MPI_Allgather( &allBoxes[6 * my_rank], 6, MPI_DOUBLE, &allBoxes_tmp[0], 6, MPI_DOUBLE,
                                  parcomm->proc_config().proc_comm() );
        allBoxes = allBoxes_tmp;
    }
#endif
    if( MPI_SUCCESS != mpi_err ) return MB_FAILURE;

    if( my_rank == 0 )
    {
        std::cout << " maximum number of vertices per cell are " << max_edges_1 << " on first mesh and " << max_edges_2
                  << " on second mesh \n";
        for( int i = 0; i < numprocs; i++ )
        {
            std::cout << "proc: " << i << " box min: " << allBoxes[6 * i] << " " << allBoxes[6 * i + 1] << " "
                      << allBoxes[6 * i + 2] << " \n";
            std::cout << "        box max: " << allBoxes[6 * i + 3] << " " << allBoxes[6 * i + 4] << " "
                      << allBoxes[6 * i + 5] << " \n";
        }
    }

    return MB_SUCCESS;
}
ErrorCode Intx2Mesh::create_departure_mesh_2nd_alg( EntityHandle& euler_set, EntityHandle& covering_lagr_set )
{
    // compute the bounding box on each proc
    assert( parcomm != NULL );

    localEnts.clear();
    ErrorCode rval = mb->get_entities_by_dimension( euler_set, 2, localEnts );ERRORR( rval, "can't get ents by dimension" );

    Tag dpTag = 0;
    std::string tag_name( "DP" );
    rval = mb->tag_get_handle( tag_name.c_str(), 3, MB_TYPE_DOUBLE, dpTag, MB_TAG_DENSE );ERRORR( rval, "can't get DP tag" );

    EntityHandle dum = 0;
    Tag corrTag;
    rval = mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE | MB_TAG_CREAT, &dum );ERRORR( rval, "can't get CORR tag" );
    // get all local verts
    Range local_verts;
    rval                = mb->get_connectivity( localEnts, local_verts );
    int num_local_verts = (int)local_verts.size();ERRORR( rval, "can't get local vertices" );

    rval = Intx2Mesh::build_processor_euler_boxes( euler_set, local_verts );ERRORR( rval, "can't build processor boxes" );

    std::vector< int > gids( num_local_verts );
    rval = mb->tag_get_data( gid, local_verts, &gids[0] );ERRORR( rval, "can't get local vertices gids" );

    // now see the departure points; to what boxes should we send them?
    std::vector< double > dep_points( 3 * num_local_verts );
    rval = mb->tag_get_data( dpTag, local_verts, (void*)&dep_points[0] );ERRORR( rval, "can't get DP tag values" );
    // ranges to send to each processor; will hold vertices and elements (quads?)
    // will look if the box of the dep quad covers box of euler mesh on proc (with tolerances)
    std::map< int, Range > Rto;
    int numprocs = parcomm->proc_config().proc_size();

    for( Range::iterator eit = localEnts.begin(); eit != localEnts.end(); ++eit )
    {
        EntityHandle q = *eit;
        const EntityHandle* conn4;
        int num_nodes;
        rval = mb->get_connectivity( q, conn4, num_nodes );ERRORR( rval, "can't get DP tag values" );
        CartVect qbmin( DBL_MAX );
        CartVect qbmax( -DBL_MAX );
        for( int i = 0; i < num_nodes; i++ )
        {
            EntityHandle v = conn4[i];
            size_t index   = local_verts.find( v ) - local_verts.begin();
            CartVect dp( &dep_points[3 * index] );  // will use constructor
            for( int j = 0; j < 3; j++ )
            {
                if( qbmin[j] > dp[j] ) qbmin[j] = dp[j];
                if( qbmax[j] < dp[j] ) qbmax[j] = dp[j];
            }
        }
        for( int p = 0; p < numprocs; p++ )
        {
            CartVect bbmin( &allBoxes[6 * p] );
            CartVect bbmax( &allBoxes[6 * p + 3] );
            if( GeomUtil::boxes_overlap( bbmin, bbmax, qbmin, qbmax, box_error ) )
            {
                Rto[p].insert( q );
            }
        }
    }

    // now, build TLv and TLq, for each p
    size_t numq = 0;
    size_t numv = 0;
    for( int p = 0; p < numprocs; p++ )
    {
        if( p == (int)my_rank ) continue;  // do not "send" it, because it is already here
        Range& range_to_P = Rto[p];
        // add the vertices to it
        if( range_to_P.empty() ) continue;  // nothing to send to proc p
        Range vertsToP;
        rval = mb->get_connectivity( range_to_P, vertsToP );ERRORR( rval, "can't get connectivity" );
        numq = numq + range_to_P.size();
        numv = numv + vertsToP.size();
        range_to_P.merge( vertsToP );
    }
    TupleList TLv;
    TupleList TLq;
    TLv.initialize( 2, 0, 0, 3, numv );  // to proc, GLOBAL ID, DP points
    TLv.enableWriteAccess();

    int sizeTuple = 2 + max_edges_1;  // determined earlier, for src, first mesh
    TLq.initialize( 2 + max_edges_1, 0, 1, 0,
                    numq );  // to proc, elem GLOBAL ID, connectivity[10] (global ID v), local eh
    TLq.enableWriteAccess();
#ifdef VERBOSE
    std::cout << "from proc " << my_rank << " send " << numv << " vertices and " << numq << " elements\n";
#endif
    for( int to_proc = 0; to_proc < numprocs; to_proc++ )
    {
        if( to_proc == (int)my_rank ) continue;
        Range& range_to_P = Rto[to_proc];
        Range V           = range_to_P.subset_by_type( MBVERTEX );

        for( Range::iterator it = V.begin(); it != V.end(); ++it )
        {
            EntityHandle v       = *it;
            unsigned int index   = local_verts.find( v ) - local_verts.begin();
            int n                = TLv.get_n();
            TLv.vi_wr[2 * n]     = to_proc;                // send to processor
            TLv.vi_wr[2 * n + 1] = gids[index];            // global id needs index in the local_verts range
            TLv.vr_wr[3 * n]     = dep_points[3 * index];  // departure position, of the node local_verts[i]
            TLv.vr_wr[3 * n + 1] = dep_points[3 * index + 1];
            TLv.vr_wr[3 * n + 2] = dep_points[3 * index + 2];
            TLv.inc_n();
        }
        // also, prep the quad for sending ...
        Range Q = range_to_P.subset_by_dimension( 2 );
        for( Range::iterator it = Q.begin(); it != Q.end(); ++it )
        {
            EntityHandle q = *it;
            int global_id;
            rval = mb->tag_get_data( gid, &q, 1, &global_id );ERRORR( rval, "can't get gid for polygon" );
            int n                        = TLq.get_n();
            TLq.vi_wr[sizeTuple * n]     = to_proc;    //
            TLq.vi_wr[sizeTuple * n + 1] = global_id;  // global id of element, used to identify it ...
            const EntityHandle* conn4;
            int num_nodes;
            rval = mb->get_connectivity( q, conn4,
                                         num_nodes );  // could be up to MAXEDGES, but it is limited by max_edges_1
            ERRORR( rval, "can't get connectivity for cell" );
            if( num_nodes > MAXEDGES ) ERRORR( MB_FAILURE, "too many nodes in a polygon" );
            for( int i = 0; i < num_nodes; i++ )
            {
                EntityHandle v                   = conn4[i];
                unsigned int index               = local_verts.find( v ) - local_verts.begin();
                TLq.vi_wr[sizeTuple * n + 2 + i] = gids[index];
            }
            for( int k = num_nodes; k < max_edges_1; k++ )
            {
                TLq.vi_wr[sizeTuple * n + 2 + k] =
                    0;  // fill the rest of node ids with 0; we know that the node ids start from 1!
            }
            TLq.vul_wr[n] = q;  // save here the entity handle, it will be communicated back
            // maybe we should forget about global ID
            TLq.inc_n();
        }
    }

    // now we are done populating the tuples; route them to the appropriate processors
    ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, TLv, 0 );
    ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, TLq, 0 );
    // the elements are already in localEnts;

    // maps from global ids to new vertex and quad handles, that are added
    std::map< int, EntityHandle > globalID_to_handle;
    /*std::map<int, EntityHandle> globalID_to_eh;*/
    globalID_to_eh.clear();  // need for next iteration
    // now, look at every TLv, and see if we have to create a vertex there or not
    int n = TLv.get_n();  // the size of the points received
    for( int i = 0; i < n; i++ )
    {
        int globalId = TLv.vi_rd[2 * i + 1];
        if( globalID_to_handle.find( globalId ) == globalID_to_handle.end() )
        {
            EntityHandle new_vert;
            double dp_pos[3] = { TLv.vr_wr[3 * i], TLv.vr_wr[3 * i + 1], TLv.vr_wr[3 * i + 2] };
            rval             = mb->create_vertex( dp_pos, new_vert );ERRORR( rval, "can't create new vertex " );
            globalID_to_handle[globalId] = new_vert;
        }
    }

    // now, all dep points should be at their place
    // look in the local list of q for this proc, and create all those quads and vertices if needed
    // it may be an overkill, but because it does not involve communication, we do it anyway
    Range& local  = Rto[my_rank];
    Range local_q = local.subset_by_dimension( 2 );
    // the local should have all the vertices in local_verts
    for( Range::iterator it = local_q.begin(); it != local_q.end(); ++it )
    {
        EntityHandle q = *it;
        int nnodes;
        const EntityHandle* conn4;
        rval = mb->get_connectivity( q, conn4, nnodes );ERRORR( rval, "can't get connectivity of local q " );
        EntityHandle new_conn[MAXEDGES];
        for( int i = 0; i < nnodes; i++ )
        {
            EntityHandle v1    = conn4[i];
            unsigned int index = local_verts.find( v1 ) - local_verts.begin();
            int globalId       = gids[index];
            if( globalID_to_handle.find( globalId ) == globalID_to_handle.end() )
            {
                // we need to create that vertex, at this position dep_points
                double dp_pos[3] = { dep_points[3 * index], dep_points[3 * index + 1], dep_points[3 * index + 2] };
                EntityHandle new_vert;
                rval = mb->create_vertex( dp_pos, new_vert );ERRORR( rval, "can't create new vertex " );
                globalID_to_handle[globalId] = new_vert;
            }
            new_conn[i] = globalID_to_handle[gids[index]];
        }
        EntityHandle new_element;
        //
        EntityType entType = MBQUAD;
        if( nnodes > 4 ) entType = MBPOLYGON;
        if( nnodes < 4 ) entType = MBTRI;

        rval = mb->create_element( entType, new_conn, nnodes, new_element );ERRORR( rval, "can't create new quad " );
        rval = mb->add_entities( covering_lagr_set, &new_element, 1 );ERRORR( rval, "can't add new element to dep set" );
        int gid_el;
        // get the global ID of the initial quad
        rval = mb->tag_get_data( gid, &q, 1, &gid_el );ERRORR( rval, "can't get element global ID " );
        globalID_to_eh[gid_el] = new_element;
        // is this redundant or not?
        rval = mb->tag_set_data( corrTag, &new_element, 1, &q );ERRORR( rval, "can't set corr tag on new el" );
        // set the global id on new elem
        rval = mb->tag_set_data( gid, &new_element, 1, &gid_el );ERRORR( rval, "can't set global id tag on new el" );
    }
    // now look at all elements received through; we do not want to duplicate them
    n = TLq.get_n();  // number of elements received by this processor
    // form the remote cells, that will be used to send the tracer info back to the originating proc
    remote_cells = new TupleList();
    remote_cells->initialize( 2, 0, 1, 0, n );  // will not have tracer data anymore
    remote_cells->enableWriteAccess();
    for( int i = 0; i < n; i++ )
    {
        int globalIdEl = TLq.vi_rd[sizeTuple * i + 1];
        int from_proc  = TLq.vi_wr[sizeTuple * i];
        // do we already have a quad with this global ID, represented?
        if( globalID_to_eh.find( globalIdEl ) == globalID_to_eh.end() )
        {
            // construct the conn quad
            EntityHandle new_conn[MAXEDGES];
            int nnodes = -1;
            for( int j = 0; j < max_edges_1; j++ )
            {
                int vgid = TLq.vi_rd[sizeTuple * i + 2 + j];  // vertex global ID
                if( vgid == 0 )
                    new_conn[j] = 0;
                else
                {
                    assert( globalID_to_handle.find( vgid ) != globalID_to_handle.end() );
                    new_conn[j] = globalID_to_handle[vgid];
                    nnodes      = j + 1;  // nodes are at the beginning, and are variable number
                }
            }
            EntityHandle new_element;
            //
            EntityType entType = MBQUAD;
            if( nnodes > 4 ) entType = MBPOLYGON;
            if( nnodes < 4 ) entType = MBTRI;
            rval = mb->create_element( entType, new_conn, nnodes, new_element );ERRORR( rval, "can't create new element " );
            globalID_to_eh[globalIdEl] = new_element;
            rval                       = mb->add_entities( covering_lagr_set, &new_element, 1 );ERRORR( rval, "can't add new element to dep set" );
            /* rval = mb->tag_set_data(corrTag, &new_element, 1, &q);ERRORR(rval, "can't set corr tag on new el");*/
            remote_cells->vi_wr[2 * i]     = from_proc;
            remote_cells->vi_wr[2 * i + 1] = globalIdEl;
            //     remote_cells->vr_wr[i] = 0.; // no contribution yet sent back
            remote_cells->vul_wr[i] = TLq.vul_rd[i];  // this is the corresponding tgt cell (arrival)
            remote_cells->inc_n();
            // set the global id on new elem
            rval = mb->tag_set_data( gid, &new_element, 1, &globalIdEl );ERRORR( rval, "can't set global id tag on new el" );
        }
    }
    // order the remote cells tuple list, with the global id, because we will search in it
    // remote_cells->print("remote_cells before sorting");
    moab::TupleList::buffer sort_buffer;
    sort_buffer.buffer_init( n );
    remote_cells->sort( 1, &sort_buffer );
    sort_buffer.reset();
    return MB_SUCCESS;
}

// this algorithm assumes lagr set is already created, and some elements will be coming from
// other procs, and populate the covering_set
// we need to keep in a tuple list the remote cells from other procs, because we need to send back
// the intersection info (like area of the intx polygon, and the current concentration) maybe total
// mass in that intx
ErrorCode Intx2Mesh::create_departure_mesh_3rd_alg( EntityHandle& lagr_set, EntityHandle& covering_set )
{
    EntityHandle dum = 0;

    Tag corrTag;
    ErrorCode rval = mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE | MB_TAG_CREAT, &dum );<--- rval is initialized
    // start copy from 2nd alg
    // compute the bounding box on each proc
    assert( parcomm != NULL );
    if( 1 == parcomm->proc_config().proc_size() )
    {
        covering_set = lagr_set;  // nothing to communicate, it must be serial
        return MB_SUCCESS;
    }

    // get all local verts
    Range local_verts;
    rval                = mb->get_connectivity( localEnts, local_verts );<--- rval is overwritten
    int num_local_verts = (int)local_verts.size();ERRORR( rval, "can't get local vertices" );

    std::vector< int > gids( num_local_verts );
    rval = mb->tag_get_data( gid, local_verts, &gids[0] );ERRORR( rval, "can't get local vertices gids" );

    Range localDepCells;
    rval = mb->get_entities_by_dimension( lagr_set, 2, localDepCells );ERRORR( rval, "can't get ents by dimension from lagr set" );

    // get all lagr verts (departure vertices)
    Range lagr_verts;
    rval = mb->get_connectivity( localDepCells, lagr_verts );  // they should be created in
    // the same order as the euler vertices
    int num_lagr_verts = (int)lagr_verts.size();ERRORR( rval, "can't get local lagr vertices" );

    // now see the departure points position; to what boxes should we send them?
    std::vector< double > dep_points( 3 * num_lagr_verts );
    rval = mb->get_coords( lagr_verts, &dep_points[0] );ERRORR( rval, "can't get departure points position" );
    // ranges to send to each processor; will hold vertices and elements (quads?)
    // will look if the box of the dep quad covers box of euler mesh on proc (with tolerances)
    std::map< int, Range > Rto;
    int numprocs = parcomm->proc_config().proc_size();

    for( Range::iterator eit = localDepCells.begin(); eit != localDepCells.end(); ++eit )
    {
        EntityHandle q = *eit;
        const EntityHandle* conn4;
        int num_nodes;
        rval = mb->get_connectivity( q, conn4, num_nodes );ERRORR( rval, "can't get DP tag values" );
        CartVect qbmin( DBL_MAX );
        CartVect qbmax( -DBL_MAX );
        for( int i = 0; i < num_nodes; i++ )
        {
            EntityHandle v = conn4[i];
            int index      = lagr_verts.index( v );
            assert( -1 != index );
            CartVect dp( &dep_points[3 * index] );  // will use constructor
            for( int j = 0; j < 3; j++ )
            {
                if( qbmin[j] > dp[j] ) qbmin[j] = dp[j];
                if( qbmax[j] < dp[j] ) qbmax[j] = dp[j];
            }
        }
        for( int p = 0; p < numprocs; p++ )
        {
            CartVect bbmin( &allBoxes[6 * p] );
            CartVect bbmax( &allBoxes[6 * p + 3] );
            if( GeomUtil::boxes_overlap( bbmin, bbmax, qbmin, qbmax, box_error ) )
            {
                Rto[p].insert( q );
            }
        }
    }

    // now, build TLv and TLq, for each p
    size_t numq = 0;
    size_t numv = 0;
    for( int p = 0; p < numprocs; p++ )
    {
        if( p == (int)my_rank ) continue;  // do not "send" it, because it is already here
        Range& range_to_P = Rto[p];
        // add the vertices to it
        if( range_to_P.empty() ) continue;  // nothing to send to proc p
        Range vertsToP;
        rval = mb->get_connectivity( range_to_P, vertsToP );ERRORR( rval, "can't get connectivity" );
        numq = numq + range_to_P.size();
        numv = numv + vertsToP.size();
        range_to_P.merge( vertsToP );
    }
    TupleList TLv;
    TupleList TLq;
    TLv.initialize( 2, 0, 0, 3, numv );  // to proc, GLOBAL ID, DP points
    TLv.enableWriteAccess();

    int sizeTuple = 2 + max_edges_1;  // max edges could be up to MAXEDGES :) for polygons
    TLq.initialize( 2 + max_edges_1, 0, 1, 0,
                    numq );  // to proc, elem GLOBAL ID, connectivity[max_edges] (global ID v)
    // send also the corresponding tgt cell it will come to
    TLq.enableWriteAccess();
#ifdef VERBOSE
    std::cout << "from proc " << my_rank << " send " << numv << " vertices and " << numq << " elements\n";
#endif

    for( int to_proc = 0; to_proc < numprocs; to_proc++ )
    {
        if( to_proc == (int)my_rank ) continue;
        Range& range_to_P = Rto[to_proc];
        Range V           = range_to_P.subset_by_type( MBVERTEX );

        for( Range::iterator it = V.begin(); it != V.end(); ++it )
        {
            EntityHandle v = *it;
            int index = lagr_verts.index( v );  // will be the same index as the corresponding vertex in euler verts
            assert( -1 != index );
            int n                = TLv.get_n();
            TLv.vi_wr[2 * n]     = to_proc;                // send to processor
            TLv.vi_wr[2 * n + 1] = gids[index];            // global id needs index in the local_verts range
            TLv.vr_wr[3 * n]     = dep_points[3 * index];  // departure position, of the node local_verts[i]
            TLv.vr_wr[3 * n + 1] = dep_points[3 * index + 1];
            TLv.vr_wr[3 * n + 2] = dep_points[3 * index + 2];
            TLv.inc_n();
        }
        // also, prep the 2d cells for sending ...
        Range Q = range_to_P.subset_by_dimension( 2 );
        for( Range::iterator it = Q.begin(); it != Q.end(); ++it )
        {
            EntityHandle q = *it;  // this is a src cell
            int global_id;
            rval = mb->tag_get_data( gid, &q, 1, &global_id );ERRORR( rval, "can't get gid for polygon" );
            int n                        = TLq.get_n();
            TLq.vi_wr[sizeTuple * n]     = to_proc;    //
            TLq.vi_wr[sizeTuple * n + 1] = global_id;  // global id of element, used to identify it ...
            const EntityHandle* conn4;
            int num_nodes;
            rval = mb->get_connectivity(<--- rval is assigned
                q, conn4, num_nodes );  // could be up to 10;ERRORR( rval, "can't get connectivity for quad" );
            if( num_nodes > MAXEDGES ) ERRORR( MB_FAILURE, "too many nodes in a polygon" );
            for( int i = 0; i < num_nodes; i++ )
            {
                EntityHandle v = conn4[i];
                int index      = lagr_verts.index( v );
                assert( -1 != index );
                TLq.vi_wr[sizeTuple * n + 2 + i] = gids[index];
            }
            for( int k = num_nodes; k < max_edges_1; k++ )
            {
                TLq.vi_wr[sizeTuple * n + 2 + k] =
                    0;  // fill the rest of node ids with 0; we know that the node ids start from 1!
            }
            EntityHandle tgtCell;
            rval = mb->tag_get_data( corrTag, &q, 1, &tgtCell );ERRORR( rval, "can't get corresponding tgt cell for dep cell" );<--- rval is overwritten
            TLq.vul_wr[n] = tgtCell;  // this will be sent to remote_cells, to be able to come back
            TLq.inc_n();
        }
    }
    // now we can route them to each processor
    // now we are done populating the tuples; route them to the appropriate processors
    ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, TLv, 0 );
    ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, TLq, 0 );
    // the elements are already in localEnts;

    // maps from global ids to new vertex and quad handles, that are added
    std::map< int, EntityHandle > globalID_to_handle;
    // we already have vertices from lagr set; they are already in the processor, even before
    // receiving other verts from neighbors
    int k = 0;
    for( Range::iterator vit = lagr_verts.begin(); vit != lagr_verts.end(); ++vit, k++ )
    {
        globalID_to_handle[gids[k]] = *vit;  // a little bit of overkill
        // we do know that the global ids between euler and lagr verts are parallel
    }
    /*std::map<int, EntityHandle> globalID_to_eh;*/  // do we need this one?
    globalID_to_eh.clear();
    // now, look at every TLv, and see if we have to create a vertex there or not
    int n = TLv.get_n();  // the size of the points received
    for( int i = 0; i < n; i++ )
    {
        int globalId = TLv.vi_rd[2 * i + 1];
        if( globalID_to_handle.find( globalId ) == globalID_to_handle.end() )
        {
            EntityHandle new_vert;
            double dp_pos[3] = { TLv.vr_wr[3 * i], TLv.vr_wr[3 * i + 1], TLv.vr_wr[3 * i + 2] };
            rval             = mb->create_vertex( dp_pos, new_vert );ERRORR( rval, "can't create new vertex " );
            globalID_to_handle[globalId] = new_vert;
        }
    }

    // now, all dep points should be at their place
    // look in the local list of 2d cells for this proc, and create all those cells if needed
    // it may be an overkill, but because it does not involve communication, we do it anyway
    Range& local  = Rto[my_rank];
    Range local_q = local.subset_by_dimension( 2 );
    // the local should have all the vertices in lagr_verts
    for( Range::iterator it = local_q.begin(); it != local_q.end(); ++it )
    {
        EntityHandle q = *it;  // these are from lagr cells, local
        int gid_el;
        rval = mb->tag_get_data( gid, &q, 1, &gid_el );ERRORR( rval, "can't get element global ID " );
        globalID_to_eh[gid_el] = q;  // do we need this? maybe to just mark the ones on this processor
        // maybe a range of global cell ids is fine?
    }
    // now look at all elements received through; we do not want to duplicate them
    n = TLq.get_n();  // number of elements received by this processor
    // a cell should be received from one proc only; so why are we so worried about duplicated
    // elements? a vertex can be received from multiple sources, that is fine

    remote_cells = new TupleList();
    remote_cells->initialize( 2, 0, 1, 0, n );  // no tracers anymore in these tuples
    remote_cells->enableWriteAccess();
    for( int i = 0; i < n; i++ )
    {
        int globalIdEl = TLq.vi_rd[sizeTuple * i + 1];
        int from_proc  = TLq.vi_rd[sizeTuple * i];
        // do we already have a quad with this global ID, represented?
        if( globalID_to_eh.find( globalIdEl ) == globalID_to_eh.end() )
        {
            // construct the conn quad
            EntityHandle new_conn[MAXEDGES];
            int nnodes = -1;
            for( int j = 0; j < max_edges_1; j++ )
            {
                int vgid = TLq.vi_rd[sizeTuple * i + 2 + j];  // vertex global ID
                if( vgid == 0 )
                    new_conn[j] = 0;
                else
                {
                    assert( globalID_to_handle.find( vgid ) != globalID_to_handle.end() );
                    new_conn[j] = globalID_to_handle[vgid];
                    nnodes      = j + 1;  // nodes are at the beginning, and are variable number
                }
            }
            EntityHandle new_element;
            //
            EntityType entType = MBQUAD;
            if( nnodes > 4 ) entType = MBPOLYGON;
            if( nnodes < 4 ) entType = MBTRI;
            rval = mb->create_element( entType, new_conn, nnodes, new_element );ERRORR( rval, "can't create new element " );
            globalID_to_eh[globalIdEl] = new_element;
            local_q.insert( new_element );
            rval = mb->tag_set_data( gid, &new_element, 1, &globalIdEl );ERRORR( rval, "can't set gid on new element " );
        }
        remote_cells->vi_wr[2 * i]     = from_proc;
        remote_cells->vi_wr[2 * i + 1] = globalIdEl;
        //    remote_cells->vr_wr[i] = 0.; will have a different tuple for communication
        remote_cells->vul_wr[i] = TLq.vul_rd[i];  // this is the corresponding tgt cell (arrival)
        remote_cells->inc_n();
    }
    // now, create a new set, covering_set
    rval = mb->create_meshset( MESHSET_SET, covering_set );ERRORR( rval, "can't create new mesh set " );
    rval = mb->add_entities( covering_set, local_q );ERRORR( rval, "can't add entities to new mesh set " );
    // order the remote cells tuple list, with the global id, because we will search in it
    // remote_cells->print("remote_cells before sorting");
    moab::TupleList::buffer sort_buffer;
    sort_buffer.buffer_init( n );
    remote_cells->sort( 1, &sort_buffer );
    sort_buffer.reset();
    return MB_SUCCESS;
    // end copy
}

ErrorCode Intx2Mesh::resolve_intersection_sharing()
{
    if( parcomm && parcomm->size() > 1 )
    {
        /*
            moab::ParallelMergeMesh pm(parcomm, epsilon_1);
            ErrorCode rval = pm.merge(outSet, false, 2); // resolve only the output set, do not skip
           local merge, use dim 2 ERRORR(rval, "can't merge intersection ");
        */
        // look at non-owned shared vertices, that could be part of original source set
        // they should be removed from intx set reference, because they might not have a
        // correspondent on the other task
        Range nonOwnedVerts;
        Range vertsInIntx;
        Range intxCells;
        ErrorCode rval = mb->get_entities_by_dimension( outSet, 2, intxCells );MB_CHK_ERR( rval );
        rval = mb->get_connectivity( intxCells, vertsInIntx );MB_CHK_ERR( rval );

        rval = parcomm->filter_pstatus( vertsInIntx, PSTATUS_NOT_OWNED, PSTATUS_AND, -1, &nonOwnedVerts );MB_CHK_ERR( rval );

        // some of these vertices can be in original set 1, which was covered, transported;
        // but they should not be "shared" from the intx point of view, because they are not shared
        // with another task they might have come from coverage as a plain vertex, so losing the
        // sharing property ?

        Range coverVerts;
        rval = mb->get_connectivity( rs1, coverVerts );MB_CHK_ERR( rval );
        // find out those that are on the interface
        Range vertsCovInterface;
        rval = parcomm->filter_pstatus( coverVerts, PSTATUS_INTERFACE, PSTATUS_AND, -1, &vertsCovInterface );MB_CHK_ERR( rval );
        // how many of these are in
        Range nodesToDuplicate = intersect( vertsCovInterface, nonOwnedVerts );
        // first, get all cells connected to these vertices, from intxCells

        Range connectedCells;
        rval = mb->get_adjacencies( nodesToDuplicate, 2, false, connectedCells, Interface::UNION );MB_CHK_ERR( rval );
        // only those in intx set:
        connectedCells = intersect( connectedCells, intxCells );
        // first duplicate vertices in question:
        std::map< EntityHandle, EntityHandle > duplicatedVerticesMap;
        for( Range::iterator vit = nodesToDuplicate.begin(); vit != nodesToDuplicate.end(); vit++ )
        {
            EntityHandle vertex = *vit;
            double coords[3];
            rval = mb->get_coords( &vertex, 1, coords );MB_CHK_ERR( rval );
            EntityHandle newVertex;
            rval = mb->create_vertex( coords, newVertex );MB_CHK_ERR( rval );
            duplicatedVerticesMap[vertex] = newVertex;
        }

        // look now at connectedCells, and change their connectivities:
        for( Range::iterator eit = connectedCells.begin(); eit != connectedCells.end(); eit++ )
        {
            EntityHandle intxCell = *eit;
            // replace connectivity
            std::vector< EntityHandle > connectivity;
            rval = mb->get_connectivity( &intxCell, 1, connectivity );MB_CHK_ERR( rval );
            for( size_t i = 0; i < connectivity.size(); i++ )
            {
                EntityHandle currentVertex                           = connectivity[i];
                std::map< EntityHandle, EntityHandle >::iterator mit = duplicatedVerticesMap.find( currentVertex );
                if( mit != duplicatedVerticesMap.end() )
                {
                    connectivity[i] = mit->second;  // replace connectivity directly
                }
            }
            int nnodes = (int)connectivity.size();
            rval       = mb->set_connectivity( intxCell, &connectivity[0], nnodes );MB_CHK_ERR( rval );
        }
    }
    return MB_SUCCESS;
}
#endif /* MOAB_HAVE_MPI */
} /* namespace moab */