HalfFacetRep.cpp 114 KB
Newer Older
1 2 3 4 5
/**
 * MOAB, a Mesh-Oriented datABase, is a software component for creating,
 * storing and accessing finite element mesh data.
 * 
 * Copyright 2004 Sandia Corporation.  Under the terms of Contract
6
 * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
7 8 9 10 11 12 13 14 15 16
 * retains certain rights in this software.
 * 
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * 
 */


17
#include "moab/HalfFacetRep.hpp"
18
#include "Internals.hpp"
19 20 21
#include <iostream>
#include <assert.h>
#include <vector>
22
#include <map>
23
#include "MBTagConventions.hpp"
nray's avatar
nray committed
24
#include "moab/ScdInterface.hpp"
25
#ifdef MOAB_HAVE_MPI
nray's avatar
nray committed
26 27
#include "moab/ParallelComm.hpp"
#endif
28 29 30

namespace moab {

31 32 33 34 35 36 37 38 39 40 41 42 43 44 45
  inline EntityHandle CREATE_HALFFACET(const unsigned lid, const EntityID id)
  {
    assert(id <= MB_END_ID && lid < 6);
    return (((HFacet)lid) << MB_ID_WIDTH)|id;
  }
  inline EntityID FID_FROM_HALFFACET(HFacet handle)
  {
    return (handle & MB_ID_MASK);
  }
  inline int LID_FROM_HALFFACET(HFacet handle)
  {
    return static_cast<int> (handle >> MB_ID_WIDTH);
  }


46 47
  HalfFacetRep::HalfFacetRep(Core *impl,   ParallelComm *comm, moab::EntityHandle rset)
    : mb(impl), pcomm(comm), _rset(rset)
48 49 50
  {
    assert(NULL != impl);
    mInitAHFmaps = false;
51 52
    chk_mixed = false;
    is_mixed = false;
53 54
  }

55 56
  HalfFacetRep::~HalfFacetRep() {}

57

58 59
  MESHTYPE HalfFacetRep::get_mesh_type(int nverts, int nedges, int nfaces, int ncells)
  {
60
    MESHTYPE mesh_type = CURVE;
61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79

    if (nverts && nedges && (!nfaces) && (!ncells))
      mesh_type = CURVE;
    else if (nverts && !nedges && nfaces && !ncells)
      mesh_type = SURFACE;
    else if (nverts && nedges && nfaces && !ncells)
      mesh_type = SURFACE_MIXED;
    else if (nverts && !nedges && !nfaces && ncells)
      mesh_type = VOLUME;
    else if (nverts && nedges && !nfaces && ncells)
      mesh_type = VOLUME_MIXED_1;
    else if (nverts && !nedges && nfaces && ncells)
      mesh_type = VOLUME_MIXED_2;
    else if (nverts && nedges && nfaces && ncells)
      mesh_type = VOLUME_MIXED;

    return mesh_type;
  }

80 81 82 83 84 85 86 87 88 89
  const HalfFacetRep::adj_matrix HalfFacetRep::adjMatrix[7] =
  {
      // Stores the adjacency matrix for each mesh type.
      //CURVE
      {{{0,1,0,0},{1,1,0,0},{0,0,0,0},{0,0,0,0}}},

      //SURFACE
      {{{0,0,1,0},{0,0,0,0},{1,0,1,0},{0,0,0,0}}},

      //SURFACE_MIXED
90
      {{{0,1,1,0},{1,1,1,0},{1,1,1,0},{0,0,0,0}}},
91 92 93 94 95

      //VOLUME
      {{{0,0,0,1},{0,0,0,0},{0,0,0,0},{1,0,0,1}}},

      //VOLUME_MIXED_1
96
      {{{0,1,0,1},{1,1,0,1},{0,0,0,0},{1,1,0,1}}},
97 98

      //VOLUME_MIXED_2
99
      {{{0,0,1,1},{0,0,0,0},{1,0,1,1},{1,0,1,1}}},
100 101

      //VOLUME_MIXED
102
      {{{0,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1}}}
103 104 105 106
  };

  int HalfFacetRep::get_index_for_meshtype(MESHTYPE mesh_type)
  {
107
      int index = 0;
nray's avatar
nray committed
108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130
      switch (mesh_type) {
        case CURVE:
          index = 0;
          break;
        case SURFACE:
          index = 1;
          break;
        case SURFACE_MIXED:
          index = 2;
          break;
        case VOLUME:
          index = 3;
          break;
        case VOLUME_MIXED_1:
          index = 4;
          break;
        case VOLUME_MIXED_2:
          index = 5;
          break;
        case VOLUME_MIXED:
          index = 6;
          break;
        }
131 132 133 134 135 136 137 138 139 140 141 142
      return index;
  }

  bool HalfFacetRep::check_mixed_entity_type()
  {
      if (!chk_mixed)
      {
          chk_mixed = true;

          ErrorCode error;
          Range felems, celems;

143
          error = mb->get_entities_by_dimension( this->_rset, 2, felems ); MB_CHK_ERR(error);
144 145 146 147 148 149 150 151 152 153 154

          if (felems.size()){
              Range tris, quad, poly;
              tris = felems.subset_by_type(MBTRI);
              quad = felems.subset_by_type(MBQUAD);
              poly = felems.subset_by_type(MBPOLYGON);
              if ((tris.size()&&quad.size())||(tris.size()&&poly.size())||(quad.size()&&poly.size()))
                  is_mixed = true;
              if (poly.size())
                  is_mixed = true;

155
              if (is_mixed) return is_mixed;
156 157
          }

nray's avatar
nray committed
158
          error = mb->get_entities_by_dimension( this->_rset, 3, celems);   MB_CHK_ERR(error);
159 160 161 162 163 164 165
          if (celems.size()){
              Range tet, pyr, prism, hex, polyhed;
              tet = celems.subset_by_type(MBTET);
              pyr = celems.subset_by_type(MBPYRAMID);
              prism = celems.subset_by_type(MBPRISM);
              hex = celems.subset_by_type(MBHEX);
              polyhed = celems.subset_by_type(MBPOLYHEDRON);
166 167 168 169 170
              if ((tet.size() && pyr.size())||(tet.size() && prism.size())||
                  (tet.size() && hex.size())||(tet.size()&&polyhed.size())||
                  (pyr.size() && prism.size())||(pyr.size() && hex.size()) ||
                  (pyr.size()&&polyhed.size())|| (prism.size() && hex.size())||
                  (prism.size()&&polyhed.size())||(hex.size()&&polyhed.size()))
171 172 173 174
                  is_mixed = true;

              if (polyhed.size())
                  is_mixed = true;
nray's avatar
nray committed
175 176 177 178 179 180 181 182 183 184
          }

          ScdInterface *scdi = NULL;
          error = mb->query_interface(scdi);MB_CHK_ERR(error);
          if (scdi){
              Range boxes;
              error = scdi->find_boxes(boxes);MB_CHK_ERR(error);

              if (!boxes.empty())
                is_mixed = true;
185 186 187 188 189
          }
      }
      return is_mixed;
  }

190 191 192 193
  /*******************************************************
   * initialize                                          *
   ******************************************************/

194 195
  ErrorCode HalfFacetRep::initialize()
  {
196 197
    ErrorCode error;

198 199
    if (!mInitAHFmaps){
        mInitAHFmaps = true;
200
#ifdef MOAB_HAVE_MPI
nray's avatar
nray committed
201
        if (pcomm){
202 203 204 205 206 207 208 209 210 211 212 213
            moab::Range _averts, _aedgs, _afacs, _acels;
            error = mb->get_entities_by_dimension(this->_rset, 0, _averts, true);MB_CHK_ERR(error);
            error = mb->get_entities_by_dimension(this->_rset, 1, _aedgs, true);MB_CHK_ERR(error);
            error = mb->get_entities_by_dimension(this->_rset, 2, _afacs, true);MB_CHK_ERR(error);
            error = mb->get_entities_by_dimension(this->_rset, 3, _acels, true);MB_CHK_ERR(error);

            // filter based on parallel status
            error = pcomm->filter_pstatus(_averts, PSTATUS_GHOST, PSTATUS_NOT, -1, &_verts);MB_CHK_ERR(error);
            error = pcomm->filter_pstatus(_aedgs, PSTATUS_GHOST, PSTATUS_NOT, -1, &_edges);MB_CHK_ERR(error);
            error = pcomm->filter_pstatus(_afacs, PSTATUS_GHOST, PSTATUS_NOT, -1, &_faces);MB_CHK_ERR(error);
            error = pcomm->filter_pstatus(_acels, PSTATUS_GHOST, PSTATUS_NOT, -1, &_cells);MB_CHK_ERR(error);
          }
nray's avatar
nray committed
214
        else {
215 216 217 218 219
            error = mb->get_entities_by_dimension( this->_rset, 0, _verts, true);MB_CHK_ERR(error);
            error = mb->get_entities_by_dimension( this->_rset, 1, _edges, true);MB_CHK_ERR(error);
            error = mb->get_entities_by_dimension( this->_rset, 2, _faces, true);MB_CHK_ERR(error);
            error = mb->get_entities_by_dimension( this->_rset, 3, _cells, true);MB_CHK_ERR(error);
          }
nray's avatar
nray committed
220 221 222 223 224 225 226
#else
        error = mb->get_entities_by_dimension( this->_rset, 0, _verts, true);MB_CHK_ERR(error);
        error = mb->get_entities_by_dimension( this->_rset, 1, _edges, true);MB_CHK_ERR(error);
        error = mb->get_entities_by_dimension( this->_rset, 2, _faces, true);MB_CHK_ERR(error);
        error = mb->get_entities_by_dimension( this->_rset, 3, _cells, true);MB_CHK_ERR(error);

#endif
227 228 229 230 231 232 233 234 235

        int nverts = _verts.size();
        int nedges = _edges.size();
        int nfaces = _faces.size();
        int ncells = _cells.size();


        MESHTYPE mesh_type = get_mesh_type(nverts, nedges, nfaces, ncells);
        thismeshtype = mesh_type;
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
        //Initialize mesh type specific maps
        if (thismeshtype == CURVE){
            error = init_curve(); MB_CHK_ERR(error);
          }
        else if (thismeshtype == SURFACE){
            error = init_surface();MB_CHK_ERR(error);
          }
        else if (thismeshtype == SURFACE_MIXED){
            error = init_curve();MB_CHK_ERR(error);
            error = init_surface();MB_CHK_ERR(error);
          }
        else if (thismeshtype == VOLUME){
            error = init_volume();MB_CHK_ERR(error);
          }
        else if (thismeshtype == VOLUME_MIXED_1){
            error = init_curve();MB_CHK_ERR(error);
            error = init_volume();MB_CHK_ERR(error);
          }
        else if (thismeshtype == VOLUME_MIXED_2){
            error = init_surface();MB_CHK_ERR(error);
            error = init_volume();MB_CHK_ERR(error);
          }
        else if (thismeshtype == VOLUME_MIXED){
            error = init_curve();MB_CHK_ERR(error);
            error = init_surface();MB_CHK_ERR(error);
            error = init_volume();MB_CHK_ERR(error);
          }
      }
265 266 267
    return MB_SUCCESS;
  }

268 269 270 271 272
  ErrorCode HalfFacetRep::deinitialize()
  {
    return MB_SUCCESS;
  }

273 274 275 276
  ErrorCode HalfFacetRep::init_curve()
  {
    ErrorCode error;

277 278
    int nv = ID_FROM_HANDLE(*(_verts.end()-1));
    int ne = ID_FROM_HANDLE(*(_edges.end()-1));
279

280 281
    v2hv.resize(nv,0);
    sibhvs.resize(ne*2,0);
282

nray's avatar
nray committed
283
    error = determine_sibling_halfverts(_verts, _edges);MB_CHK_ERR(error);
284
    error = determine_incident_halfverts( _edges);MB_CHK_ERR(error);
285 286 287 288 289 290 291

    return MB_SUCCESS;
  }

  ErrorCode HalfFacetRep::init_surface()
  {
    ErrorCode error;
292 293
    EntityType ftype = mb->type_from_handle(*_faces.begin());
    int nepf = lConnMap2D[ftype-2].num_verts_in_face;
nray's avatar
nray committed
294

295 296
    int nv = ID_FROM_HANDLE(*(_verts.end()-1));
    int nf = ID_FROM_HANDLE(*(_faces.end()-1));
nray's avatar
nray committed
297

298 299
    v2he.resize(nv, 0);
    sibhes.resize(nf*nepf,0);
300 301

    // Construct ahf maps
nray's avatar
nray committed
302 303
    error = determine_sibling_halfedges(_faces);MB_CHK_ERR(error);
    error = determine_incident_halfedges(_faces);MB_CHK_ERR(error);
304

305 306 307 308 309
    //Initialize queues for storing face and local id's during local search
    for (int i = 0; i< MAXSIZE; i++)
      {
        queue_fid[i] = 0;
        queue_lid[i] = 0;
nray's avatar
nray committed
310
        trackfaces[i] = 0;
311 312
      }

313 314 315 316 317 318 319
    return MB_SUCCESS;
  }

  ErrorCode HalfFacetRep::init_volume()
  {
    ErrorCode error;

320 321 322 323 324 325
    //Initialize std::map between cell-types and their index in lConnMap3D
    cell_index[MBTET] = 0;
    cell_index[MBPYRAMID] = 1;
    cell_index[MBPRISM] = 2;
    cell_index[MBHEX] = 3;

nray's avatar
nray committed
326
    int index = get_index_in_lmap(*_cells.begin());
327
    int nfpc = lConnMap3D[index].num_faces_in_cell;
328 329
    int nv = ID_FROM_HANDLE(*(_verts.end()-1));
    int nc = ID_FROM_HANDLE(*(_cells.end()-1));;
330

331 332
    v2hf.resize(nv, 0);
    sibhfs.resize(nc*nfpc, 0);
333

334
    //Construct the maps
nray's avatar
nray committed
335 336
    error = determine_sibling_halffaces(_cells);MB_CHK_ERR(error);
    error = determine_incident_halffaces(_cells);MB_CHK_ERR(error);
337

nray's avatar
nray committed
338 339 340 341 342 343 344
    //Initialize queues for storing face and local id's during local search
    for (int i = 0; i< MAXSIZE; i++)
      {
        Stkcells[i] = 0;
        cellq[i] = 0;
        trackcells[i] = 0;
      }
345

346 347 348 349
    return MB_SUCCESS;
  }

   //////////////////////////////////////////////////
350
   ErrorCode HalfFacetRep::print_tags(int dim)
351
   {
352 353 354 355 356 357 358
     if (dim==1){
         EntityHandle start_edge = *_edges.begin();
         std::cout<<"start_edge = "<<start_edge<<std::endl;
         std::cout<<"<SIBHVS_EID,SIBHVS_LVID>"<<std::endl;

         for (Range::iterator i = _edges.begin(); i != _edges.end(); ++i){
             EntityHandle eid[2];  int lvid[2];
359
             int eidx = ID_FROM_HANDLE(*i)-1;
360 361 362 363 364 365
             HFacet hf1 = sibhvs[2*eidx];
             HFacet hf2 = sibhvs[2*eidx+1];
             eid[0] = fid_from_halfacet(hf1,MBEDGE); eid[1] = fid_from_halfacet(hf2,MBEDGE);
             lvid[0] = lid_from_halffacet(hf1); lvid[1] = lid_from_halffacet(hf2);
             std::cout<<"Entity = "<<*i<<" :: <"<<eid[0]<<","<<lvid[0]<<">"<<"      "<<"<"<<eid[1]<<","<<lvid[1]<<">"<<std::endl;
           }
366

367
         std::cout<<"<V2HV_EID, V2HV_LVID>"<<std::endl;
368

369
         for (Range::iterator i = _verts.begin(); i != _verts.end(); ++i){
370
             int vidx = ID_FROM_HANDLE(*i)-1;
371 372 373 374 375
             HFacet hf = v2hv[vidx];
             EntityHandle eid = fid_from_halfacet(hf, MBEDGE);
             int lvid = lid_from_halffacet(hf);
             std::cout<<"Vertex = "<<*i<<" :: <"<<eid<<","<<lvid<<">"<<std::endl;
           }
376
       }
377 378 379 380 381 382 383 384
     else if (dim==2){
         EntityType ftype = mb->type_from_handle(*_faces.begin());
         int nepf = lConnMap2D[ftype-2].num_verts_in_face;
         EntityHandle start_face = *_faces.begin();
         std::cout<<"start_face = "<<start_face<<std::endl;
         std::cout<<"<SIBHES_FID,SIBHES_LEID>"<<std::endl;

         for (Range::iterator i = _faces.begin(); i != _faces.end(); ++i){
385
             int fidx = ID_FROM_HANDLE(*i)-1;
386 387 388 389 390 391 392 393 394
             std::cout<<"Entity = "<<*i;
             for (int j=0; j<nepf; j++){
                 HFacet hf = sibhes[nepf*fidx+j];
                 EntityHandle sib = fid_from_halfacet(hf, ftype);
                 int lid = lid_from_halffacet(hf);
                 std::cout<<" :: <"<<sib<<","<<lid<<">"<<"       ";
               }
             std::cout<<std::endl;
           }
395

396
         std::cout<<"<V2HE_FID, V2HE_LEID>"<<std::endl;
397

398
         for (Range::iterator i = _verts.begin(); i != _verts.end(); ++i){
399
             int vidx = ID_FROM_HANDLE(*i)-1;
400 401
             HFacet hf = v2he[vidx];
             EntityHandle fid = fid_from_halfacet(hf, ftype);
402
             int lid = lid_from_halffacet(hf);
403
             std::cout<<"Vertex = "<<*i<<" :: <"<<fid<<","<<lid<<">"<<std::endl;
404
           }
405
       }
406 407 408 409 410 411 412 413 414 415
     else if (dim==3){
         EntityType ctype = mb->type_from_handle(*_cells.begin());

         int index = get_index_in_lmap(*_cells.begin());
         int nfpc = lConnMap3D[index].num_faces_in_cell;
         EntityHandle start_cell = *_cells.begin();
         std::cout<<"start_cell = "<<start_cell<<std::endl;
         std::cout<<"<SIBHES_CID,SIBHES_LFID>"<<std::endl;

         for (Range::iterator i = _cells.begin(); i != _cells.end(); ++i){
416
             int cidx = ID_FROM_HANDLE(*i)-1;
417 418 419 420 421 422 423 424 425
             std::cout<<"Entity = "<<*i;
             for (int j=0; j<nfpc; j++){
                 HFacet hf = sibhfs[nfpc*cidx+j];
                 EntityHandle sib = fid_from_halfacet(hf, ctype);
                 int lid = lid_from_halffacet(hf);
                 std::cout<<" :: <"<<sib<<","<<lid<<">"<<"       ";
               }
             std::cout<<std::endl;
           }
426

427
         std::cout<<"<V2HF_CID, V2HF_LFID>"<<std::endl;
428 429
         EntityHandle cid;
         int lid;
430

431
         for (Range::iterator i = _verts.begin(); i != _verts.end(); ++i){
432
             int vidx = ID_FROM_HANDLE(*i)-1;
433
             HFacet hf = v2hf[vidx];
434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452

             if (hf ==0 && (v2hfs.find(*i) != v2hfs.end()))
               {
                 std::pair <std::multimap<EntityHandle, HFacet>::iterator, std::multimap<EntityHandle, HFacet>::iterator> it_hfs;
                 it_hfs = v2hfs.equal_range(*i);

                 for (std::multimap<EntityHandle, HFacet>::iterator it = it_hfs.first; it != it_hfs.second; ++it)
                   {
                     cid = fid_from_halfacet(it->second, ctype);
                     lid = lid_from_halffacet(hf);

                     std::cout<<"Vertex = "<<*i<<" :: <"<<cid<<","<<lid<<">"<<std::endl;
                   }
               }
             else {
                 cid = fid_from_halfacet(hf, ctype);
                 lid = lid_from_halffacet(hf);
                 std::cout<<"Vertex = "<<*i<<" :: <"<<cid<<","<<lid<<">"<<std::endl;
               }
453
           }
454

455
       }
456 457 458 459 460 461 462
     return MB_SUCCESS;
   }

   /**********************************************************
   *      User interface for adjacency functions            *
   ********************************************************/

463 464 465 466 467 468 469 470
  ErrorCode HalfFacetRep::get_adjacencies(const EntityHandle source_entity,
                                          const unsigned int target_dimension,
                                          std::vector<EntityHandle> &target_entities)
  {

      ErrorCode error;

      unsigned int source_dimension = mb->dimension_from_handle(source_entity);
471
      assert((source_dimension <= target_dimension) || (source_dimension > target_dimension));
472 473 474

      if (mInitAHFmaps == false)
      {
nray's avatar
nray committed
475
          error = initialize(); MB_CHK_ERR(error);
476 477
      }

478 479 480 481
      int mindex = get_index_for_meshtype(thismeshtype);
      int adj_possible = adjMatrix[mindex].val[source_dimension][target_dimension];

      if (adj_possible)
482
      {
483
          if (source_dimension < target_dimension)
484
          {
nray's avatar
nray committed
485
              error = get_up_adjacencies(source_entity, target_dimension, target_entities); MB_CHK_ERR(error);
486
          }
487
          else if (source_dimension == target_dimension)
488
          {
nray's avatar
nray committed
489
              error = get_neighbor_adjacencies(source_entity, target_entities);  MB_CHK_ERR(error);
490
          }
491
          else
492
          {           
nray's avatar
nray committed
493
              error = get_down_adjacencies(source_entity, target_dimension, target_entities); MB_CHK_ERR(error);
494
          }
495
      }
496 497 498
      else
          return MB_SUCCESS;

499 500 501 502 503 504 505 506
      return MB_SUCCESS;
  }


   ErrorCode HalfFacetRep::get_up_adjacencies(EntityHandle ent,
                                              int out_dim,
                                              std::vector<EntityHandle> &adjents,
                                              std::vector<int> * lids)
507 508 509
   {
    ErrorCode error;
    int in_dim = mb->dimension_from_handle(ent);
510
    assert((in_dim >=0 && in_dim <= 2) && (out_dim > in_dim));
511

512
    if (in_dim == 0)
513
      {
514 515
        if (out_dim == 1)
        {
nray's avatar
nray committed
516
            error = get_up_adjacencies_1d(ent, adjents, lids);  MB_CHK_ERR(error);
517 518 519
        }
        else if (out_dim == 2)
        {
nray's avatar
nray committed
520
            error = get_up_adjacencies_vert_2d(ent, adjents);  MB_CHK_ERR(error);
521 522 523
        }
        else if (out_dim == 3)
        {
nray's avatar
nray committed
524
            error = get_up_adjacencies_vert_3d(ent, adjents); MB_CHK_ERR(error);
525
        }
526 527 528 529
      }

    else if ((in_dim == 1) && (out_dim == 2))
      {
nray's avatar
nray committed
530
        error = get_up_adjacencies_2d(ent, adjents, lids); MB_CHK_ERR(error);
531 532 533
      }
    else if ((in_dim == 1) && (out_dim == 3))
      {
nray's avatar
nray committed
534
        error = get_up_adjacencies_edg_3d(ent, adjents, lids); MB_CHK_ERR(error);
535 536 537
      }
    else if ((in_dim == 2) && (out_dim ==3))
      {
nray's avatar
nray committed
538
        error = get_up_adjacencies_face_3d(ent, adjents, lids); MB_CHK_ERR(error);
539 540 541 542
      }
    return MB_SUCCESS;
   }

543 544
   ErrorCode HalfFacetRep::get_neighbor_adjacencies(EntityHandle ent,
                                                    std::vector<EntityHandle> &adjents)
545 546 547
   {
     ErrorCode error;
     int in_dim = mb->dimension_from_handle(ent);
548
     assert(in_dim >=1 && in_dim <= 3);
549 550 551

     if (in_dim == 1)
       {
nray's avatar
nray committed
552
         error = get_neighbor_adjacencies_1d(ent, adjents); MB_CHK_ERR(error);
553 554 555 556
       }

     else if (in_dim == 2)
       {
nray's avatar
nray committed
557
         error = get_neighbor_adjacencies_2d(ent, adjents); MB_CHK_ERR(error);
558 559 560
       }
     else if (in_dim == 3)
       {
nray's avatar
nray committed
561
         error = get_neighbor_adjacencies_3d(ent, adjents); MB_CHK_ERR(error);
562 563 564 565
       }
     return MB_SUCCESS;
   }

566 567 568 569
   ErrorCode HalfFacetRep::get_down_adjacencies(EntityHandle ent, int out_dim, std::vector<EntityHandle> &adjents)
   {
       ErrorCode error;
       int in_dim = mb->dimension_from_handle(ent);
570 571
       assert((in_dim >=2 && in_dim <= 3) && (out_dim < in_dim));

572 573
       if ((in_dim == 2)&&(out_dim == 1))
       {
nray's avatar
nray committed
574
           error = get_down_adjacencies_2d(ent, adjents);  MB_CHK_ERR(error);
575 576 577
       }
       else if ((in_dim == 3)&&(out_dim == 1))
       {
nray's avatar
nray committed
578
           error = get_down_adjacencies_edg_3d(ent, adjents); MB_CHK_ERR(error);
579 580 581
       }
       else if ((in_dim == 3)&&(out_dim == 2))
       {
nray's avatar
nray committed
582
           error = get_down_adjacencies_face_3d(ent, adjents); MB_CHK_ERR(error);
583 584 585 586
       }
       return MB_SUCCESS;
   }

587
   ErrorCode HalfFacetRep::count_subentities(Range &edges, Range &faces, Range &cells, int *nedges, int *nfaces)
588 589
   {
     ErrorCode error;
nray's avatar
nray committed
590
     if (edges.size() && !faces.size() && !cells.size())
591
       {
nray's avatar
nray committed
592
         nedges[0] = edges.size();
593 594 595 596 597 598 599 600 601
         nfaces[0] = 0;
       }
     else if (faces.size() && !cells.size())
       {
         nedges[0] = find_total_edges_2d(faces);
         nfaces[0] = 0;
       }
     else if (cells.size())
       {
nray's avatar
nray committed
602
         error = find_total_edges_faces_3d(cells, nedges, nfaces); MB_CHK_ERR(error);
603
       }
604
     return MB_SUCCESS;
605 606
   }

607 608 609
  /******************************************************** 
  * 1D: sibhvs, v2hv, incident and neighborhood queries   *
  *********************************************************/
nray's avatar
nray committed
610
  ErrorCode HalfFacetRep::determine_sibling_halfverts(Range &verts, Range &edges)
611 612 613 614
  {
    ErrorCode error;

    //Step 1: Create an index list storing the starting position for each vertex
615
    int nv = verts.size();
616
    std::vector<int> is_index(nv+1);
617 618 619 620 621
    for (int i =0; i<nv+1; i++)
      is_index[i] = 0;

    for (Range::iterator eid = edges.begin(); eid != edges.end(); ++eid)
      {
622 623
        const EntityHandle* conn;
        int num_conn = 0;
624
        error = mb->get_connectivity(*eid, conn, num_conn, true);MB_CHK_ERR(error);
625

626
        int index = verts.index(conn[0]);
627
        is_index[index+1] += 1;
628
        index = verts.index(conn[1]);
629
        is_index[index+1] += 1;
630 631 632 633 634 635 636
      }
    is_index[0] = 0;

    for (int i=0; i<nv; i++)
      is_index[i+1] = is_index[i] + is_index[i+1];

    //Step 2: Define two arrays v2hv_eid, v2hv_lvid storing every half-facet on a vertex
637 638
    std::vector<EntityHandle> v2hv_map_eid(2*edges.size());
    std::vector<int> v2hv_map_lvid(2*edges.size());
639 640 641

    for (Range::iterator eid = edges.begin(); eid != edges.end(); ++eid)
      {
642 643
        const EntityHandle* conn;
        int num_conn = 0;
644
        error = mb->get_connectivity(*eid, conn, num_conn, true);MB_CHK_ERR(error);
645 646 647

        for (int j = 0; j< 2; j++)
          {
648
            int v = verts.index(conn[j]);
649 650 651 652 653 654 655 656 657 658 659
            v2hv_map_eid[is_index[v]] = *eid;
            v2hv_map_lvid[is_index[v]] = j;
            is_index[v] += 1;
          }
      }

    for (int i=nv-2; i>=0; i--)
      is_index[i+1] = is_index[i];
    is_index[0] = 0;

    //Step 3: Fill up sibling half-verts map
660
    for (Range::iterator vid = verts.begin(); vid != verts.end(); ++vid)
661
      {
662
        int v = verts.index(*vid);
nray's avatar
nray committed
663 664
        int last = is_index[v+1] - 1;
        if (last > is_index[v])
665
          {
nray's avatar
nray committed
666 667
            EntityHandle prev_eid = v2hv_map_eid[last];
            int prev_lvid = v2hv_map_lvid[last];
668

nray's avatar
nray committed
669
            for (int i=is_index[v]; i<=last; i++)
670
              {
nray's avatar
nray committed
671 672
                EntityHandle cur_eid = v2hv_map_eid[i];
                int cur_lvid = v2hv_map_lvid[i];
673

674
                int pidx = ID_FROM_HANDLE(prev_eid)-1;
675
                sibhvs[2*pidx+prev_lvid] = create_halffacet(cur_eid, cur_lvid);
676

nray's avatar
nray committed
677 678
                prev_eid = cur_eid;
                prev_lvid = cur_lvid;
679

680 681 682 683 684 685 686 687
              }
          }
      }

    return MB_SUCCESS;
  }

  ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
688
  ErrorCode HalfFacetRep::determine_incident_halfverts( Range &edges){
689 690
    ErrorCode error;

691 692
    for (Range::iterator e_it = edges.begin(); e_it != edges.end(); ++e_it){
        EntityHandle cur_eid = *e_it;
693 694
        const EntityHandle* conn;
        int num_conn = 0;
695
        error = mb->get_connectivity(*e_it, conn, num_conn, true);MB_CHK_ERR(error);
696 697

        for(int i=0; i<2; ++i){
698
            EntityHandle v = conn[i];
699
            int vidx = ID_FROM_HANDLE(v)-1;
nray's avatar
nray committed
700

701 702
            HFacet hf = v2hv[vidx];
            EntityHandle eid = fid_from_halfacet(hf, MBEDGE);
703
            if (eid==0){
704
                v2hv[vidx] = create_halffacet(cur_eid,i);
705 706 707
              }
          }
      }
708 709 710 711

    return MB_SUCCESS;
  }
  ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
712 713
  ErrorCode  HalfFacetRep::get_up_adjacencies_1d( EntityHandle vid,
                                                  std::vector< EntityHandle > &adjents,
714 715
                                                  std::vector<int> * lvids)
  {
716
    adjents.clear();
717
    adjents.reserve(20);
718 719

    if (lvids != NULL)
nray's avatar
nray committed
720
      lvids->reserve(20);
nray's avatar
nray committed
721

722
    int vidx = ID_FROM_HANDLE(vid)-1;
723
    HFacet hf = v2hv[vidx];
724

nray's avatar
nray committed
725 726 727 728 729
    EntityHandle  start_eid = fid_from_halfacet(hf, MBEDGE);
    int start_lid = lid_from_halffacet(hf);

    EntityHandle eid = start_eid;
    int lid = start_lid;
730 731 732

    if (eid != 0){
      adjents.push_back(eid);
nray's avatar
nray committed
733
      if (lvids != NULL)
734 735
        lvids->push_back(lid);

736
      while (eid !=0) {
737
	  int eidx = ID_FROM_HANDLE(eid)-1;
738 739 740
	  HFacet shf = sibhvs[2*eidx+lid];
	  eid = fid_from_halfacet(shf, MBEDGE);
	  lid = lid_from_halffacet(shf);
nray's avatar
nray committed
741

742 743
	  if ((!eid)||(eid == start_eid))
	    break;
nray's avatar
nray committed
744

745
	  adjents.push_back(eid);
nray's avatar
nray committed
746
	  if (lvids != NULL)
747 748 749 750 751 752 753
	    lvids->push_back(lid);
      }
    }

    return MB_SUCCESS;
  }
  //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
754
  ErrorCode  HalfFacetRep::get_neighbor_adjacencies_1d( EntityHandle eid,
nray's avatar
nray committed
755 756
                                                        std::vector<EntityHandle> &adjents)
  {
757
    adjents.clear();
758 759
    adjents.reserve(20);

nray's avatar
nray committed
760 761
    EntityHandle sibhv_eid;
    int sibhv_lid;
762
    int eidx = ID_FROM_HANDLE(eid)-1;
763 764

    for (int lid = 0;lid < 2; ++lid){
765 766 767
        HFacet shf = sibhvs[2*eidx+lid];
        sibhv_eid = fid_from_halfacet(shf, MBEDGE);
        sibhv_lid = lid_from_halffacet(shf);
768 769 770

      if (sibhv_eid != 0){
        adjents.push_back(sibhv_eid);
nray's avatar
nray committed
771

772
        eidx = ID_FROM_HANDLE(sibhv_eid)-1;
773 774 775
        HFacet nhf = sibhvs[2*eidx+sibhv_lid];
        EntityHandle hv_eid = fid_from_halfacet(nhf,MBEDGE);
        int hv_lid = lid_from_halffacet(nhf);
776 777 778 779 780
	
	while (hv_eid != 0){	    
	  if (hv_eid != eid)
	    adjents.push_back(hv_eid);

781
	  eidx = ID_FROM_HANDLE(hv_eid)-1;
782 783 784
	  HFacet hf = sibhvs[2*eidx+hv_lid];
	  EntityHandle edge = fid_from_halfacet(hf, MBEDGE);
	  if (edge == sibhv_eid)
785
	    break;
nray's avatar
nray committed
786

787 788
	  hv_eid = edge;
	  hv_lid = lid_from_halffacet(hf);
789 790 791 792 793 794 795 796 797 798
	}
      }
    } 

    return MB_SUCCESS;   
  }
  
  /*******************************************************
  * 2D: sibhes, v2he, incident and neighborhood queries  *
  ********************************************************/
799 800 801 802 803 804
  const HalfFacetRep::LocalMaps2D HalfFacetRep::lConnMap2D[2] = {
    //Triangle
    {3, {1,2,0}, {2,0,1}},
    //Quad
    {4,{1,2,3,0},{3,0,1,2}}
  };
805 806 807 808 809
  ///////////////////////////////////////////////////////////////////////////////////////////////////////////////
   ErrorCode HalfFacetRep::determine_sibling_halfedges( Range &faces)
  {
    ErrorCode error;
    EntityHandle start_face = *faces.begin();
810 811 812
    EntityType ftype = mb->type_from_handle(start_face);
    int nfaces = faces.size();
    int nepf = lConnMap2D[ftype-2].num_verts_in_face;
813 814 815

    //Step 1: Create an index list storing the starting position for each vertex
    int nv = _verts.size();
816
    std::vector<int> is_index(nv+1);
817 818 819
    for (int i =0; i<nv+1; i++)
      is_index[i] = 0;

820
    int index;
nray's avatar
nray committed
821

822 823
    for (Range::iterator fid = faces.begin(); fid != faces.end(); ++fid)
       {
824
        const EntityHandle* conn;
825
        error = mb->get_connectivity(*fid, conn, nepf, true);MB_CHK_ERR(error);
826 827 828

         for (int i = 0; i<nepf; i++)
           {
829 830
             index = _verts.index(conn[i]);
             is_index[index+1] += 1;
831 832 833 834 835 836 837 838
           }
       }
     is_index[0] = 0;

     for (int i=0; i<nv; i++)
       is_index[i+1] = is_index[i] + is_index[i+1];

     //Step 2: Define two arrays v2hv_eid, v2hv_lvid storing every half-facet on a vertex
839 840 841
     std::vector<EntityHandle> v2nv(nepf*nfaces);
     std::vector<EntityHandle> v2he_map_fid(nepf*nfaces);
     std::vector<int> v2he_map_leid(nepf*nfaces);
842 843 844

     for (Range::iterator fid = faces.begin(); fid != faces.end(); ++fid)
       {
845
         const EntityHandle* conn;
846
         error = mb->get_connectivity(*fid, conn, nepf, true);MB_CHK_ERR(error);
847 848 849

         for (int j = 0; j< nepf; j++)
           {
850
             int v = _verts.index(conn[j]);
851 852
             int nidx = lConnMap2D[ftype-2].next[j];
             v2nv[is_index[v]] = conn[nidx];
853 854 855 856 857 858 859 860 861 862
             v2he_map_fid[is_index[v]] = *fid;
             v2he_map_leid[is_index[v]] = j;
             is_index[v] += 1;
           }
       }

     for (int i=nv-2; i>=0; i--)
       is_index[i+1] = is_index[i];
     is_index[0] = 0;

863

864 865 866
     //Step 3: Fill up sibling half-verts map
     for (Range::iterator fid = faces.begin(); fid != faces.end(); ++fid)
       {
867
         const EntityHandle* conn;
868
         error = mb->get_connectivity(*fid, conn, nepf, true);MB_CHK_ERR(error);
869

870
         int fidx = ID_FROM_HANDLE(*fid)-1;
871 872
         for (int k =0; k<nepf; k++)
           {
873 874
             HFacet hf = sibhes[nepf*fidx+k];
             EntityHandle sibfid = fid_from_halfacet(hf, ftype);
875 876

             if (sibfid != 0)
877 878
               continue;

879
             int nidx = lConnMap2D[ftype-2].next[k];
880
             int v = _verts.index(conn[k]);
881
             int vn = _verts.index(conn[nidx]);
882 883 884 885 886 887 888

             EntityHandle first_fid = *fid;
             int first_leid = k;

             EntityHandle prev_fid = *fid;
             int prev_leid = k;

889
             for (index = is_index[vn]; index <= is_index[vn+1]-1; index++)
890 891 892 893 894
               {
                 if (v2nv[index] == conn[k])
                   {
                     EntityHandle cur_fid = v2he_map_fid[index];
                     int cur_leid = v2he_map_leid[index];
895

896
                     int pidx = ID_FROM_HANDLE(prev_fid)-1;
897
                     sibhes[nepf*pidx+prev_leid] = create_halffacet(cur_fid, cur_leid);
nray's avatar
nray committed
898

899 900 901 902 903
                     prev_fid = cur_fid;
                     prev_leid = cur_leid;
                   }
               }

904
             for (index = is_index[v]; index <= is_index[v+1]-1; index++)
905
               {
906
                 if ((v2nv[index] == conn[nidx])&&(v2he_map_fid[index] != *fid))
907
                   {
908

909 910
                     EntityHandle cur_fid = v2he_map_fid[index];
                     int cur_leid = v2he_map_leid[index];
911

912
                     int pidx = ID_FROM_HANDLE(prev_fid)-1;
913
                     sibhes[nepf*pidx+prev_leid] = create_halffacet(cur_fid, cur_leid);
nray's avatar
nray committed
914

915 916
                     prev_fid = cur_fid;
                     prev_leid = cur_leid;
917

918 919 920 921
                   }
               }

             if (prev_fid != first_fid){
nray's avatar
nray committed
922

923
                 int pidx = ID_FROM_HANDLE(prev_fid)-1;
924
                 sibhes[nepf*pidx+prev_leid] = create_halffacet(first_fid, first_leid);
925
             }
926 927 928 929 930 931 932 933 934
           }
       }

     return MB_SUCCESS;
  }
  ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  ErrorCode HalfFacetRep::determine_incident_halfedges( Range &faces)
  {
    ErrorCode error;    
935 936
    EntityType ftype = mb->type_from_handle(*faces.begin());
    int nepf = lConnMap2D[ftype-2].num_verts_in_face;
937

nray's avatar
nray committed
938
    std::vector<char> markEdges(nepf*faces.size(), 0);
939

940
    for (Range::iterator it = faces.begin(); it != faces.end(); ++it){      
941
        EntityHandle fid = *it;
942
        const EntityHandle* conn;
943
        error = mb->get_connectivity(fid, conn, nepf, true);MB_CHK_ERR(error);
944

945 946
        for(int i=0; i<nepf; ++i){
            EntityHandle v = conn[i];
947
            int vidx = ID_FROM_HANDLE(v)-1;
948
            HFacet hf = v2he[vidx];
949 950 951

            if (hf == 0 && (v2hes.empty() || (v2hes.find(v) == v2hes.end())))
              {
952
                //This is the first time a half-facet is assigned to a vertex.
953
                HFacet nwhf=0;
954
                error = mark_halfedges(v, fid, i, faces, markEdges,nwhf);MB_CHK_ERR(error);
955 956 957 958 959 960

                if (nwhf ==0)
                  nwhf = create_halffacet(fid, i);

                v2he[vidx] = nwhf;
              }
961
            else if (hf != 0 && !markEdges[nepf*faces.index(fid)+i])
962 963 964 965
              {
                //This is the first time a non-manifold vertex is encountered. Copy the existing he in v2he[v] to the multimap.
                v2hes.insert(std::pair<EntityHandle,HFacet>(v,hf));
                HFacet nwhf = 0;
966
                error = mark_halfedges(v, fid, i, faces, markEdges, nwhf);MB_CHK_ERR(error);
967 968 969 970 971 972 973

                if (nwhf == 0)
                  nwhf = create_halffacet(fid, i);

                v2hes.insert(std::pair<EntityHandle,HFacet>(v, nwhf));
                v2he[vidx] = 0;
              }
974
            else if (hf == 0 && (!v2hes.empty()) && (v2hes.find(v) != v2hes.end()) && !markEdges[nepf*faces.index(fid)+i])
975
              {
976
                //This is check if reached if the vertex is non-manifold and has encountered a half-facet to a new component.
977
                HFacet nwhf = 0;
978
                error = mark_halfedges(v, fid, i, faces, markEdges, nwhf);MB_CHK_ERR(error);
979 980 981

                if (nwhf == 0)
                  nwhf = create_halffacet(fid, i);
982

983 984
                v2hes.insert(std::pair<EntityHandle,HFacet>(v, nwhf));
              }
985
          }
986 987
      }

988 989
   // error = print_tags(2);

990 991
    return MB_SUCCESS;
  }
992 993

  ///////////////////////////////////////////////////////////////////
nray's avatar
nray committed
994
  ErrorCode HalfFacetRep::mark_halfedges(EntityHandle vid, EntityHandle he_fid, int he_lid, Range &faces, std::vector<char> &markHEdgs, HFacet &bnd_hf)
995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013
  {
    ErrorCode error;
    EntityType ftype = mb->type_from_handle(he_fid);
    int nepf = lConnMap2D[ftype-2].num_verts_in_face;

    int qsize = 0, count = -1;
    int num_qvals = 0;

    error = gather_halfedges(vid, he_fid, he_lid, &qsize, &count);MB_CHK_ERR(error);

    while (num_qvals < qsize)
      {
        EntityHandle curfid = queue_fid[num_qvals];
        int curlid = queue_lid[num_qvals];
        num_qvals += 1;

        int fidx = ID_FROM_HANDLE(curfid)-1;

        const EntityHandle* conn;
1014
        error = mb->get_connectivity(curfid, conn, nepf, true);MB_CHK_ERR(error);
1015

nray's avatar
nray committed
1016 1017
        if (!markHEdgs[nepf*faces.index(curfid)+curlid] && (conn[curlid]==vid)){
            markHEdgs[nepf*faces.index(curfid)+curlid] = 1;
1018 1019 1020 1021 1022 1023 1024 1025 1026
            HFacet hf = sibhes[nepf*fidx+curlid];
            EntityHandle sibfid = fid_from_halfacet(hf, ftype);
            if (sibfid == 0)
              bnd_hf = create_halffacet(curfid, curlid);
          }

        EntityHandle he2_fid = 0; int he2_lid = 0;
        error = another_halfedge(vid, curfid, curlid, &he2_fid, &he2_lid);  MB_CHK_ERR(error);

nray's avatar
nray committed
1027
        if (!markHEdgs[nepf*faces.index(curfid)+he2_lid] && (conn[he2_lid]==vid)){
1028
            markHEdgs[nepf*faces.index(curfid)+he2_lid] = 1;
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
            HFacet hf = sibhes[nepf*fidx+he2_lid];
            EntityHandle sibfid = fid_from_halfacet(hf, ftype);
            if (sibfid == 0)
              bnd_hf = create_halffacet(he2_fid, he2_lid);
          }

        bool val = find_match_in_array(he2_fid, trackfaces, count);

        if (val)
          continue;

        count += 1;
        trackfaces[count] = he2_fid;

        error = get_up_adjacencies_2d(he2_fid, he2_lid, &qsize, &count);MB_CHK_ERR(error);
      }

    //Change the visited faces to false, also empty the queue
    for (int i = 0; i<=qsize; i++)
      {
        queue_fid[i] = 0;
        queue_lid[i] = 0;
      }

    for (int i = 0; i<=count; i++)
      trackfaces[i] = 0;
    return MB_SUCCESS;
  }

1058
  ///////////////////////////////////////////////////////////////////
1059
  ErrorCode HalfFacetRep::get_up_adjacencies_vert_2d(EntityHandle vid, std::vector<EntityHandle> &adjents)
1060 1061
  {
    ErrorCode error;
1062
    EntityType ftype = mb->type_from_handle(*_faces.begin());
nray's avatar
nray committed
1063

1064
    int vidx = ID_FROM_HANDLE(vid)-1;
1065
    HFacet hf = v2he[vidx];
1066

1067 1068
    std::vector<EntityHandle> start_fids;
    std::vector<int> start_lids;
1069

1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088
    if (hf == 0 && (v2hes.find(vid) != v2hes.end()))
      {
        std::pair <std::multimap<EntityHandle, HFacet>::iterator, std::multimap<EntityHandle, HFacet>::iterator> it_hes;
        it_hes = v2hes.equal_range(vid);

        for (std::multimap<EntityHandle, HFacet>::iterator it = it_hes.first; it != it_hes.second; ++it)
          {
            start_fids.push_back(fid_from_halfacet(it->second, ftype));
            start_lids.push_back(lid_from_halffacet(it->second));
          }
      }
    else if (hf != 0)
      {
        start_fids.push_back( fid_from_halfacet(hf, ftype));
        start_lids.push_back(lid_from_halffacet(hf));
      }

    if (start_fids.empty())
      return MB_SUCCESS;
1089

1090 1091
    int qsize = 0, count = -1;
    int num_qvals = 0;
1092

1093
    adjents.reserve((int)start_fids.size());
1094 1095 1096

    for (int i=0; i<(int)start_fids.size(); i++)
      {
1097 1098
        adjents.push_back(start_fids[i]);
        error = gather_halfedges(vid, start_fids[i], start_lids[i], &qsize, &count);MB_CHK_ERR(error);
1099
      }
1100

1101 1102 1103 1104 1105
    while (num_qvals < qsize)
      {
        EntityHandle curfid = queue_fid[num_qvals];
        int curlid = queue_lid[num_qvals];
        num_qvals += 1;
1106

1107
        EntityHandle he2_fid = 0; int he2_lid = 0;
1108
        error = another_halfedge(vid, curfid, curlid, &he2_fid, &he2_lid);MB_CHK_ERR(error);
1109

1110
        bool val = find_match_in_array(he2_fid, trackfaces, count);
1111

1112 1113
        if (val)
          continue;
1114

1115 1116
        count += 1;
        trackfaces[count] = he2_fid;
1117

nray's avatar
nray committed
1118
        error = get_up_adjacencies_2d(he2_fid, he2_lid, &qsize, &count);MB_CHK_ERR(error);
1119 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131

        adjents.push_back(he2_fid);
      }

    //Change the visited faces to false, also empty the queue
    for (int i = 0; i<=qsize; i++)
      {
        queue_fid[i] = 0;
        queue_lid[i] = 0;
      }

    for (int i = 0; i<=count; i++)
      trackfaces[i] = 0;
nray's avatar
nray committed
1132

1133
    return MB_SUCCESS;
1134
  }
1135

1136
  //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1137 1138 1139
  ErrorCode  HalfFacetRep::get_up_adjacencies_2d( EntityHandle eid,
                                                  std::vector<EntityHandle> &adjents,
                                                  std::vector<int>  * leids)
1140 1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151
{

  // Given an explicit edge eid, find the incident faces.
    ErrorCode error;
    EntityHandle he_fid=0; int he_lid=0;

    // Step 1: Given an explicit edge, find a corresponding half-edge from the surface mesh.
    bool found = find_matching_halfedge(eid, &he_fid, &he_lid);

    // Step 2: If there is a corresponding half-edge, collect all sibling half-edges and store the incident faces.
    if (found)
      { 
nray's avatar
nray committed
1152
        error = get_up_adjacencies_2d(he_fid, he_lid, true, adjents, leids);MB_CHK_ERR(error);
1153 1154 1155 1156 1157 1158
      }

    return MB_SUCCESS;
  }

  /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1159
  ErrorCode HalfFacetRep::get_up_adjacencies_2d(EntityHandle fid,
1160 1161
                                                 int leid,
                                                 bool add_inent,
1162 1163
                                                 std::vector<EntityHandle> &adj_ents,
                                                 std::vector<int> *adj_leids, std::vector<int> *adj_orients)
1164 1165 1166
  {
    // Given an implicit half-edge <fid, leid>, find the incident half-edges.
    ErrorCode error;
1167

1168 1169
    EntityType ftype = mb->type_from_handle(fid);
    int nepf = lConnMap2D[ftype-2].num_verts_in_face;
1170 1171 1172

    if (!fid)
      return MB_FAILURE;
1173
    adj_ents.reserve(20);
1174

1175 1176
    bool local_id =false;
    bool orient = false;
1177
    if (adj_leids != NULL)
1178
      {
1179
        local_id = true;
1180 1181
        adj_leids->reserve(20);
      }
1182
    if (adj_orients != NULL)
1183 1184 1185 1186
      {
        orient = true;
        adj_orients->reserve(20);
      }
1187

1188 1189
    if (add_inent)
      {
1190
        adj_ents.push_back(fid);
1191
        if (local_id)
1192
          adj_leids->push_back(leid);
1193 1194
      }

1195
    EntityHandle fedge[2] = {0,0};
1196 1197 1198 1199

    if (orient)
      {
        //get connectivity and match their directions
1200
        const EntityHandle* fid_conn;
1201
        error = mb->get_connectivity(fid, fid_conn, nepf, true);MB_CHK_ERR(error);
1202

1203
        int nidx = lConnMap2D[ftype-2].next[leid];
1204
        fedge[0] = fid_conn[leid];
1205
        fedge[1] = fid_conn[nidx];
1206 1207
      }

1208
    int fidx = ID_FROM_HANDLE(fid)-1;
1209 1210 1211
    HFacet hf = sibhes[nepf*fidx+leid];
    EntityHandle curfid = fid_from_halfacet(hf, ftype);
    int curlid = lid_from_halffacet(hf);
1212
    
1213
    while ((curfid != fid)&&(curfid != 0)){//Should not go into the loop when no sibling exists
1214
        adj_ents.push_back(curfid);
1215

1216
        if (local_id)
1217
          adj_leids->push_back(curlid);
1218

1219 1220 1221
        if (orient)
          {
            //get connectivity and match their directions
1222
            const EntityHandle* conn;
1223
            error = mb->get_connectivity(curfid, conn, nepf, true);MB_CHK_ERR(error);
1224

1225
            int nidx = lConnMap2D[ftype-2].next[curlid];
1226

1227
            if ((fedge[0] == conn[curlid])&&(fedge[1] == conn[nidx]))
1228
              adj_orients->push_back(1);
1229
            else if ((fedge[1] == conn[curlid])&&(fedge[0] == conn[nidx]))
1230 1231 1232
              adj_orients->push_back(0);
          }

1233
        int cidx = ID_FROM_HANDLE(curfid)-1;
1234 1235 1236
        hf = sibhes[nepf*cidx+curlid];